BLACK_HOLES

The module for black holes includes methods for BH seeding/merging, movement/positioning, accretion, and feedback.

The majority of these aspects were originally developed for the “Illustris model”. Subsequent additions (particularly to the feedback model) were developed for the “Auriga model” as well as (particularly to the low-state feedback model) for the “IllustrisTNG model”.

Usage

The BLACK_HOLES master switch enables the module as a whole. It is possible to use only specific subsets of the functionality described herein.

Note

In general, BLACK_HOLES is independent of GFM, and the two do not need to be active together. The only exception is GFM_AGN_RADIATION.

Warning

As the number of configuration options is large, not all possible permutations and combinations have been tested. There are also a significant number of non-functional and/or depreciated pieces of code within BLACK_HOLES.

Additional Parameters

Seeding and accretion:

BlackHoleAccretionFactor the multiplicative “boost” parameter which modifies the otherwise default Bondi accretion rate, i.e. the variable \(\alpha\) of Eqn. 1 of Sijacki+ 2007. Illustris value: 100. TNG value: 1.
BlackHoleEddingtonFactor fraction, or factor, of the Eddington rate to enforce as the maximum accretion rate. If unity, then black hole accretion is Eddington limited. Only allowed if BH_BONDI_DENSITY is not set. Illustris/TNG value: 1.
SeedBlackHoleMass the initial seed mass that newly created black holes are given. Code mass units. Illustris value: 1e-5. TNG value: 8.0e-5.
MinFoFMassForNewSeed the friends-of-friends halo mass threshold for on the fly black hole seeding. If the FoF halo finder identifies a halo above this minimum mass which does not already have a black hole, then a central particle/cell of that halo is converted into a black hole. Code mass units. Illustris/TNG value: 5.0.
DesNumNgbBlackHole the number of (weighted) gas cell neighbors, used for accretion and possibly also feedback. That is, \(n_{\rm ngb}\) of Eqn. 1 of Weinberger+ 2017. It therefore defines the set of gas cells from which local gas properties are derived, such as density, pressure, and/or velocity. This value is used to set BH_Hsml for each black hole particle, which is then also used to identify mergers, as well as e.g. the region for feedback energy injection (depending on model). As a function of numerical resolution, if one wants e.g. a constant volume, this parameter would need to scale linearly with the total number of particles. In TNG we adopted a sub-linear scaling. TNG100-1 value: 256, TNG100-2 value: 128, TNG100-3 value: 64.
BlackHoleMaxAccretionRadius sets the maximum allowed value for BH_Hsml, therefore potentially affecting not only the accretion radius, but also e.g. the feedback injection region, as noted above. However, actually using this parameter (setting it to a small value) is not recommended, until the implementation improves. At present, if a sufficient number of neighboring cells cannot be found within the restricted radius, the code will terminate. TNG value: 1e20.

Feedback:

BlackHoleFeedbackFactor this is the main efficiency parameter for the high state, i.e. \(\epsilon_{\rm f,high}\) of Eqn. 7 of Weinberger+ 2017. Physically, this is the fraction of the radiatively released energy which couples to the surrounding gas. TNG value: 0.1, Illustris value: 0.05.
BlackHoleRadiativeEfficiency this is the \(\epsilon_{\rm r}\) factor, the radiatively efficiency applied to the SMBH mass accretion rate in order to determine the available energy, hypothetically released in the form of radiation. See Eqn. 7 of Weinberger+ 2017. Note that this same parameter is one factor to determine the value of \(\epsilon_{f,kin}\) of Eqn. 9 of Weinberger+ 2017, i.e. \(\epsilon_{f,kin} =\) BlackHoleRadiativeEfficiency \(\times\) RadioFeedbackFactor (with potential additional factors depending on Config-flags). TNG value: 0.2.

If the ADIOS low-state model is enabled:

