Implement the reduction of the accretion rate of the BHs based on the local...

Implement the reduction of the accretion rate of the BHs based on the local viscuous time-scale in the EAGLE model.

examples/EAGLE_ICs/EAGLE_12/eagle_12.yml

@@ -165,6 +165,7 @@ EAGLEFeedback:
 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.
+  viscuous_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.

examples/EAGLE_ICs/EAGLE_25/eagle_25.yml

@@ -166,6 +166,7 @@ EAGLEFeedback:
 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.
+  viscuous_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.

examples/EAGLE_ICs/EAGLE_50/eagle_50.yml

@@ -166,6 +166,7 @@ EAGLEFeedback:
 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.
+  viscuous_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.

examples/EAGLE_low_z/EAGLE_12/eagle_12.yml

@@ -157,6 +157,7 @@ EAGLEFeedback:
 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.
+  viscuous_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.

examples/EAGLE_low_z/EAGLE_25/eagle_25.yml

@@ -164,6 +164,7 @@ EAGLEFeedback:
 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.
+  viscuous_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.

examples/EAGLE_low_z/EAGLE_50/eagle_50.yml

@@ -156,6 +156,7 @@ EAGLEFeedback:
 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.
+  viscuous_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.

examples/EAGLE_low_z/EAGLE_6/eagle_6.yml

@@ -166,6 +166,7 @@ EAGLEFeedback:
 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.
+  viscuous_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.

examples/parameter_example.yml

@@ -396,6 +396,7 @@ EAGLEFeedback:
 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.
+  viscuous_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.

src/black_holes/EAGLE/black_holes.h

@@ -82,6 +82,9 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart(
   bp->velocity_gas[0] = 0.f;
   bp->velocity_gas[1] = 0.f;
   bp->velocity_gas[2] = 0.f;
+  bp->circular_velocity_gas[0] = 0.f;
+  bp->circular_velocity_gas[1] = 0.f;
+  bp->circular_velocity_gas[2] = 0.f;
   bp->ngb_mass = 0.f;
   bp->num_ngbs = 0;
 }
