Adaptive AGN heating temperature jump

examples/EAGLE_ICs/EAGLE_100/eagle_100.yml

@@ -217,7 +217,16 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.1        # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               1          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  1          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_norm:              3e8        # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_reference:         1e8        # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_exponent:          0.666667.  # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1).
+  AGN_delta_T_crit_factor:            1.0        # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_background_factor:      0.0        # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_min:                    1e7        # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_max:                    3e9        # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                1e20       # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition).
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_ICs/EAGLE_12/eagle_12.yml

@@ -216,7 +216,16 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.1        # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               1          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  1          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_norm:              3e8        # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_reference:         1e8        # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_exponent:          0.666667.  # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1).
+  AGN_delta_T_crit_factor:            1.0        # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_background_factor:      0.0        # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_min:                    1e7        # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_max:                    3e9        # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                1e20       # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition).
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_ICs/EAGLE_25/eagle_25.yml

@@ -216,7 +216,16 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.1        # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               1          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  1          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_norm:              3e8        # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_reference:         1e8        # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_exponent:          0.666667.  # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1).
+  AGN_delta_T_crit_factor:            1.0        # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_background_factor:      0.0        # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_min:                    1e7        # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_max:                    3e9        # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                1e20       # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition).
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_ICs/EAGLE_50/eagle_50.yml

@@ -214,7 +214,16 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.1        # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               1          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  1          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_norm:              3e8        # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_reference:         1e8        # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_exponent:          0.666667.  # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1).
+  AGN_delta_T_crit_factor:            1.0        # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_background_factor:      0.0        # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_min:                    1e7        # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_max:                    3e9        # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                1e20       # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition).
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml

@@ -213,7 +213,16 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.1        # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               1          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  1          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_norm:              3e8        # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_reference:         1e8        # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_exponent:          0.666667.  # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1).
+  AGN_delta_T_crit_factor:            1.0        # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_background_factor:      0.0        # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_min:                    1e7        # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_max:                    3e9        # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                1e20       # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition).
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_ICs/EAGLE_6/eagle_6.yml

@@ -216,7 +216,16 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.1        # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               1          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  1          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_norm:              3e8        # Normalisation temperature of AGN dT scaling with BH subgrid mass [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_reference:         1e8        # BH subgrid mass at which the normalisation temperature set above applies [M_Sun] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_mass_exponent:          0.666667.  # Power-law index of AGN dT scaling with BH subgrid mass (only used if use_variable_delta_T is 1).
+  AGN_delta_T_crit_factor:            1.0        # Multiple of critical dT for numerical efficiency (Dalla Vecchia & Schaye 2012) to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_background_factor:      0.0        # Multiple of local gas temperature to use as dT floor (only used if use_variable_delta_T and AGN_with_locally_adaptive_delta_T are both 1).
+  AGN_delta_T_min:                    1e7        # Minimum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_max:                    3e9        # Maximum allowed value of AGN dT [K] (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                1e20       # Maximal BH mass considered for BH repositioning in solar masses (large number implies we always reposition).
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_ICs/README

@@ -73,7 +73,12 @@ the following changes have been made (at standard resolution):
  - The angular momentum term in the BH accretion model of
    Rosas-Guevara et al. (2015) is now switched off.
  - The coupling efficency of the BH feedback is now 0.1 (was 0.15).
-
+ - The AGN feedback temperature jump is now a function of the BH
+   subgrid mass (was a constant at 10^8.5K). The temperature varies
+   between 10^7 and 3*10^9 K with a power-law of the subgrid mass
+   exponent of 2/3 and a reference point of 10^8 Msun. The temperature
+   jump also has to be larger than the T_crit for efficiency given
+   by the Dalla Vecchia & Schaye 2012 criterion.
 
 The scripts in this directory download the tables required to
 run the EAGLE model. Plotting scripts are also provided

examples/EAGLE_low_z/EAGLE_12/eagle_12.yml

@@ -200,7 +200,9 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               0          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  0          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                2e8        # Maximal BH mass considered for BH repositioning in solar masses.
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_low_z/EAGLE_25/eagle_25.yml

@@ -208,7 +208,9 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               0          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  0          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                2e8        # Maximal BH mass considered for BH repositioning in solar masses.
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_low_z/EAGLE_50/eagle_50.yml