QuasarThreshold is required only if BH_ADIOS_WIND_WITH_QUASARTHRESHOLD is enabled. This is \(\chi_0\) of Eqn. 5 of Weinberger+ 2017. Note that the corresponding value of \(\beta = 2\) from the same equation is hard-coded, and is not a parameter. TNG value: 0.002
RadioFeedbackFactor is required only if BH_ADIOS_WIND or BH_BUBBLES is enabled. In the case of BH_ADIOS_WIND, this factor multiplies the available feedback energy in the low state, i.e. it is part of the efficiency factor in Eqn. 8 of Weinberger+ 2017 with \(\epsilon_{f,kin} =\) BlackHoleRadiativeEfficiency \(\times\) RadioFeedbackFactor. Note that this parameter also affects the mode choice, if and only if BH_ADIOS_WIND_WITH_QUASARTHRESHOLD is not enabled. This parameter also affects the AGN luminosity estimate, if and only if it is less than 1.0, if GFM_AGN_RADIATION is enabled. This parameter also affects the energy division between modes if and only if BH_CONTINOUS_MODE_SWITCH is enabled. TNG value: 1.0, Illustris value: 0.35.
RadioFeedbackReiorientationFactor is required only if BH_ADIOS_RANDOMIZED is enabled. Controls the frequency (i.e. accumulation time) of kinetic wind energy injection events. This is \(f_{\rm re}\) of Eqn. 13 of Weinberger+ 2017. TNG value: 20
RadioFeedbackMinDensityFactor is required only if BH_ADIOS_ONLY_ABOVE_MINIMUM_DENSITY is enabled. This corresponds to the \(f_{\rm thresh}\) parameter from Eqn. 9 of Weinberger+ 2017. In the implementation, however, the equation is \(\epsilon_{k,\rm kin} =\) BlackHoleRadiativeEfficiency \(\times\) min(\(\rho_{\rm ngb}\) / (RadioFeedbackMinDensityFactor \(\times\) PhysDensThresh)), i.e., it has a pre-factor of BlackHoleRadiativeEfficiency which was pulled into the minimum in Weinberger+ 2017. TNG value: 0.01

Only if BH_ADIOS_DENS_DEP_EFFICIANCY is enabled, then the following four parameters set the functional form of the \(\epsilon_{\rm r}(n)\) scaling (unclear if this was an experimental implementation, i.e. never used in production, or if the functional form in the code reflects a published equation/paper) and its maximum value:

RadioFeedbackFactorPivotMass the parameter \(M_{\rm BH,pivot}\) in the scaling equation.
RadioFeedbackFactorSlope the parameter \(\gamma\) in the scaling equation.
RadioFeedbackFactorRefDensityFactor the parameter \(f_{\rm dens}\)
RadioFeedbackFactorMaxEfficiency the maximum value allowed for \(\epsilon_{\rm r}(n)\).

If the bubble low-state (radio mode) is enabled through BH_BUBBLES, then:

BubbleDistance normalization (in kpc/h) of the scaling-relation for the distance of the AGN-driven bubbles. This is Eqn. 5 of Sijacki+ 2007. Illustris value: 200.
BubbleRadius normalization (in kpc/h) of the scaling-relation for the radius of the AGN-driven bubbles. This is Eqn. 5 of Sijacki+ 2007. Illustris value: 50.
BubbleEnergy reference energy (erg) in the scaling-relations for the distance and radius of the AGN-driven bubbles. This is Eqn. 5 of Sijacki+ 2007. Illustris value: 1e60.
BlackHoleRadioTriggeringFactor the duty cycle of the radio-mode (burstyness), \(\delta_{\rm BH}\). Illustris value: 1.15.
DefaultICMDensity reference gas density (\(10^{10} \rm{M}_\odot h^2 / \rm{kpc}^3\)) in the scaling relations for the distance and radius of the AGN-driven bubbles (Eqn. 5 of Sijacki+ 2007). Illustris value: 1e-6.

If BH_NEW_CENTERING is enabled:

BlackHoleCenteringMassMultiplier the number of neighbors (within plus/minus five percent) which are considered to find the local gravitational potential minimum, to which the black hole is subsequently moved. TNG value: 1000.

If the GFM_AGN_RADIATION option is enabled:

TreecoolFileAGN path to a Treecool-like text file, giving the AGN radiation field heating rates as a function of AGN bolometric intensity. Each row of the file should have seven float columns. The first is the AGN bolometric intensity (same units as the PartType0/GFM_AGNRadiation snapshot field). The last six are the same as in the Treecool file: Gamma_HI, Gamma_HeI, Gamma_HeII, Qdot_HI, Qdot_HeI, Qdot_HeII, where ‘Gamma’ values are the photoionization rate (in units of \(s^{-1}\)) and ‘Qdot’ values are the photoheating rates (in units of \(\rm{erg s^{-1}}\)).
SelfShieldingDensity the hydrogen number density limit for shielding. Gas above this density is not affected by the AGN radiation field. Illustris/TNG value: 0.1295.
ObscurationFactor the normalization factor \(\omega_1\) for the strength of the bolometric luminosity dependent obscuration model (Eqn. 28 of Vogelsberger+ 2013.) Illustris/TNG value: 0.3.
ObscurationSlope the exponent factor \(\omega_2\) for the scaling of the bolometric luminosity dependent obscuration model (Eqn. 28 of Vogelsberger+ 2013.) Illustris/TNG value: 0.07.

If BH_SPIN_EVOLUTION is enabled:

BHInitialSpin is the initial spin of black holes. The recommended default is 0.0.
ShakuraSunyaevParameter is the \(\alpha\) parameter of the disk. The recommended default is 0.1.

Misc parameters:

DesNumNgbBlackHoleRelDeviationFactor is a factor which, multiplied by DesNumNgbBlackHole, sets the maximum allowed deviation on the number of black hole neighbors. Only if BH_RELATIVE_NGB_DEVIATION is enabled.
BlackholeDeltaTemp only required if BH_THERMALFEEDBACK_ACC is enabled.
BlackholeDeltaTime only required if BH_THERMALFEEDBACK_ACC is enabled.
MagneticEnergyFraction only required if BH_MAGNETIC_BUBBLES is enabled.
RadioThreshold only required if UNIFIED_FEEDBACK is enabled. Gives the threshold, in units of the black hole mass accretion rate normalized by the Eddington rate, below which the feedback mode switches from thermal to radio/bubble.
BHFrictionCoefficient only required if BH_FRICTION is enabled. Coefficient which multiplies the relative velocity over the dynamical time in order to compute the friction acceleration. Here, relative velocity is the motion of the black hole with respect to the minimum potential particle.
BHFrictionAvgTime only required if BH_FRICTION is enabled. Averaging time, in units of the cosmological scalefactor, over which to accumulate i.e. average the minimum potential particle velocity.

Additional Config.sh Options

BLACK_HOLES: The master switch for enabling black holes.

Seeding and Accretion

DRAINGAS=1: Non-stochastic smooth accretion. Options are: 1 = enabled, and drain gas from the single parent gas cell only. 2 = enabled, drain from single parent cell, using the gas density of that cell instead of the neighbor region BH_Density value in the accretion formula. 3 = enabled, drain gas from all cells within the (weighted) neighbor region.
BH_EXACT_INTEGRATION: If on, the process of mass accretion will be integrated analytically. Conversely, if not activated, accreted mass will be determined through the calculation of mdot multiplied by dt. Note: do not use this switch in conjunction with BH_AD_SPIN when using the CAMBS models.
BH_BONDI_DEFAULT: Default Bondi accretion prescription. Do not use this switch in conjunction with any of the three alternative Bondi options listed below. As described in Eqn. 1 of Sijacki+ 2007.
BH_BONDI_DENSITY: Density dependent Bondi accretion prescription. Do not use this switch in conjunction with BH_BONDI_DEFAULT, BH_BONDI_CAPTURE, BH_BONDI_DISK_VORTICITY. Note: also cannot be used together with BH_PRESSURE_CRITERION.
BH_BONDI_CAPTURE: Mass dependent Bondi accretion prescription. Do not use this switch in conjunction with BH_BONDI_DEFAULT, BH_BONDI_DENSITY, BH_BONDI_DISK_VORTICITY.
BH_BONDI_DISK_VORTICITY: Vorticity dependent Bondi accretion prescription. Do not use this switch in conjunction with BH_BONDI_DEFAULT, BH_BONDI_DENSITY, BH_BONDI_CAPTURE.
BH_DO_NOT_PREVENT_MERGERS: If enabled, there will be no limitations imposed on the merger of two black holes in terms of their relative velocity. It is imperative to activate this switch, as the code may not function properly if the velocities of particles in the tree are not imported.
BH_USE_GASVEL_IN_BONDI: When enabled, the velocity of the surrounding gas will be taken into consideration in addition to the sound speed when determining the Bondi accretion rate.
BH_USE_ALFVEN_SPEED_IN_BONDI: Take into account the alfven speed of the plasma in the BH accretion rate calculation. Also, the total gas pressure around the BH includes the magnetic contribution when compared the reference pressure in BH_PRESSURE_CRITERION (requires MHD).
MASSIVE_SEEDS: When a BH is seeded, acquire mass from a predefined number of gas neighbors (DesNumNgbSeed parameter) out to a maximum distance (SeedMaxAccretionRadaius parameter), which is then transfered to the BH seed dynamical mass. Given a sufficiently high mass resolution and a sufficiently massive BH seed mass, this option would ideally remove the need for repositioning, while mass conservation is ensured.
MASSIVE_SEEDS_MERGER: If enabled, this switch modifies the criteria for constraining the relative velocity of two black hole merger candidates. If BH_DO_NOT_PREVENT_MERGERS is activated, this option does not change the relative velocity criterion. However, this option always changes the search radius for considering BH-BH mergers, from BH_Hsml to five percent of the (approximate) virial radius of the host halo as determined by the on-the-fly FoF.