@@ -136,6 +139,9 @@ __attribute__((always_inline)) INLINE static void black_holes_end_density(
   bp->velocity_gas[0] *= h_inv_dim;
   bp->velocity_gas[1] *= h_inv_dim;
   bp->velocity_gas[2] *= h_inv_dim;
+  bp->circular_velocity_gas[0] *= h_inv_dim;
+  bp->circular_velocity_gas[1] *= h_inv_dim;
+  bp->circular_velocity_gas[2] *= h_inv_dim;
 
   const float rho_inv = 1.f / bp->rho_gas;
 
@@ -144,6 +150,9 @@ __attribute__((always_inline)) INLINE static void black_holes_end_density(
   bp->velocity_gas[0] *= rho_inv;
   bp->velocity_gas[1] *= rho_inv;
   bp->velocity_gas[2] *= rho_inv;
+  bp->circular_velocity_gas[0] *= rho_inv * h_inv;
+  bp->circular_velocity_gas[1] *= rho_inv * h_inv;
+  bp->circular_velocity_gas[2] *= rho_inv * h_inv;
 }
 
 /**
@@ -290,12 +299,13 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
   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;
 
   /* (Subgrid) mass of the BH (internal units) */
   const double BH_mass = bp->subgrid_mass;
 
   /* Convert the quantities we gathered to physical frame (all internal units)
-   */
+   * Note: for the velocities this means peculiar velocities */
   const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
   const double gas_c_phys = bp->sound_speed_gas * cosmo->a_factor_sound_speed;
   const double gas_v_peculiar[3] = {bp->velocity_gas[0] * cosmo->a_inv,
@@ -306,6 +316,11 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
                                    bp->v[1] * cosmo->a_inv,
                                    bp->v[2] * cosmo->a_inv};
 
+  const double gas_v_circular[3] = {
+      bp->circular_velocity_gas[0] * cosmo->a_inv,
+      bp->circular_velocity_gas[1] * cosmo->a_inv,
+      bp->circular_velocity_gas[2] * cosmo->a_inv};
+
   /* Difference in peculiar velocity between the gas and the BH
    * Note that there is no need for a Hubble flow term here. We are
    * computing the gas velocity at the position of the black hole. */
@@ -316,13 +331,31 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
                               v_diff_peculiar[1] * v_diff_peculiar[1] +
                               v_diff_peculiar[2] * v_diff_peculiar[2];
 
+  /* Norm of the cirecular velocity of the gas around the BH */
+  const double tangential_velocity2 = gas_v_circular[0] * gas_v_circular[0] +
+                                      gas_v_circular[0] * gas_v_circular[0] +
+                                      gas_v_circular[0] * gas_v_circular[0];
+  const double tangential_velocity = sqrt(tangential_velocity2);
+
   /* We can now compute the Bondi accretion rate (internal units) */
   const double gas_c_phys2 = gas_c_phys * gas_c_phys;
   const double denominator2 = v_diff_norm2 + gas_c_phys2;
   const double denominator_inv = 1. / sqrt(denominator2);
-  const double Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass *
-                            gas_rho_phys * denominator_inv * denominator_inv *
-                            denominator_inv;
+  double Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass * gas_rho_phys *
+                      denominator_inv * denominator_inv * denominator_inv;
+
+  /* Compute the reduction factor from Rosas-Guevara et al. (2015) */
+  const double Bondi_radius = G * BH_mass / gas_c_phys2;
+  const double Bondi_time = Bondi_radius / gas_c_phys;
+  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 f_visc = max(Bondi_time / viscous_time, 1.);
+
+  /* Limit the Bondi rate by the Bondi viscuous time ratio */
+  Bondi_rate *= f_visc;
 
   /* Compute the Eddington rate (internal units) */
   const double Eddington_rate =

src/black_holes/EAGLE/black_holes_iact.h

@@ -84,6 +84,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
   bi->velocity_gas[1] += mj * vj[1] * wi;
   bi->velocity_gas[2] += mj * vj[2] * wi;
 
+  /* Contribution to the circular valocity */
+  const float dv[3] = {bi->v[0] - vj[0], bi->v[1] - vj[1], bi->v[2] - vj[2]};
+  bi->circular_velocity_gas[0] += mj * wi * (dx[1] * dv[2] - dx[2] * dv[1]);
+  bi->circular_velocity_gas[1] += mj * wi * (dx[0] * dv[2] - dx[2] * dv[0]);
+  bi->circular_velocity_gas[2] += mj * wi * (dx[0] * dv[1] - dx[1] * dv[0]);
+
 #ifdef DEBUG_INTERACTIONS_BH
   /* Update ngb counters */
   if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_BH)

src/black_holes/EAGLE/black_holes_part.h

@@ -94,6 +94,9 @@ struct bpart {
   /*! Smoothed velocity (peculiar) of the gas surrounding the black hole */
   float velocity_gas[3];
 
+  /*! Curl of the velocity field around the black hole */
+  float circular_velocity_gas[3];
+
   /*! Total mass of the gas neighbours. */
   float ngb_mass;
 

src/black_holes/EAGLE/black_holes_properties.h

@@ -63,6 +63,9 @@ struct black_holes_props {
   /*! Feedback coupling efficiency of the black holes. */
   float epsilon_f;
 
+  /*! Normalisation of the viscuous angular momentum accretion reduction */
+  float alpha_visc;
+
   /* ---- Properties of the feedback model ------- */
 
   /*! Temperature increase induced by AGN feedback (Kelvin) */
@@ -144,6 +147,7 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
       parser_get_param_float(params, "EAGLEAGN:radiative_efficiency");
   bp->epsilon_f =
       parser_get_param_float(params, "EAGLEAGN:coupling_efficiency");
+  bp->alpha_visc = parser_get_param_float(params, "EAGLEAGN:viscuous_alpha");
 
   /* Feedback parameters ---------------------------------- */