@@ -199,7 +199,9 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               0          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  0          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                2e8        # Maximal BH mass considered for BH repositioning in solar masses.
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/EAGLE_low_z/EAGLE_6/eagle_6.yml

@@ -211,7 +211,9 @@ EAGLEAGN:
   max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
   eddington_fraction_for_recording:   0.1        # Record the last time BHs reached an Eddington ratio above this threshold.
   coupling_efficiency:                0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
+  use_variable_delta_T:               0          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  0          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
   AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.
   max_reposition_mass:                2e8        # Maximal BH mass considered for BH repositioning in solar masses.
   max_reposition_distance_ratio:      3.0        # Maximal distance a BH can be repositioned, in units of the softening length.

examples/SubgridTests/BlackHoleSwallowing/swallowing.yml

@@ -94,10 +94,12 @@ EAGLECooling:
 
 # EAGLE AGN model
 EAGLEAGN:
-  subgrid_seed_mass_Msun:           1.5e5      # Black hole subgrid mass at creation time in solar masses.
-  max_eddington_fraction:           1.         # Maximal allowed accretion rate in units of the Eddington rate.
-  viscous_alpha:                    1e6        # Normalisation constant of the Bondi viscuous time-scale accretion reduction term.
-  radiative_efficiency:             0.1        # Fraction of the accreted mass that gets radiated.
-  coupling_efficiency:              0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
-  AGN_delta_T_K:                    3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
-  AGN_num_ngb_to_heat:              1.         # Target number of gas neighbours to heat in an AGN feedback event.
+  subgrid_seed_mass_Msun:             1.5e5      # Black hole subgrid mass at creation time in solar masses.
+  max_eddington_fraction:             1.         # Maximal allowed accretion rate in units of the Eddington rate.
+  viscous_alpha:                      1e6        # Normalisation constant of the Bondi viscuous time-scale accretion reduction term.
+  radiative_efficiency:               0.1        # Fraction of the accreted mass that gets radiated.
+  coupling_efficiency:                0.15       # Fraction of the radiated energy that couples to the gas in feedback events.
+  use_variable_delta_T:               0          # Switch to enable adaptive calculation of AGN dT [1], rather than using a constant value [0].
+  AGN_with_locally_adaptive_delta_T:  0          # Switch to enable additional dependence of AGN dT on local gas density and temperature (only used if use_variable_delta_T is 1).
+  AGN_delta_T_K:                      3.16228e8  # Change in temperature to apply to the gas particle in an AGN feedback event [K] (only used if use_variable_delta_T is 0).
+  AGN_num_ngb_to_heat:                1.         # Target number of gas neighbours to heat in an AGN feedback event.

src/black_holes/EAGLE/black_holes.h

@@ -460,6 +460,63 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_bpart(
   bpi->number_of_mergers++;
 }
 