Positioning

BH_NEW_CENTERING: An alternative to the BH_FRICTION and REPOSITION_ON_POTMIN switches. This is still a repositioning scheme, to the local gravitational potential minimum. The parameter BlackHoleCenteringMassMultiplier is required. TNG value: enabled.
BH_FRICTION: Estimates the local DM density around BH, and applies a friction force to the relative velocity, meant as a replacement for REPOSITION_ON_POTMIN. Experimental, never used in production simulations or described in a published paper. In this case, two additional parameters are required: BHFrictionCoefficient and BHFrictionAvgTime.
BH_FRICTION_AGGRESSIVE: Requires BH_FRICTION. If enabled, a second minimum potential search is carried out, on a larger scale than the primary search, and if this minimum is smaller, it is used instead.
BH_HARMONIC_OSCILLATOR_FORCE: Requires BH_FRICTION. If enabled, the gravitational acceleration of black holes is set to \(-(4 \pi / 3) G \rho (\rm{Pos} - \rm{PosMinPot})\) in the sense of a harmonic oscillator.
BH_DRAG: Include drag on black-holes due to accretion. Note: The current implementation simply double-accounts the accretion momentum transfer and, therefore, should not be used. Cannot be used with BH_FRICTION.
REPOSITION_ON_POTMIN: Repositions the black hole on potential minimum. Requires EVALPOTENTIAL. Appears depreciated in practice, in favor of BH_NEW_CENTERING.
BH_REPOSITION_POTMIN_TRUST_THRESHOLD: If not specified, defaults to a value of 0.9. Only relevant if BH_FRICTION or REPOSITION_ON_POTMIN is enabled. In both cases, potential minima are only considered if the total active mass in the black hole region is this factor times the total mass in the same region, or larger.

Misc Options

BH_PRESSURE_CRITERION: If enabled, modifies the accretion rates of black holes in the quasar mode (only). If the total pressure (thermal, plus magnetic if BH_USE_ALFVEN_SPEED_IN_BONDI is enabled) is below a reference pressure, then Mdot is multiplied by this pressure ratio squared. This switch is described in Section 2.6.2 and Eqn. 23 of Vogelsberger+ 2013.
BH_RELATIVE_NGB_DEVIATION: If enabled, the maximum NGB number deviation will be calculated relative to the total number of neighbours. In this case the DesNumNgbBlackHoleRelDeviationFactor parameter is required.
OUTPUT_BLACK_HOLE_TIMESTEP: If enabled, output the 3 time-steps for black hole particles in the snapshot. These are the timestep requirements from the gravitational acceleration, gas accretion rate, and hydro/Courant.
BH_RECOIL_KICK: Applies a recoil kick to the black hole remnant just after a merger event. If BH_SPIN_EVOLUTION is also enabled, then the contribution to the kick from black hole spin is included, by kicking along the angular momentum direction, using the spin components on the orbital plane. Otherwise, the contribution from the mass asymmetry of the two black holes sets the kick. Experimental/never used in production.
BH_INFLOW_RATE: Attempt to measure gas inflow rate around black holes. Old/experimental code, will not actually be executed.