+/**
+ * @brief Computes the temperature increase delta_T for black hole feedback.
+ *
+ * This is calculated as delta_T = min(max(dT_crit, T_gas), dT_num, dT_max):
+ * dT_crit = f_crit * critical temperature for suppressing numerical losses
+ * T_gas = f_gas * temperature of ambient gas
+ * dT_num = temperature increase affordable if N particles should be heated
+ * dT_max = maximum allowed energy increase.
+ *
+ * @param bp The #bpart.
+ * @param props The properties of the black hole model.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static double black_hole_feedback_delta_T(
+    const struct bpart* bp, const struct black_holes_props* props,
+    const struct cosmology* cosmo) {
+
+  /* If we do not want a variable delta T, we can stop right here. */
+  if (!props->use_variable_delta_T) return props->AGN_delta_T_desired;
+
+  if (bp->internal_energy_gas < 0)
+    error("Attempting to compute feedback energy for BH without neighbours.");
+
+  /* Black hole properties */
+  const double n_gas_phys = bp->rho_gas * cosmo->a3_inv * props->rho_to_n_cgs;
+  const double mean_ngb_mass = bp->ngb_mass / ((double)bp->num_ngbs);
+  const double T_gas = bp->internal_energy_gas *
+                       cosmo->a_factor_internal_energy /
+                       props->temp_to_u_factor;
+
+  /* Calculate base line delta T from BH subgrid mass. The assumption is that
+   * the BH mass scales (via halo mass) with the virial temperature, so that
+   * this aims for delta T > T_vir. */
+  double delta_T = props->AGN_delta_T_mass_norm *
+                   pow((bp->subgrid_mass / props->AGN_delta_T_mass_reference),
+                       props->AGN_delta_T_mass_exponent);
+
+  /* If desired, also make sure that delta T is not below the numerically
+   * critical temperature or that of the ambient gas */
+  if (props->AGN_with_locally_adaptive_delta_T) {
+
+    /* Critical temperature for numerical efficiency, based on equation 18
+     * of Dalla Vecchia & Schaye (2012) */
+    const double T_crit =
+        3.162e7 * pow(n_gas_phys * 0.1, 0.6666667) *
+        pow(mean_ngb_mass * props->mass_to_solar_mass * 1e-6, 0.33333333);
+    delta_T = max(delta_T, T_crit * props->AGN_delta_T_crit_factor);
+
+    /* Delta_T should be (at least) some multiple of the local gas T */
+    delta_T = max(delta_T, T_gas * props->AGN_delta_T_background_factor);
+  }
+
+  /* Respect the limits */
+  delta_T = max(delta_T, props->AGN_delta_T_min);
+  return min(delta_T, props->AGN_delta_T_max);
+}
+
 /**
  * @brief Compute the accretion rate of the black hole and all the quantites
  * required for the feedback loop.
@@ -498,9 +555,6 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
   const double epsilon_r = props->epsilon_r;
   const double epsilon_f = props->epsilon_f;
   const double num_ngbs_to_heat = props->num_ngbs_to_heat;
-  const double delta_T = props->AGN_delta_T_desired;
-  const double delta_u = delta_T * props->temp_to_u_factor;
-  const double alpha_visc = props->alpha_visc;
   const int with_angmom_limiter = props->with_angmom_limiter;
 
   /* (Subgrid) mass of the BH (internal units) */
@@ -644,8 +698,9 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
     const double r_times_v_tang = Bondi_radius * tangential_velocity;
     const double r_times_v_tang_3 =
         r_times_v_tang * r_times_v_tang * r_times_v_tang;
-    const double viscous_time = 2. * M_PI * r_times_v_tang_3 /
-                                (1e-6 * alpha_visc * G * G * BH_mass * BH_mass);
+    const double viscous_time =
+        2. * M_PI * r_times_v_tang_3 /
+        (1e-6 * props->alpha_visc * G * G * BH_mass * BH_mass);
 
     const double f_visc = min(Bondi_time / viscous_time, 1.);
     bp->f_visc = f_visc;
@@ -706,6 +761,11 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
   bp->accreted_angular_momentum[2] +=
       bp->circular_velocity_gas[2] * mass_rate * dt / bp->h;
 
+  /* Now find the temperature increase for a possible feedback event */
+  const double delta_T = black_hole_feedback_delta_T(bp, props, cosmo);
+  bp->AGN_delta_T = delta_T;
+  const double delta_u = delta_T * props->temp_to_u_factor;
+
   /* Energy required to have a feedback event
    * Note that we have subtracted the particles we swallowed from the ngb_mass
    * and num_ngbs accumulators. */

src/black_holes/EAGLE/black_holes_io.h