Thermal/Quasar Feedback

Feedback routines which inject thermal energy, either locally around the SMBH (i.e. the thermal/quasar mode of Illustris/TNG), as well as non-locally in inflated bubbles (i.e. the radio mode of Illustris/Auriga/Fable).

BH_THERMALFEEDBACK: Quasar-mode black hole feedback. It couples a fraction of the BH luminosity into the surrounding gas cells. Illustris/TNG/Auriga: enabled.
BH_THERMALFEEDBACK_ACC: Bursty quasar-mode black hole feedback that accumulate thermal energy. If enabled, then the two parameters BlackholeDeltaTemp and BlackholeDeltaTime are required. Illustris/TNG value: disabled.
BH_NF_RADIO: Radio-mode model based on the model of Nulsen&Fabian 2000. Broadly speaking, this injects enough bubble energy to directly offset the cooling of the hot halo gas. It does so with a large number of small and relatively gentle bubbles. Auriga: enabled.
UNIFIED_FEEDBACK: This is an experimental switch and is not recommended/used. It activates thermal mode feedback at high Mdot and radio/bubble mode feedback at low Mdot, the split determined by the ratio of Mdot to the Eddington accretion rate, compared to a threshold value given by the RadioThreshold parameter.

Radio/Bubble Feedback

BH_BUBBLES: The “bubble” low-state (radio mode) feedback model. It calculates the bubble energy directly from the black hole accretion rate. Illustris value: enabled.
BH_MAGNETIC_BUBBLES: If enabled, inject part of the bubble energy as magnetic energy. Requires MHD. If enabled the MagneticEnergyFraction parameter is required. Never used in production, not recommended.
BH_MAGNETIC_DIPOLAR_BUBBLES: If enabled, inject part of the bubble energy as magnetic energy, arranging the field as a dipole with random orientation. Never used in production, not recommended.
BH_NEW_RADIO: Unclear what this is, or was. In the past the documentation has written that both UNIFIED_FEEDBACK and BH_BUBBLES are obsolete and replaced by BH_NEW_RADIO, but this option does not exist.

Kinetic Wind Feedback

The “kinetic wind” model of TNG, as originally presented in Weinberger+ 2017. Designed to operate at low accretion rates. Gives kinetic energy kicks to gas cells near to the SMBH. Not compatible with several other feedback options: do NOT use together with BH_BUBBLES, BH_MAGNETIC_BUBBLES, BH_MAGNETIC_DIPOLAR_BUBBLES, UNIFIED_FEEDBACK, BH_NF_RADIO.

Note

This model can produce large, low density regions around black holes, in particular in cosmological simulations with low resolution. It can also create quite fast velocities, leading to small timesteps and a large timebin hierarchy.