@@ -139,7 +139,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
                                                int with_cosmology) {
 
   /* Say how much we want to write */
-  *num_fields = 32;
+  *num_fields = 33;
 
   /* List what we want to write */
   list[0] = io_make_output_field_convert_bpart(
@@ -347,6 +347,12 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "Integer number of gas neighbour particles within the black hole "
       "kernels.");
 
+  list[32] = io_make_output_field(
+      "FeedbackDeltaT", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, bparts,
+      AGN_delta_T,
+      "Temperature by which gas particles have been heated by the black hole "
+      "particles in the most recent feedback event.");
+
 #ifdef DEBUG_INTERACTIONS_BLACK_HOLES
 
   list += *num_fields;

src/black_holes/EAGLE/black_holes_part.h

@@ -160,6 +160,9 @@ struct bpart {
   /*! Total (physical) angular momentum accumulated from subgrid accretion */
   float accreted_angular_momentum[3];
 
+  /*! Instantaneous temperature increase for feedback */
+  float AGN_delta_T;
+
   /*! Union for the last high Eddington ratio point in time */
   union {
 

src/black_holes/EAGLE/black_holes_properties.h

@@ -93,9 +93,35 @@ struct black_holes_props {
   /*! Feedback coupling efficiency of the black holes. */
   float epsilon_f;
 
-  /*! Temperature increase induced by AGN feedback (Kelvin) */
+  /*! (Constant) temperature increase induced by AGN feedback [Kelvin] */
   float AGN_delta_T_desired;
 
+  /*! Switch on adaptive heating temperature scheme? */
+  int use_variable_delta_T;
+
+  /*! If we use variable delta_T, should we scale with local gas properties
+   *  in addition to BH mass? */
+  int AGN_with_locally_adaptive_delta_T;
+
+  /*! Normalisation for dT scaling with BH mass */
+  float AGN_delta_T_mass_norm;
+
+  /*! Reference BH mass for dT scaling [M_Sun] */
+  float AGN_delta_T_mass_reference;
+
+  /*! Exponent for dT scaling with BH mass */
+  float AGN_delta_T_mass_exponent;
+
+  /*! Buffer factor for numerical efficiency temperature */
+  float AGN_delta_T_crit_factor;
+
+  /*! Buffer factor for background temperature */
+  float AGN_delta_T_background_factor;
+
+  /*! Max/min temperature increase induced by AGN feedback [Kelvin] */
+  float AGN_delta_T_max;
+  float AGN_delta_T_min;
+
   /*! Number of gas neighbours to heat in a feedback event */
   float num_ngbs_to_heat;
 
@@ -150,6 +176,12 @@ struct black_holes_props {
 
   /*! Conversion factor from temperature to internal energy */
   float temp_to_u_factor;
+
+  /*! Conversion factor from physical density to n_H [cgs] */
+  float rho_to_n_cgs;
+
+  /*! Conversion factor from internal mass to solar masses */
+  float mass_to_solar_mass;
 };
 
 /**
@@ -257,8 +289,38 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
   bp->epsilon_f =
       parser_get_param_float(params, "EAGLEAGN:coupling_efficiency");
 
-  bp->AGN_delta_T_desired =
-      parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_K");
+  const double T_K_to_int =
+      1. / units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+  bp->use_variable_delta_T =
+      parser_get_param_int(params, "EAGLEAGN:use_variable_delta_T");
+  if (bp->use_variable_delta_T) {
+    bp->AGN_delta_T_mass_norm =
+        parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_mass_norm") *
+        T_K_to_int;
+    bp->AGN_delta_T_mass_reference =
+        parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_mass_reference") *
+        phys_const->const_solar_mass;
+    bp->AGN_delta_T_mass_exponent =
+        parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_mass_exponent");
+
+    bp->AGN_with_locally_adaptive_delta_T = parser_get_param_int(
+        params, "EAGLEAGN:AGN_with_locally_adaptive_delta_T");
+    if (bp->AGN_with_locally_adaptive_delta_T) {
+      bp->AGN_delta_T_crit_factor =
+          parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_crit_factor");
+      bp->AGN_delta_T_background_factor = parser_get_param_float(
+          params, "EAGLEAGN:AGN_delta_T_background_factor");
+    }
+
+    bp->AGN_delta_T_max =
+        parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_max") * T_K_to_int;
+    bp->AGN_delta_T_min =
+        parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_min") * T_K_to_int;
+  } else {
+    bp->AGN_delta_T_desired =
+        parser_get_param_float(params, "EAGLEAGN:AGN_delta_T_K");
+  }
 
   bp->num_ngbs_to_heat =
       parser_get_param_float(params, "EAGLEAGN:AGN_num_ngb_to_heat");
@@ -345,6 +407,15 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
   const double m_p = phys_const->const_proton_mass;
   const double mu = hydro_props->mu_ionised;
   bp->temp_to_u_factor = k_B / (mu * hydro_gamma_minus_one * m_p);
+
+  /* Calculate conversion factor from rho to n_H.
+   * Note this assumes primoridal abundance */
+  const double X_H = hydro_props->hydrogen_mass_fraction;
+  bp->rho_to_n_cgs =
+      (X_H / m_p) * units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY);
+
+  /* Conversion factor for internal mass to M_solar */
+  bp->mass_to_solar_mass = 1. / phys_const->const_solar_mass;
 }
 
 /**