BH_ADIOS_WIND: Master switch for kinetic feedback. Without BH_ADIOS_WIND_DIRECTIONAL or BH_ADIOS_RANDOMIZED, this will kick the gas cells within a given radius around the black hole radially outward in a momentum conserving way. TNG setting: enabled.
BH_ADIOS_WIND_DIRECTIONAL: Introduces a preferred direction for the wind, randomly chosen at every injection event. In this case, the weights which are applied to nearby gas cells are modified, from their usual radial distance only dependence, to a dependence which also includes the relative angle between the chosen wind direction and the cell direction relative to the SMBH. The weight is modified by \(cos^2(dx*dir_x + dy*dir_y + dz*dir_z)\). Requires BH_ADIOS_WIND. Do not use together with BH_ADIOS_RANDOMIZED.
BH_ADIOS_RANDOMIZED: Makes the wind kicks occur in one direction, again randomly chosen at every injection event. This method preserves momentum only in a time-average sense, but behaves better in the injection regions. The update of gas cell kinetic energy is given by Eqn. 10 of Weinberger+ 2017. In this case, the injection only happens once the available energy exceeds amount, which is quantified in terms of some fraction of the dark matter velocity dispersion times enclosed gas mass, as set by the RadioFeedbackReiorientationFactor parameter. Requires BH_ADIOS_WIND. Do not use together with BH_ADIOS_WIND_DIRECTIONAL. TNG setting: enabled.
BH_ADIOS_DENS_DEP_EFFICIANCY: Make the radiative efficiency density dependent. If this is disabled, the radiative efficiency is constant and set by the BlackHoleRadiativeEfficiency parameter. If this is enabled, the radiative efficiency scales as \(\epsilon_r(\rho) = \epsilon_{r,0} * \rho_{\rm ngb} / (f_{\rm dens} * \rho{\rm thresh}) * (M_{\rm BH}/M_{\rm BH,pivot}^{\gamma})\), up to a maximum value as set by the RadioFeedbackFactorMaxEfficiency parameter. Here \(f_{\rm dens}\) is the parameter RadioFeedbackFactorRefDensityFactor, \(M_{\rm BH,pivot}\) is the parameter RadioFeedbackFactorPivotMass, and \(\gamma\) is the parameter RadioFeedbackFactorSlope. TNG setting: disabled.
BH_ADIOS_WIND_WITH_QUASARTHRESHOLD: Decide the mode switch between quasar-mode feedback and ADIOS wind mode feedback based on the ratio of the SMBH accretion rate relative to the Eddington rate. The value of this ratio is given by the QuasarThreshold parameter. This is Eqn. 1 of Weinberger+ 2017. TNG setting: enabled.
BH_ADIOS_WIND_WITH_VARIABLE_QUASARTHRESHOLD: Scale the threshold (see above) with SMBH mass as \((M_{\rm BH}/M_{\rm ref})^2\), where the value of \(M_{\rm ref} = 10^8 \rm{M}_\odot\) is hard-coded in the code. This is turning on Eqn. 5 of Weinberger+ 2017 for \(\chi(M_{\rm BH})\). Note that the final value of \(\chi\) is capped at 0.1, a constant which is hard-coded in the code. TNG setting: enabled.
BH_ADIOS_ONLY_ABOVE_MINIMUM_DENSITY: “Disable” the wind if the gas density around the SMBH drops below a certain fraction of the star formation threshold density, as set by the RadioFeedbackMinDensityFactor parameter. In fact, the wind is never disabled, instead the efficiency is lowered linearly if the surrounding density is below this chosen density. This only occurs if the SMBH is in the low state. TNG setting: enabled.
BH_CONTINOUS_MODE_SWITCH: Normally, a SMBH is either in the low-state or the high-state mode, and feedback energy is only available in that mode, at a given time. If this option is enabled, then a fraction of the total feedback energy is available to each mode, and both will release energy at a given time. The two fractions always sum to unity. The kinetic (low state) fraction is derived as \(f_k = [(f_{\rm low} \dot{M}_{\rm low})^4 / (f_{\rm high} \dot{M}_{\rm bondi})^4 + (f_{\rm low} \dot{M}_{\rm low})^4] * [1-(\dot{M}_{\rm max} / \dot{M}_{\rm Edd})^4]\) where \(\dot{M}_{\rm low}\) is the Bondi accretion rate divided by the parameter BlackHoleAccretionFactor, \(f_{\rm high}\) is the parameter BlackHoleFeedbackFactor, and \(f_{\rm low}\) is the parameter RadioFeedbackFactor. Here \(\dot{M}_{\rm max}\) is the larger of the two, capped at Eddington. The origin of this equation for \(f_k\) is unclear. TNG setting: disabled.

Spin Evolution and Sub-grid Accretion Disk

A model for sub-grid accretion disks, spin tracking, and spin evolution, as originally presented in Bustamante+ 2019.

BH_SPIN_EVOLUTION: Master switch to enable the SMBH spin evolution model.
BH_SPIN_MODEL: Select the spin model to use. 0 = Prolonged spin model, 1 = Chaotic spin model, 2 - Mass dependent spin model, or 3 = Self-gravity dependent spin model.

Authors

Volker Springel
Mark Vogelsberger
Debora Sijacki
Rainer Weinberger
Sebastian Bustamante
and others

Usage Policy and Citation

Please contact the authors before using this code for a new project. Co-authorship on papers may be requested.

To use aspects of BLACK_HOLES for a new simulation project, including un-modified versions of the “Illustris”, “Auriga”, or “IllustrisTNG” models un-modified, please contact together (at a minimum): Mark, Volker, Lars, and Annalisa.

Papers to cite:

Depends on application.