diff --git a/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst b/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst
index a4ac82e8de3e71c8dbc9d98c247a5ead87cd3cd3..821097f0f4302a1f00ba6ccf6a6f12ef951ad095 100644
--- a/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst
+++ b/doc/RTD/source/SubgridModels/AGNSpinJets/params.rst
@@ -49,36 +49,36 @@ Below we give an example of parameter choices applicable for e.g. a 50 Mpc box.
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
include_jets: 1 # Global switch whether to include jet feedback [1] or not [0].
turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another.
- alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between thin and thick disk, as dot(m) = 0.2 * alpha^2.
- mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
+ alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2.
+ mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1.
AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported.
v_jet_km_p_s: 3160. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s.
- v_jet_cs_ratio: 10. # This sets the jet velocity to v_jet_cs_ratio times the sound speed of the hot gas of the parent halo the black hole is in. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
- v_jet_BH_mass_scaling_reference_mass_Msun: 3.4e3 # The reference mass used in the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
- v_jet_BH_mass_scaling_slope: 0.65 # The slope of the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
- v_jet_mass_loading: 400. # The constant mass loading to use if 'AGN_jet_velocity_model' is 'MassLoading'.
- v_jet_xi: 1. # The numerical multiplier by which the jet velocity formula is scaled, if 'AGN_jet_velocity_model' is 'Local'. If a value of 1 is used, the formulas are used as derived, i.e. they are not rescaled.
+ v_jet_BH_mass_scaling_reference_mass_Msun: 1e9 # The reference mass used in the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
+ v_jet_BH_mass_scaling_slope: 0.5 # The slope of the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
v_jet_min_km_p_s: 500 # The minimal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
- v_jet_max_km_p_s: 1e5 # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
- temperature_hot_gas_min_K: 1e4 # The minimum temperature (in Kelvin) of hot gas to include in the calculation of the hot gas sound speed around the BH. This is only used if 'AGN_jet_velocity_model' is 'Local'.
+ v_jet_max_km_p_s: 1e4 # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector).
- AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
+ AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
- fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests.
+ fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0].
jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector.
- accretion_efficiency: 0.1 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc.
+ accretion_efficiency_mode: Constant # How the accretion efficiencies are calculated for the thick accretion disc. If 'Constant', the value of 'accretion_efficiency_thick' will be used. If 'Vari
+ able', the accretion efficiency will scale with Eddington ratio.
+ accretion_efficiency_thick: 0.01 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid ADIOS winds that take away most of the mass flowing through the accretion disc.
+ accretion_efficiency_slim: 1 # The constant accretion efficiency to use in the slim disc, at super-Eddington rates.
fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0].
radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate.
- TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
+ TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
include_GRMHD_spindown: 1 # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH.
- include_ADIOS_suppression: 0 # Whether to suppress the accretion rate in the fully thick disc regime [1] (Eddington rate below 0.2alpha^2) by the amount expected to be taken away by isotropic kinetic disk winds.
- turn_off_secondary_feedback: 1 # If set to 1, there will be only radiative (thermal) feedback in the thin disk mode, and only jets in the thick disk mode.
- jet_h_r_slope: 1. # The slope of the dependence of jet efficiency on aspect ratio of the subgrid accretion disk, H/R. Default value is 1, and another reasonable value is 0 (same jet efficiency for all disks). Reality could be anything in between. This parameter is only used if turn_off_secondary_feedback is set to 0.
delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0.
- include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
+ include_slim_disk: 1 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
+ use_jets_in_thin_disc: 1 # Whether to use jets alongside radiation in the thin disc at moderate Eddington ratios.
+ use_ADIOS_winds: 1 # Whether to include ADIOS winds in the thick disc as thermal isotropic feedback (same channel as thin disc quasar feedback, but with a different efficiency).
+ slim_disc_wind_factor: 1 # The relative efficiency of slim disc winds at super-Eddington rates. If '1', full winds will be used, while '0' will lead to no winds. Any value in between
+ those can also be used. The wind is implemented in the thermal isotropic feedback channel.
Most of these parameters should work well generally, and should not be changed except for tests. We will discuss only some of the more important ones. You can choose whether to have only the thick and thin disk (low and high BH accretion rates, respectively), or you can also include the slim disk at super-Eddington rates with ``include_slim_disk``. You can control what type of feedback you (do not) want with ``include_jets`` and ``turn_off_radiative_feedback``. If you choose to turn off jets, everything will be modeled as a thin disk (regardless of accretion rate), since jets go hand-in-hand with the thick and the slim disk. Similarly, if you turn off radiation, everything will be treated as a thick disk.
diff --git a/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml b/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml
index 517ae155f55e12f95a63862a7226c31a48ce0928..1f487d9267f5b6eea8fb9979b7672759c89fe0aa 100644
--- a/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml
+++ b/examples/IdealisedCluster/IdealisedCluster_M13/idealised_cluster_M13.yml
@@ -283,24 +283,31 @@ SPINJETAGN:
merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity').
merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length.
- include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1.
- turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1.
- alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range.
- mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
- delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5.
- TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
- N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation.
- jet_velocity_model: Constant # Which jet velocity model to use.
- v_jet_km_p_s: 10000. # Jet velocity if 'jet_velocity_model' is set to 'Constant'.
- opening_angle_in_degrees: 10 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
- eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
- fix_jet_efficiency: 1 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests.
- jet_efficiency: 0.05 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
- AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
- include_slim_disk: 0 # Whether to include additional super-Eddington accretion modeling and jets.
- turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates.
- fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis.
- fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value.
- include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula.
- accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc.
-
+ include_jets: 1 # Global switch whether to include jet feedback [1] or not [0].
+ turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another.
+ alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2.
+ mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
+ seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1.
+ AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported.
+ v_jet_km_p_s: 5000. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s.
+ opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
+ N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector).
+ AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
+ eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
+ fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0].
+ jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
+ fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector.
+ accretion_efficiency_mode: Constant # How the accretion efficiencies are calculated for the thick accretion disc. If 'Constant', the value of 'accretion_efficiency_thick' will be used. If 'Vari
+ able', the accretion efficiency will scale with Eddington ratio.
+ accretion_efficiency_thick: 0.01 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid ADIOS winds that take away most of the mass flowing through the accretion disc.
+ accretion_efficiency_slim: 1 # The constant accretion efficiency to use in the slim disc, at super-Eddington rates.
+ fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0].
+ radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate.
+ TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
+ include_GRMHD_spindown: 1 # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH.
+ delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0.
+ include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
+ use_jets_in_thin_disc: 1 # Whether to use jets alongside radiation in the thin disc at moderate Eddington ratios.
+ use_ADIOS_winds: 0 # Whether to include ADIOS winds in the thick disc as thermal isotropic feedback (same channel as thin disc quasar feedback, but with a different efficiency).
+ slim_disc_wind_factor: 0 # The relative efficiency of slim disc winds at super-Eddington rates. If '1', full winds will be used, while '0' will lead to no winds. Any value in between
+ those can also be used. The wind is implemented in the thermal isotropic feedback channel.
diff --git a/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml b/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml
index 263165007e0078152559f73f47070ce7ed4b7193..0fa12db35be874ca00b3adcedc3621b804824c0b 100644
--- a/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml
+++ b/examples/IdealisedCluster/IdealisedCluster_M135/idealised_cluster_M135.yml
@@ -283,24 +283,31 @@ SPINJETAGN:
merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity').
merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length.
- include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1.
- turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1.
- alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range.
- mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
- delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5.
- TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
- N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation.
- jet_velocity_model: Constant # Which jet velocity model to use.
- v_jet_km_p_s: 10000. # Jet velocity if 'jet_velocity_model' is set to 'Constant'.
- opening_angle_in_degrees: 10 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
- eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
- fix_jet_efficiency: 1 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests.
- jet_efficiency: 0.05 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
- AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
- include_slim_disk: 0 # Whether to include additional super-Eddington accretion modeling and jets.
- turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates.
- fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis.
- fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value.
- include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula.
- accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc.
-
+ include_jets: 1 # Global switch whether to include jet feedback [1] or not [0].
+ turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another.
+ alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2.
+ mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
+ seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1.
+ AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported.
+ v_jet_km_p_s: 5000. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s.
+ opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
+ N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector).
+ AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
+ eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
+ fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0].
+ jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
+ fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector.
+ accretion_efficiency_mode: Constant # How the accretion efficiencies are calculated for the thick accretion disc. If 'Constant', the value of 'accretion_efficiency_thick' will be used. If 'Vari
+ able', the accretion efficiency will scale with Eddington ratio.
+ accretion_efficiency_thick: 0.01 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid ADIOS winds that take away most of the mass flowing through the accretion disc.
+ accretion_efficiency_slim: 1 # The constant accretion efficiency to use in the slim disc, at super-Eddington rates.
+ fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0].
+ radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate.
+ TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
+ include_GRMHD_spindown: 1 # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH.
+ delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0.
+ include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
+ use_jets_in_thin_disc: 1 # Whether to use jets alongside radiation in the thin disc at moderate Eddington ratios.
+ use_ADIOS_winds: 0 # Whether to include ADIOS winds in the thick disc as thermal isotropic feedback (same channel as thin disc quasar feedback, but with a different efficiency).
+ slim_disc_wind_factor: 0 # The relative efficiency of slim disc winds at super-Eddington rates. If '1', full winds will be used, while '0' will lead to no winds. Any value in between
+ those can also be used. The wind is implemented in the thermal isotropic feedback channel.
diff --git a/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml b/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml
index f6c85884538b1eda5d846a6db56936615f448680..3cbd96e9e9dee77114fce299fc93843064c36794 100644
--- a/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml
+++ b/examples/IdealisedCluster/IdealisedCluster_M14/idealised_cluster_M14.yml
@@ -283,23 +283,31 @@ SPINJETAGN:
merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity').
merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length.
- include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1.
- turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1.
- alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range.
- mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
- delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5.
- TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
- N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation.
- jet_velocity_model: Constant # Which jet velocity model to use.
- v_jet_km_p_s: 10000. # Jet velocity if 'jet_velocity_model' is set to 'Constant'.
- opening_angle_in_degrees: 10 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
- eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
- fix_jet_efficiency: 1 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests.
- jet_efficiency: 0.05 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
- AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
- include_slim_disk: 0 # Whether to include additional super-Eddington accretion modeling and jets.
- turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates.
- fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis.
- fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value.
- include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula.
- accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc.
+ include_jets: 1 # Global switch whether to include jet feedback [1] or not [0].
+ turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another.
+ alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2.
+ mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
+ seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1.
+ AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported.
+ v_jet_km_p_s: 5000. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s.
+ opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
+ N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector).
+ AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
+ eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
+ fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0].
+ jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
+ fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector.
+ accretion_efficiency_mode: Constant # How the accretion efficiencies are calculated for the thick accretion disc. If 'Constant', the value of 'accretion_efficiency_thick' will be used. If 'Vari
+ able', the accretion efficiency will scale with Eddington ratio.
+ accretion_efficiency_thick: 0.01 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid ADIOS winds that take away most of the mass flowing through the accretion disc.
+ accretion_efficiency_slim: 1 # The constant accretion efficiency to use in the slim disc, at super-Eddington rates.
+ fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0].
+ radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate.
+ TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
+ include_GRMHD_spindown: 1 # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH.
+ delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0.
+ include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
+ use_jets_in_thin_disc: 1 # Whether to use jets alongside radiation in the thin disc at moderate Eddington ratios.
+ use_ADIOS_winds: 0 # Whether to include ADIOS winds in the thick disc as thermal isotropic feedback (same channel as thin disc quasar feedback, but with a different efficiency).
+ slim_disc_wind_factor: 0 # The relative efficiency of slim disc winds at super-Eddington rates. If '1', full winds will be used, while '0' will lead to no winds. Any value in between
+ those can also be used. The wind is implemented in the thermal isotropic feedback channel.
\ No newline at end of file
diff --git a/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml b/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml
index dba209578c1865b626b67ffaf50007963fb1a4b2..d0367c9cdd330739da26780521d19b5e9e726262 100644
--- a/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml
+++ b/examples/IdealisedCluster/IdealisedCluster_M15/idealised_cluster_M15.yml
@@ -283,23 +283,31 @@ SPINJETAGN:
merger_threshold_type: DynamicalEscapeVelocity # Type of velocity threshold for BH mergers ('CircularVelocity', 'EscapeVelocity', 'DynamicalEscapeVelocity').
merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
minimum_timestep_Myr: 0.1 # Minimum of the accretion-limited time-step length.
- include_jets: 1 # Global switch whether to include jet feedback [1] or not [0]. Requires the 'include_spin' parameter to also be set to 1.
- turn_off_radiative_feedback: 1 # Global switch whether to turn off radiative feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1.
- alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Should be within the 0.1-0.3 range.
- mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
- delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5.
- TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
- N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation.
- jet_velocity_model: Constant # Which jet velocity model to use.
- v_jet_km_p_s: 10000. # Jet velocity if 'jet_velocity_model' is set to 'Constant'.
- opening_angle_in_degrees: 10 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
- eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
- fix_jet_efficiency: 1 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0]. If used, jets will be launched exclusively along the z axis. Should be set to 1 only for tests.
- jet_efficiency: 0.05 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
- AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
- include_slim_disk: 0 # Whether to include additional super-Eddington accretion modeling and jets.
- turn_off_secondary_feedback: 1 # Whether to turn off radiative feedback at very low accretion rates and jets at moderate accretion rates.
- fix_jet_direction: 1 # Whether to fix the jet direction to be along the z-axis.
- fix_radiative_efficiency: 0 # Whether to fix the radiative efficiency to a constant value.
- include_GRMHD_spindown: 0 # Use GRMHD simulation-derived formulas for spindown of BHs, or an analytical formula.
- accretion_efficiency: 1. # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc.
+ include_jets: 1 # Global switch whether to include jet feedback [1] or not [0].
+ turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another.
+ alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2.
+ mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
+ seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1.
+ AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported.
+ v_jet_km_p_s: 5000. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s.
+ opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
+ N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector).
+ AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
+ eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
+ fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0].
+ jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
+ fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector.
+ accretion_efficiency_mode: Constant # How the accretion efficiencies are calculated for the thick accretion disc. If 'Constant', the value of 'accretion_efficiency_thick' will be used. If 'Vari
+ able', the accretion efficiency will scale with Eddington ratio.
+ accretion_efficiency_thick: 0.01 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid ADIOS winds that take away most of the mass flowing through the accretion disc.
+ accretion_efficiency_slim: 1 # The constant accretion efficiency to use in the slim disc, at super-Eddington rates.
+ fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0].
+ radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate.
+ TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
+ include_GRMHD_spindown: 1 # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH.
+ delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0.
+ include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
+ use_jets_in_thin_disc: 1 # Whether to use jets alongside radiation in the thin disc at moderate Eddington ratios.
+ use_ADIOS_winds: 0 # Whether to include ADIOS winds in the thick disc as thermal isotropic feedback (same channel as thin disc quasar feedback, but with a different efficiency).
+ slim_disc_wind_factor: 0 # The relative efficiency of slim disc winds at super-Eddington rates. If '1', full winds will be used, while '0' will lead to no winds. Any value in between
+ those can also be used. The wind is implemented in the thermal isotropic feedback channel.
\ No newline at end of file
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 232eed341614884e687a62ddcc8fecdd42b24b7a..9b5074491c6e6b9401f11dce675570958aefb0ac 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -800,35 +800,35 @@ SPINJETAGN:
include_jets: 1 # Global switch whether to include jet feedback [1] or not [0].
turn_off_radiative_feedback: 0 # Global switch whether to turn off radiative (thermal) feedback [1] or not [0]. This should only be used if 'include_jets' is set to 1, since we want feedback in some form or another.
alpha_acc: 0.2 # Viscous alpha of the subgrid accretion disks. Likely to be within the 0.1-0.3 range. The main effect is that it sets the transition accretion rate between the thin and thick disk, as dot(m) = 0.2 * alpha^2.
- mdot_crit_ADAF: 0.03 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
+ mdot_crit_ADAF: 0.01 # The transition normalized accretion rate (Eddington ratio) at which the disc goes from thick (low accretion rates) to thin (high accretion rates). The feedback also changes from kinetic jets to thermal isotropic, respectively.
seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1.
AGN_jet_velocity_model: Constant # How AGN jet velocities are calculated. If 'Constant', a single value is used. If 'BlackHoleMass', then an empirical relation between halo mass and black hole mass is used to calculate jet velocities. 'HaloMass' is currently not supported.
v_jet_km_p_s: 3160. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s.
- v_jet_cs_ratio: 10. # This sets the jet velocity to v_jet_cs_ratio times the sound speed of the hot gas of the parent halo the black hole is in. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
- v_jet_BH_mass_scaling_reference_mass_Msun: 3.4e3 # The reference mass used in the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
- v_jet_BH_mass_scaling_slope: 0.65 # The slope of the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
- v_jet_mass_loading: 400. # The constant mass loading to use if 'AGN_jet_velocity_model' is 'MassLoading'.
- v_jet_xi: 1. # The numerical multiplier by which the jet velocity formula is scaled, if 'AGN_jet_velocity_model' is 'Local'. If a value of 1 is used, the formulas are used as derived, i.e. they are not rescaled.
+ v_jet_BH_mass_scaling_reference_mass_Msun: 1e9 # The reference mass used in the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
+ v_jet_BH_mass_scaling_slope: 0.5 # The slope of the relation between halo mass and BH mass used to calculate jet velocities. Only used if 'AGN_jet_velocity_model' is 'BlackHoleMass'.
v_jet_min_km_p_s: 500 # The minimal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
- v_jet_max_km_p_s: 1e5 # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
- temperature_hot_gas_min_K: 1e4 # The minimum temperature (in Kelvin) of hot gas to include in the calculation of the hot gas sound speed around the BH. This is only used if 'AGN_jet_velocity_model' or 'AGN_heating_temperature_model' is 'Local'.
+ v_jet_max_km_p_s: 1e4 # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests.
N_jet: 2 # Target number of particles to kick as part of a single jet feedback event. Should be a multiple of 2 to ensure approximate momentum conservation (we always kick particles in pairs, one from each 'side' of the BH, relative to the spin vector).
- AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
+ AGN_jet_feedback_model: MinimumDistance # Which particles to kick from the black hole smoothing kernels. Should be 'SpinAxis', 'MinimumDistance', 'MaximumDistance' or 'MinimumDensity'
eps_f_jet: 1. # Coupling efficiency for jet feedback. No reason to expect this to be less than 1.
fix_jet_efficiency: 0 # Global switch whether to fix jet efficiency to a particular value [1], or use a spin-dependant formula [0].
jet_efficiency: 0.1 # The constant jet efficiency used if 'fix_jet_efficiency' is set to 1.
fix_jet_direction: 0 # Global switch whether to fix the jet direction to be along the z-axis, instead of along the spin vector.
- accretion_efficiency: 0.1 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid winds that take away most of the mass flowing through the accretion disc.
+ accretion_efficiency_mode: Constant # How the accretion efficiencies are calculated for the thick accretion disc. If 'Constant', the value of 'accretion_efficiency_thick' will be used. If 'Vari
+able', the accretion efficiency will scale with Eddington ratio.
+ accretion_efficiency_thick: 0.01 # The accretion efficiency (suppression factor of the accretion rate) to use in the thick disc (ADAF), to represent the effects of subgrid ADIOS winds that take away most of the mass flowing through the accretion disc.
+ accretion_efficiency_slim: 1 # The constant accretion efficiency to use in the slim disc, at super-Eddington rates.
fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0].
radiative_efficiency: 0.1 # The constant jet efficiency used if 'fix_radiative_efficiency' is set to 1. Otherwise, this value is used to define the Eddington accretion rate.
TD_region: B # How to treat the subgrid accretion disk if it is thin, according to the Shakura & Sunyaev (1973) model. If set to B, region b will be used. If set to C, region c will be used.
include_GRMHD_spindown: 1 # Whether to include high jet spindown rates from GRMHD simulations [1], or use an analytical formula that assumes extraction of energy from the rotational mass/energy of the BH.
- turn_off_secondary_feedback: 1 # If set to 1, there will be only radiative (thermal) feedback in the thin disk mode, and only jets in the thick disk mode.
- jet_h_r_slope: 1. # The slope of the dependence of jet efficiency on aspect ratio of the subgrid accretion disk, H/R. Default value is 1, and another reasonable value is 0 (same jet efficiency for all disks). Reality could be anything in between. This parameter is only used if turn_off_secondary_feedback is set to 0.
delta_ADAF: 0.2 # Electron heating parameter, which controls the strength of radiative feedback in thick disks. Should be between 0.1 and 0.5. This parameter is only used if turn_off_secondary_feedback is set to 0.
- include_slim_disk: 0 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
-
+ include_slim_disk: 1 # Global switch whether to include super-Eddington accretion, modeled as the slim disk. If set to 0, disks will be considered thin even at very large accretion rates.
+ use_jets_in_thin_disc: 1 # Whether to use jets alongside radiation in the thin disc at moderate Eddington ratios.
+ use_ADIOS_winds: 1 # Whether to include ADIOS winds in the thick disc as thermal isotropic feedback (same channel as thin disc quasar feedback, but with a different efficiency).
+ slim_disc_wind_factor: 1 # The relative efficiency of slim disc winds at super-Eddington rates. If '1', full winds will be used, while '0' will lead to no winds. Any value in between
+ those can also be used. The wind is implemented in the thermal isotropic feedback channel.
# Parameters related to the neutrinos --------------------------------------------
Neutrino:
diff --git a/src/black_holes/SPIN_JET/black_holes.h b/src/black_holes/SPIN_JET/black_holes.h
index 6b43417bf51017e75084377cc69a19222f4e300f..9ebe99a3dea025bd34bf038a2355e6991a7aae90 100644
--- a/src/black_holes/SPIN_JET/black_holes.h
+++ b/src/black_holes/SPIN_JET/black_holes.h
@@ -58,8 +58,9 @@ __attribute__((always_inline)) INLINE static float black_holes_compute_timestep(
/* Compute instantaneous energy supply rate to the BH energy reservoir
* which is proportional to the BH mass accretion rate */
- const double Energy_rate = bp->radiative_efficiency * props->epsilon_f *
- bp->accretion_rate * c * c;
+ const double Energy_rate =
+ (bp->radiative_efficiency * props->epsilon_f + bp->wind_efficiency) *
+ bp->accretion_rate * c * c;
/* Compute instantaneous jet energy supply rate to the BH jet reservoir
* which is proportional to the BH mass accretion rate */
@@ -150,7 +151,6 @@ __attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
bp->AGN_number_of_AGN_events = 0;
bp->AGN_number_of_energy_injections = 0;
bp->AGN_cumulative_energy = 0.f;
- bp->aspect_ratio = 0.01f;
bp->jet_efficiency = 0.1f;
bp->radiative_efficiency = 0.1f;
bp->accretion_disk_angle = 0.01f;
@@ -158,17 +158,24 @@ __attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
bp->eddington_fraction = 0.01f;
bp->jet_reservoir = 0.f;
bp->total_jet_energy = 0.f;
+ bp->wind_efficiency = 0.f;
+ bp->wind_energy = 0.f;
+ bp->radiated_energy = 0.f;
+ bp->accretion_efficiency = 1.f;
bp->dt_jet = 0.f;
bp->dt_ang_mom = 0.f;
bp->AGN_number_of_AGN_jet_events = 0;
bp->AGN_number_of_jet_injections = 0;
- bp->group_mass = 0.f;
for (int i = 0; i < BH_accretion_modes_count; ++i)
bp->accreted_mass_by_mode[i] = 0.f;
for (int i = 0; i < BH_accretion_modes_count; ++i)
bp->thermal_energy_by_mode[i] = 0.f;
for (int i = 0; i < BH_accretion_modes_count; ++i)
bp->jet_energy_by_mode[i] = 0.f;
+ for (int i = 0; i < BH_accretion_modes_count; ++i)
+ bp->wind_energy_by_mode[i] = 0.f;
+ for (int i = 0; i < BH_accretion_modes_count; ++i)
+ bp->radiated_energy_by_mode[i] = 0.f;
bp->jet_direction[0] = bp->angular_momentum_direction[0];
bp->jet_direction[1] = bp->angular_momentum_direction[1];
bp->jet_direction[2] = bp->angular_momentum_direction[2];
@@ -558,7 +565,7 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_bpart(
}
/* Evolve the black hole spin. */
- merger_spin_evolve(bpi, bpj, constants);
+ black_hole_merger_spin_evolve(bpi, bpj, constants);
/* Increase the masses of the BH. */
bpi->mass += bpj->mass;
@@ -585,19 +592,19 @@ __attribute__((always_inline)) INLINE static void black_holes_swallow_bpart(
dot_product = 0.;
}
- if (j_BH(bpi, constants) * dot_product <
- -0.5 * j_warp(bpi, constants, props)) {
+ if (black_hole_angular_momentum_magnitude(bpi, constants) * dot_product <
+ -0.5 * black_hole_warp_angular_momentum(bpi, constants, props)) {
bpi->spin = -1. * fabsf(bpi->spin);
} else {
bpi->spin = fabsf(bpi->spin);
}
/* Update various quantities with new spin */
- decide_mode(bpi, props);
- bpi->aspect_ratio = aspect_ratio(bpi, constants, props);
+ black_hole_select_accretion_mode(bpi, props);
bpi->accretion_disk_angle = dot_product;
- bpi->radiative_efficiency = rad_efficiency(bpi, props);
- bpi->jet_efficiency = jet_efficiency(bpi, props);
+ bpi->radiative_efficiency = black_hole_radiative_efficiency(bpi, props);
+ bpi->jet_efficiency = black_hole_jet_efficiency(bpi, props);
+ bpi->wind_efficiency = black_hole_wind_efficiency(bpi, props);
/* Collect the swallowed angular momentum */
bpi->swallowed_angular_momentum[0] += bpj->swallowed_angular_momentum[0];
@@ -818,10 +825,12 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
bp->f_visc = 1.0;
}
- /* Compute the Eddington rate (internal units) */
+ /* Compute the Eddington rate (internal units). The radiative efficiency
+ is taken to be equal to the Novikov-Thorne (1973) one, and thus varies
+ with spin between 4% and 40%. */
const double Eddington_rate =
4. * M_PI * G * BH_mass * proton_mass /
- (props->radiative_efficiency * c * sigma_Thomson);
+ (eps_Novikov_Thorne(bp->spin) * c * sigma_Thomson);
/* Should we record this time as the most recent high accretion rate? */
if (accr_rate > props->f_Edd_recording * Eddington_rate) {
@@ -833,17 +842,24 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
}
/* Limit the accretion rate to a fraction of the Eddington rate */
- accr_rate = min(accr_rate, props->f_Edd * Eddington_rate);
- bp->accretion_rate = accr_rate;
bp->eddington_fraction = accr_rate / Eddington_rate;
+ accr_rate = min(accr_rate, props->f_Edd * Eddington_rate);
- /* If we are in the thick disc mode, suppress the accretion rate by the
- * accretion efficiency. */
- if (bp->accretion_mode == BH_thick_disc) {
- accr_rate *= props->accretion_efficiency;
- bp->accretion_rate *= props->accretion_efficiency;
- bp->eddington_fraction *= props->accretion_efficiency;
- }
+ /* Include the effects of the accretion efficiency */
+ bp->accretion_rate = accr_rate * bp->accretion_efficiency;
+ bp->eddington_fraction *= bp->accretion_efficiency;
+
+ /* Change mode based on new accretion rates if necessary, and update
+ all important quantities. */
+ black_hole_select_accretion_mode(bp, props);
+ bp->accretion_efficiency =
+ black_hole_accretion_efficiency(bp, props, constants, cosmo);
+ bp->accretion_rate = accr_rate * bp->accretion_efficiency;
+ bp->eddington_fraction = bp->accretion_rate / Eddington_rate;
+
+ bp->radiative_efficiency = black_hole_radiative_efficiency(bp, props);
+ bp->jet_efficiency = black_hole_jet_efficiency(bp, props);
+ bp->wind_efficiency = black_hole_wind_efficiency(bp, props);
/* Define feedback-related quantities that we will update and need later on */
double luminosity = 0.;
@@ -874,8 +890,8 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
/* Decide if accretion is prograde (spin positive) or retrograde
(spin negative) based on condition from King et al. (2005) */
- if ((j_BH(bp, constants) * dot_product <
- -0.5 * j_warp(bp, constants, props)) &&
+ if ((black_hole_angular_momentum_magnitude(bp, constants) * dot_product <
+ -0.5 * black_hole_warp_angular_momentum(bp, constants, props)) &&
(fabsf(bp->spin) > 0.01)) {
bp->spin = -1. * fabsf(bp->spin);
} else {
@@ -886,7 +902,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
step */
double n_i = 0.;
if (bp->accretion_rate > 0.) {
- n_i = delta_m_0 / m_warp(bp, constants, props);
+ n_i = delta_m_0 / black_hole_warp_mass(bp, constants, props);
}
/* Update the angular momentum vector of the BH based on how many
@@ -903,16 +919,19 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
bp->angular_momentum_direction[2] =
bp->spec_angular_momentum_gas[2] / spec_ang_mom_norm;
} else {
+ const double j_warp =
+ black_hole_warp_angular_momentum(bp, constants, props);
+ const double j_BH = black_hole_angular_momentum_magnitude(bp, constants);
const double ang_mom_total[3] = {
- bp->angular_momentum_direction[0] * j_BH(bp, constants) +
+ bp->angular_momentum_direction[0] * j_BH +
n_i * bp->spec_angular_momentum_gas[0] / spec_ang_mom_norm *
- j_warp(bp, constants, props),
- bp->angular_momentum_direction[1] * j_BH(bp, constants) +
+ j_warp,
+ bp->angular_momentum_direction[1] * j_BH +
n_i * bp->spec_angular_momentum_gas[1] / spec_ang_mom_norm *
- j_warp(bp, constants, props),
- bp->angular_momentum_direction[2] * j_BH(bp, constants) +
+ j_warp,
+ bp->angular_momentum_direction[2] * j_BH +
n_i * bp->spec_angular_momentum_gas[2] / spec_ang_mom_norm *
- j_warp(bp, constants, props)};
+ j_warp};
/* Modulus of the new J_BH */
const double modulus = sqrt(ang_mom_total[0] * ang_mom_total[0] +
@@ -934,11 +953,15 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
bp->angular_momentum_direction[2] = 1;
}
+ bp->jet_direction[0] = bp->angular_momentum_direction[0];
+ bp->jet_direction[1] = bp->angular_momentum_direction[1];
+ bp->jet_direction[2] = bp->angular_momentum_direction[2];
+
float spin_final = -1.;
/* Calculate the change in the BH spin */
if (bp->subgrid_mass > 0.) {
spin_final = bp->spin + delta_m_0 / bp->subgrid_mass *
- da_dln_mbh_0(bp, constants, props);
+ black_hole_spinup_rate(bp, constants, props);
} else {
error(
"Black hole with id %lld tried to evolve spin with zero "
@@ -969,28 +992,32 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
/* Update other quantities */
bp->accretion_disk_angle = dot_product;
- bp->aspect_ratio = aspect_ratio(bp, constants, props);
/* Update efficiencies given the new spin */
- bp->radiative_efficiency = rad_efficiency(bp, props);
- bp->jet_efficiency = jet_efficiency(bp, props);
+ bp->radiative_efficiency = black_hole_radiative_efficiency(bp, props);
+ bp->jet_efficiency = black_hole_jet_efficiency(bp, props);
+ bp->wind_efficiency = black_hole_wind_efficiency(bp, props);
/* Final jet power at the end of the step */
- jet_power = bp->jet_efficiency * accr_rate * c * c;
+ jet_power = bp->jet_efficiency * bp->accretion_rate * c * c;
/* Final luminosity at the end of the step */
- luminosity = bp->radiative_efficiency * accr_rate * c * c;
+ luminosity = bp->radiative_efficiency * bp->accretion_rate * c * c;
/* The amount of mass the BH is actually swallowing, including the
effects of the updated efficiencies. */
const double delta_m_real =
- delta_m_0 * (1. - rad_efficiency(bp, props) - jet_efficiency(bp, props));
+ delta_m_0 * (1. - black_hole_radiative_efficiency(bp, props) -
+ black_hole_jet_efficiency(bp, props) -
+ black_hole_wind_efficiency(bp, props));
/* Increase the reservoir */
- bp->jet_reservoir +=
- delta_m_0 * c * c * props->eps_f_jet * jet_efficiency(bp, props);
+ bp->jet_reservoir += delta_m_0 * c * c * props->eps_f_jet *
+ black_hole_jet_efficiency(bp, props);
bp->energy_reservoir +=
- delta_m_0 * c * c * props->epsilon_f * rad_efficiency(bp, props);
+ delta_m_0 * c * c *
+ (props->epsilon_f * black_hole_radiative_efficiency(bp, props) +
+ black_hole_wind_efficiency(bp, props));
/* Update the masses */
bp->subgrid_mass += delta_m_real;
@@ -999,6 +1026,14 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
/* Update the total accreted masses split by accretion mode of the BHs */
bp->accreted_mass_by_mode[bp->accretion_mode] += delta_m_real;
+ /* Update total energies launched as radiation/winds, and by mode */
+ bp->radiated_energy += delta_m_0 * c * c * bp->radiative_efficiency;
+ bp->radiated_energy_by_mode[bp->accretion_mode] +=
+ delta_m_0 * c * c * bp->radiative_efficiency;
+ bp->wind_energy_by_mode[bp->accretion_mode] +=
+ delta_m_0 * c * c * bp->wind_efficiency;
+ bp->wind_energy += delta_m_0 * c * c * bp->wind_efficiency;
+
if (bp->subgrid_mass < 0.) {
warning(
"Black hole %lld has reached a negative mass (%f) due"
@@ -1023,7 +1058,8 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
* in proprtion to its current accretion rate. We do not account for this
* in the swallowing approach, however. */
bp->mass -=
- (bp->radiative_efficiency + bp->jet_efficiency) * accr_rate * dt;
+ (bp->radiative_efficiency + bp->jet_efficiency + bp->wind_efficiency) *
+ bp->accretion_rate * dt;
if (bp->mass < 0)
error(
"Black hole %lld has reached a negative mass (%f). This is "
@@ -1311,13 +1347,18 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
/* Decide the accretion mode of the BH, based on the new spin and Eddington
* fraction */
- decide_mode(bp, props);
+ black_hole_select_accretion_mode(bp, props);
+ bp->accretion_efficiency =
+ black_hole_accretion_efficiency(bp, props, constants, cosmo);
+ bp->accretion_rate = accr_rate * bp->accretion_efficiency;
+ bp->eddington_fraction = bp->accretion_rate / Eddington_rate;
/* The accretion/feedback mode is now possibly different, so the feedback
efficiencies need to be updated. This is important for computing the
next time-step of the BH. */
- bp->radiative_efficiency = rad_efficiency(bp, props);
- bp->jet_efficiency = jet_efficiency(bp, props);
+ bp->radiative_efficiency = black_hole_radiative_efficiency(bp, props);
+ bp->jet_efficiency = black_hole_jet_efficiency(bp, props);
+ bp->wind_efficiency = black_hole_wind_efficiency(bp, props);
/* Calculate a BH angular momentum evolution time step. Two conditions are
used, one ensures that the BH spin changes by a small amount over the
@@ -1336,7 +1377,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
epsilon_spin = 0.01 * spin_magnitude;
}
float dt_ang_mom = epsilon_spin /
- fabsf(da_dln_mbh_0(bp, constants, props)) *
+ fabsf(black_hole_spinup_rate(bp, constants, props)) *
bp->subgrid_mass / bp->accretion_rate;
/* We now compute the angular-momentum (direction) time-step. We allow
@@ -1348,20 +1389,25 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
/* Angular momentum direction of gas in kernel along the direction of
the current spin vector */
- const float cosine = (bp->spec_angular_momentum_gas[0] *
- bp->angular_momentum_direction[0] +
- bp->spec_angular_momentum_gas[1] *
- bp->angular_momentum_direction[1] +
- bp->spec_angular_momentum_gas[2] *
- bp->angular_momentum_direction[2]) /
- spec_ang_mom_norm;
- /* Compute sine, i.e. the componenent perpendicular to that. */
- const float sine = fmaxf(0., sqrtf(1. - cosine * cosine));
-
- const float dt_redirection = 0.1 * m_warp(bp, constants, props) /
- bp->accretion_rate * j_BH(bp, constants) /
- j_warp(bp, constants, props) / sine;
- dt_ang_mom = min(dt_ang_mom, dt_redirection);
+
+ if (spec_ang_mom_norm > 0.) {
+ const float cosine = (bp->spec_angular_momentum_gas[0] *
+ bp->angular_momentum_direction[0] +
+ bp->spec_angular_momentum_gas[1] *
+ bp->angular_momentum_direction[1] +
+ bp->spec_angular_momentum_gas[2] *
+ bp->angular_momentum_direction[2]) /
+ spec_ang_mom_norm;
+ /* Compute sine, i.e. the componenent perpendicular to that. */
+ const float sine = fmaxf(0., sqrtf(1. - cosine * cosine));
+
+ const float dt_redirection =
+ 0.1 * black_hole_warp_mass(bp, constants, props) *
+ black_hole_angular_momentum_magnitude(bp, constants) /
+ (bp->accretion_rate *
+ black_hole_warp_angular_momentum(bp, constants, props) * sine);
+ dt_ang_mom = min(dt_ang_mom, dt_redirection);
+ }
}
bp->dt_ang_mom = dt_ang_mom;
@@ -1576,7 +1622,6 @@ INLINE static void black_holes_create_from_gas(
bp->jet_direction[1] = bp->angular_momentum_direction[1];
bp->jet_direction[2] = bp->angular_momentum_direction[2];
- bp->aspect_ratio = 0.01f;
bp->jet_efficiency = 0.1f;
bp->radiative_efficiency = 0.1f;
bp->accretion_disk_angle = 0.01f;
@@ -1584,6 +1629,10 @@ INLINE static void black_holes_create_from_gas(
bp->eddington_fraction = 0.01f;
bp->jet_reservoir = 0.f;
bp->total_jet_energy = 0.f;
+ bp->wind_efficiency = 0.f;
+ bp->wind_energy = 0.f;
+ bp->radiated_energy = 0.f;
+ bp->accretion_efficiency = 1.f;
bp->dt_jet = 0.f;
bp->dt_ang_mom = 0.f;
bp->delta_T = black_hole_feedback_delta_T(bp, props, cosmo, constants);
@@ -1596,6 +1645,10 @@ INLINE static void black_holes_create_from_gas(
bp->thermal_energy_by_mode[i] = 0.f;
for (int i = 0; i < BH_accretion_modes_count; ++i)
bp->jet_energy_by_mode[i] = 0.f;
+ for (int i = 0; i < BH_accretion_modes_count; ++i)
+ bp->wind_energy_by_mode[i] = 0.f;
+ for (int i = 0; i < BH_accretion_modes_count; ++i)
+ bp->radiated_energy_by_mode[i] = 0.f;
/* We haven't accreted anything yet */
bp->total_accreted_mass = 0.f;
@@ -1636,16 +1689,4 @@ INLINE static void black_holes_create_from_gas(
black_holes_mark_bpart_as_not_swallowed(&bp->merger_data);
}
-/**
- * @brief Store the halo mass in the fof algorithm for the black
- * hole particle.
- *
- * @param p_data The black hole particle data.
- * @param halo_mass The halo mass to update.
- */
-__attribute__((always_inline)) INLINE static void black_holes_update_halo_mass(
- struct bpart* bp, float halo_mass) {
- bp->group_mass = halo_mass;
-}
-
#endif /* SWIFT_SPIN_JET_BLACK_HOLES_H */
diff --git a/src/black_holes/SPIN_JET/black_holes_iact.h b/src/black_holes/SPIN_JET/black_holes_iact.h
index 8fb0e0d6815e22af0b5e29ebd2a113799fcfb31f..3c1db7894b9acca76d2cc3243cdc589f9d58e91d 100644
--- a/src/black_holes/SPIN_JET/black_holes_iact.h
+++ b/src/black_holes/SPIN_JET/black_holes_iact.h
@@ -1062,10 +1062,11 @@ runner_iact_nonsym_bh_gas_feedback(
bi->id, pj->id, xpj->v_full[0], xpj->v_full[1], xpj->v_full[2]);
#endif
- /* Store the jet energy */
+ /* Store the jet energy and other variables of interest */
const double delta_energy_jet = delta_u_jet * hydro_get_mass(pj);
tracers_after_jet_feedback(pj, xpj, with_cosmology, cosmo->a, time,
- delta_energy_jet, vel_kick);
+ delta_energy_jet, vel_kick, bi->accretion_mode,
+ bi->id);
/* Impose maximal viscosity */
hydro_diffusive_feedback_reset(pj);
diff --git a/src/black_holes/SPIN_JET/black_holes_io.h b/src/black_holes/SPIN_JET/black_holes_io.h
index 56851bf06b33f834f14c550401e15c398b033eee..825ddddf08c3075e6ee469a66a27859102476115 100644
--- a/src/black_holes/SPIN_JET/black_holes_io.h
+++ b/src/black_holes/SPIN_JET/black_holes_io.h
@@ -176,7 +176,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
const int with_cosmology) {
/* Say how much we want to write */
- *num_fields = 56;
+ *num_fields = 59;
/* List what we want to write */
list[0] = io_make_output_field_convert_bpart(
@@ -416,7 +416,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
"AGNTotalInjectedEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, bparts,
AGN_cumulative_energy,
"Total (cumulative) physical energies injected into gas particles "
- "in AGN feedback.");
+ "in AGN feedback, including the effects of both radiation and winds.");
list[36] = io_make_output_field_convert_bpart(
"Potentials", FLOAT, 1, UNIT_CONV_POTENTIAL, -1.f, bparts,
@@ -433,102 +433,122 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
"Direction of the black hole spin vector, normalised to unity.");
list[39] = io_make_output_field(
- "AccretionDiscAspectRatios", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
- aspect_ratio,
- "The aspect ratio, h/r, of the subgrid accretion disc "
- "around the black hole.");
-
- list[40] = io_make_output_field(
"JetEfficiencies", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
jet_efficiency, "Jet power divided by accretion rate.");
- list[41] = io_make_output_field(
+ list[40] = io_make_output_field(
"RadiativeEfficiencies", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
radiative_efficiency, "AGN luminosity divided by accretion rate.");
- list[42] = io_make_output_field("CosAccretionDiskAngle", FLOAT, 1,
+ list[41] = io_make_output_field("CosAccretionDiskAngle", FLOAT, 1,
UNIT_CONV_NO_UNITS, 0.f, bparts,
accretion_disk_angle,
"Cosine of the angle between the spin vector "
"and the accreting gas angular momentum.");
- list[43] = io_make_output_field(
+ list[42] = io_make_output_field(
"AccretionModes", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts, accretion_mode,
"Accretion flow regime. 0 - Thick disk, 1 - Thin disk, 2 - Slim disk");
- list[44] = io_make_output_field(
+ list[43] = io_make_output_field(
"JetReservoir", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, bparts, jet_reservoir,
"Total jet energy waiting to be released (once it "
"grows large enough to kick a single particle).");
- list[45] = io_make_output_field(
+ list[44] = io_make_output_field(
"InjectedJetEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, bparts,
total_jet_energy, "Total jet energy injected into AGN surroundings.");
- list[46] = io_make_output_field(
+ list[45] = io_make_output_field(
"JetTimeSteps", FLOAT, 1, UNIT_CONV_TIME, 0.f, bparts, dt_jet,
"Jet-launching-limited time-steps of black holes.");
- list[47] = io_make_output_field(
+ list[46] = io_make_output_field(
"NumberOfJetParticlesLaunched", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
AGN_number_of_jet_injections,
"Integer number of (kinetic) energy injections the black hole has had "
"so far");
- list[48] = io_make_output_field(
+ list[47] = io_make_output_field(
"NumberOfAGNJetEvents", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
AGN_number_of_AGN_jet_events,
"Integer number of AGN jet launching events the black hole has had"
" (the number of times the BH did AGN jet feedback)");
if (with_cosmology) {
- list[49] = io_make_output_field(
+ list[48] = io_make_output_field(
"LastAGNJetScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
last_AGN_jet_event_scale_factor,
"Scale-factors at which the black holes last had an AGN jet event.");
} else {
- list[49] = io_make_output_field(
+ list[48] = io_make_output_field(
"LastAGNJetTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f, bparts,
last_AGN_jet_event_time,
"Times at which the black holes last had an AGN jet event.");
}
- list[50] = io_make_output_field(
+ list[49] = io_make_output_field(
"EddingtonFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
eddington_fraction,
"Accretion rates of black holes in units of their Eddington rates. "
"This is based on the unlimited accretion rates, so these fractions "
"can be above the limiting fEdd.");
- list[51] = io_make_output_field(
- "FOFGroupMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f, bparts, group_mass,
- "Parent halo masses of the black holes, as determined from the FOF "
- "algorithm.");
-
- list[52] =
+ list[50] =
io_make_output_field("JetVelocities", FLOAT, 1, UNIT_CONV_VELOCITY, 0.f,
bparts, v_jet, "The current jet velocities.");
- list[53] = io_make_output_field(
- "TotalAccretedMassesByMode", FLOAT, 3, UNIT_CONV_MASS, 0.f, bparts,
- accreted_mass_by_mode,
+ list[51] = io_make_output_field(
+ "TotalAccretedMassesByMode", FLOAT, BH_accretion_modes_count,
+ UNIT_CONV_MASS, 0.f, bparts, accreted_mass_by_mode,
"The total accreted mass in each accretion mode. The components to the "
"mass accreted in the thick, thin and slim disc modes, respectively.");
- list[54] = io_make_output_field(
- "AGNTotalInjectedEnergiesByMode", FLOAT, 3, UNIT_CONV_ENERGY, 0.f, bparts,
- thermal_energy_by_mode,
- "The total energy injected in the thermal AGN feedback mode, split by "
+ list[52] = io_make_output_field(
+ "AGNTotalInjectedEnergiesByMode", FLOAT, BH_accretion_modes_count,
+ UNIT_CONV_ENERGY, 0.f, bparts, thermal_energy_by_mode,
+ "The total energy injected in the thermal AGN feedback mode, including "
+ "the contributions of both radiation and wind feedback, split by "
"accretion mode. The components correspond to the thermal energy dumped "
"in the thick, thin and slim disc modes, respectively.");
- list[55] = io_make_output_field(
- "InjectedJetEnergiesByMode", FLOAT, 3, UNIT_CONV_ENERGY, 0.f, bparts,
- jet_energy_by_mode,
+ list[53] = io_make_output_field(
+ "InjectedJetEnergiesByMode", FLOAT, BH_accretion_modes_count,
+ UNIT_CONV_ENERGY, 0.f, bparts, jet_energy_by_mode,
"The total energy injected in the kinetic jet AGN feedback mode, split "
- "by accretion mode. The components correspond to the thermal energy "
+ "by accretion mode. The components correspond to the jet energy "
"dumped in the thick, thin and slim disc modes, respectively.");
+ list[54] = io_make_output_field(
+ "WindEfficiencies", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+ wind_efficiency, "The wind efficiencies of the black holes.");
+
+ list[55] = io_make_output_field(
+ "TotalRadiatedEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, bparts,
+ radiated_energy,
+ "The total energy launched into radiation by the black holes, "
+ "in all accretion modes. ");
+
+ list[56] = io_make_output_field(
+ "RadiatedEnergiesByMode", FLOAT, BH_accretion_modes_count,
+ UNIT_CONV_ENERGY, 0.f, bparts, radiated_energy_by_mode,
+ "The total energy launched into radiation by the black holes, split "
+ "by accretion mode. The components correspond to the radiative energy "
+ "dumped in the thick, thin and slim disc modes, respectively.");
+
+ list[57] = io_make_output_field(
+ "TotalWindEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, bparts, wind_energy,
+ "The total energy launched into accretion disc winds by the black "
+ "holes, in all accretion modes. ");
+
+ list[58] = io_make_output_field(
+ "WindEnergiesByMode", FLOAT, BH_accretion_modes_count, UNIT_CONV_ENERGY,
+ 0.f, bparts, wind_energy_by_mode,
+ "The total energy launched into accretion disc winds by the black "
+ "holes, split by accretion mode. The components correspond to the "
+ "radiative energy dumped in the thick, thin and slim disc modes, "
+ "respectively.");
+
#ifdef DEBUG_INTERACTIONS_BLACK_HOLES
list += *num_fields;
diff --git a/src/black_holes/SPIN_JET/black_holes_part.h b/src/black_holes/SPIN_JET/black_holes_part.h
index b26ee1740389b9ed8964ec99d1913af0c15f45b3..7fda793927eb4be44f67eb29a39e5eb28feea100 100644
--- a/src/black_holes/SPIN_JET/black_holes_part.h
+++ b/src/black_holes/SPIN_JET/black_holes_part.h
@@ -30,12 +30,12 @@
#include "timeline.h"
/*! The possible accretion modes every black hole can take. */
-enum BH_accretion_modes {
+typedef enum BH_accretion_modes {
BH_thick_disc = 0, /* At low Eddington ratios */
BH_thin_disc, /* At moderate Eddington ratios */
BH_slim_disc, /* Super-Eddington accretion */
BH_accretion_modes_count /* Number of possible accretion modes */
-};
+} BH_accretion_modes;
/**
* @brief Particle fields for the black hole particles.
@@ -230,9 +230,6 @@ struct bpart {
of the gas in the smoothing kernel */
float accretion_disk_angle;
- /*! Aspect ratio of the subgrid accretion disk */
- float aspect_ratio;
-
/*! Which type is the subgrid accretion disk (thick, thin or slim) */
enum BH_accretion_modes accretion_mode;
@@ -245,17 +242,31 @@ struct bpart {
/*! The current jet kick velocity to be applied */
float v_jet;
- /*! The energ in the jet reservoir */
+ /*! The energy in the jet reservoir */
float jet_reservoir;
/*! Total jet energy launched so far */
float total_jet_energy;
+ /*! Efficiency of wind launching */
+ float wind_efficiency;
+
+ /*! Energy launched into radiation */
+ float radiated_energy;
+
+ /*! Energy launched as winds */
+ float wind_energy;
+
+ /*! Accretion efficiency of this BH */
+ float accretion_efficiency;
+
/*! Total accreted masses, radiated energies and jet energies launched
by BHs, split by accretion mode */
float accreted_mass_by_mode[BH_accretion_modes_count];
float thermal_energy_by_mode[BH_accretion_modes_count];
float jet_energy_by_mode[BH_accretion_modes_count];
+ float wind_energy_by_mode[BH_accretion_modes_count];
+ float radiated_energy_by_mode[BH_accretion_modes_count];
/*! Total number of jet kicks */
int AGN_number_of_jet_injections;
@@ -263,9 +274,6 @@ struct bpart {
/*! Total number of jet kicking events */
int AGN_number_of_AGN_jet_events;
- /*! Halo mass the black hole is assigned to */
- float group_mass;
-
union {
/*! Last time a jet event occurred */
diff --git a/src/black_holes/SPIN_JET/black_holes_properties.h b/src/black_holes/SPIN_JET/black_holes_properties.h
index d61cc7ecdf8a15a894ebbe86af8247e4de874b0e..fe661a90a84846cbc9f291c819b7e79752fe4507 100644
--- a/src/black_holes/SPIN_JET/black_holes_properties.h
+++ b/src/black_holes/SPIN_JET/black_holes_properties.h
@@ -65,6 +65,11 @@ enum thin_disc_regions {
TD_region_C /*< Region C from Shakura & Sunyaev (1973) */
};
+enum accretion_efficiency_modes {
+ BH_accretion_efficiency_constant, /*< Single number */
+ BH_accretion_efficiency_variable /*< Scaling with Eddington ratio */
+};
+
/**
* @brief Properties of black holes and AGN feedback in the EAGEL model.
*/
@@ -332,16 +337,13 @@ struct black_holes_props {
*/
int turn_off_radiative_feedback;
- /*! Global switch for whether to turn off radiation in the thick disk and
- jets in the thin disk [1] or not [0] */
- int turn_off_secondary_feedback;
-
/* Whether we want to include super-Eddington accretion, modeled as the slim
disk */
int include_slim_disk;
- /* Whether to use GRMHD fits for the spindown rate due to jets */
- int include_GRMHD_spindown;
+ /* Whether or not to use jets from the thin disc regime (at moderate
+ * Eddington ratios. */
+ int use_jets_in_thin_disc;
/*! Whether to fix the radiative efficiency to some value [1] or not [0]. */
int fix_radiative_efficiency;
@@ -360,10 +362,6 @@ struct black_holes_props {
/*! Jet velocity if the constant velocity model is used */
float v_jet;
- /*! Parameters of the scaling between AGN jet velocity and BH mass */
- float v_jet_BH_mass_scaling_reference_mass;
- float v_jet_BH_mass_scaling_slope;
-
/*! Sets the launching velocity of the jet to v_jet_cs_ratio times the
sound speed of the hot gas in the halo, assuming it is at virial
temperature. This is used if the launching model is BH_mass or
@@ -378,6 +376,14 @@ struct black_holes_props {
sound_speed launching models are used. */
float v_jet_xi;
+ /*! The reference BH mass to use in the case that we employ a BH mass scaling
+ * for the jet velocity. */
+ float v_jet_BH_mass_scaling_reference_mass;
+
+ /*! The power law slope to use in the case that we employ a BH mass scaling
+ * for the jet velocity. */
+ float v_jet_BH_mass_scaling_slope;
+
/*! The minimal jet velocity to use in the variable-velocity models */
float v_jet_min;
@@ -391,11 +397,6 @@ struct black_holes_props {
/*! The effective (half-)opening angle of the jet. */
float opening_angle;
- /*! The slope of the dependence of jet efficiency on aspect ratio of the
- subgrid accretion disk, H/R. Default value is 1, and another reasonable
- value is 0 (same jet efficiency for all disks). */
- float jet_h_r_slope;
-
/*! The coupling efficiency for jet feedback. */
float eps_f_jet;
@@ -410,11 +411,31 @@ struct black_holes_props {
/*! The accretion efficiency (suppression of accretion rate) to use in
* the thick disc regime (at low Eddington ratios). */
- float accretion_efficiency;
+ float accretion_efficiency_thick;
+
+ /*! The accretion efficiency (suppression of accretion rate) to use in
+ * the slim disc regime (at super-Eddington ratios). */
+ float accretion_efficiency_slim;
+
+ /*! Expontent to use for scaling of accretion efficiency with transition
+ * radius in the thick disc. */
+ float ADIOS_s;
+
+ /* Whether or not we want to use wind feedback in the ADAF/ADIOS regime
+ (at low Eddington ratios). */
+ int use_ADIOS_winds;
+
+ /* The factor by which we multiply the slim disc wind efficiency - 0 meaning
+ no winds, and 1 meaning full winds. */
+ float slim_disc_wind_factor;
/*! The jet launching scheme to use: minimum distance,
maximum distance, closest to spin axis or minimum density. */
enum AGN_jet_feedback_models jet_feedback_model;
+
+ /*! The accretion efficiency mode to use: constant or variable
+ * (Eddington-ratio dependent) . */
+ enum accretion_efficiency_modes accretion_efficiency_mode;
};
/**
@@ -843,17 +864,6 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
"least one of the two feedback modes must be turned on.");
}
- bp->turn_off_secondary_feedback =
- parser_get_param_int(params, "SPINJETAGN:turn_off_secondary_feedback");
-
- if ((bp->turn_off_secondary_feedback != 0) &&
- (bp->turn_off_secondary_feedback != 1)) {
- error(
- "The turn_off_secondary_feedback parameter must be either 0 or 1, "
- "not %d",
- bp->turn_off_secondary_feedback);
- }
-
bp->include_slim_disk =
parser_get_param_int(params, "SPINJETAGN:include_slim_disk");
@@ -864,14 +874,14 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->include_slim_disk);
}
- bp->include_GRMHD_spindown =
- parser_get_param_int(params, "SPINJETAGN:include_GRMHD_spindown");
+ bp->use_jets_in_thin_disc =
+ parser_get_param_int(params, "SPINJETAGN:use_jets_in_thin_disc");
- if ((bp->include_GRMHD_spindown != 0) && (bp->include_GRMHD_spindown != 1)) {
+ if ((bp->use_jets_in_thin_disc != 0) && (bp->use_jets_in_thin_disc != 1)) {
error(
- "The include_GRMHD_spindown parameter must be either 0 or 1, "
+ "The use_jets_in_thin_disc parameter must be either 0 or 1, "
"not %d",
- bp->include_GRMHD_spindown);
+ bp->use_jets_in_thin_disc);
}
bp->N_jet = parser_get_param_float(params, "SPINJETAGN:N_jet");
@@ -902,13 +912,14 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->v_jet_BH_mass_scaling_slope = parser_get_param_float(
params, "SPINJETAGN:v_jet_BH_mass_scaling_slope");
- bp->v_jet_cs_ratio =
- parser_get_param_float(params, "SPINJETAGN:v_jet_cs_ratio");
-
bp->v_jet_min =
parser_get_param_float(params, "SPINJETAGN:v_jet_min_km_p_s");
bp->v_jet_min *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs));
+ bp->v_jet_max =
+ parser_get_param_float(params, "SPINJETAGN:v_jet_max_km_p_s");
+ bp->v_jet_max *= (1e5 / (us->UnitLength_in_cgs / us->UnitTime_in_cgs));
+
} else if (strcmp(temp5, "MassLoading") == 0) {
bp->AGN_jet_velocity_model = AGN_jet_velocity_mass_loading;
@@ -960,9 +971,6 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
parser_get_param_float(params, "SPINJETAGN:opening_angle_in_degrees");
bp->opening_angle = bp->opening_angle * M_PI / 180.;
- bp->jet_h_r_slope =
- parser_get_param_float(params, "SPINJETAGN:jet_h_r_slope");
-
bp->eps_f_jet = parser_get_param_float(params, "SPINJETAGN:eps_f_jet");
if ((bp->eps_f_jet <= 0.) || (bp->eps_f_jet > 1.)) {
@@ -1000,14 +1008,6 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->fix_jet_direction);
}
- bp->accretion_efficiency =
- parser_get_param_float(params, "SPINJETAGN:accretion_efficiency");
-
- if (bp->accretion_efficiency <= 0.) {
- error("The accretion efficiency must be larger than 0., not %f",
- bp->accretion_efficiency);
- }
-
bp->fix_radiative_efficiency =
parser_get_param_int(params, "SPINJETAGN:fix_radiative_efficiency");
@@ -1029,6 +1029,26 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
bp->radiative_efficiency);
}
+ bp->use_ADIOS_winds =
+ parser_get_param_int(params, "SPINJETAGN:use_ADIOS_winds");
+
+ if ((bp->use_ADIOS_winds != 0) && (bp->use_ADIOS_winds != 1)) {
+ error(
+ "The use_ADIOS_winds parameter must be either 0 or 1, "
+ "not %d",
+ bp->use_ADIOS_winds);
+ }
+
+ bp->slim_disc_wind_factor =
+ parser_get_param_float(params, "SPINJETAGN:slim_disc_wind_factor");
+
+ if ((bp->slim_disc_wind_factor < 0) || (bp->slim_disc_wind_factor > 1)) {
+ error(
+ "The slim_disc_wind_factor parameter must be between 0 and 1, "
+ "(inclusive), not %f",
+ bp->slim_disc_wind_factor);
+ }
+
char temp6[60];
parser_get_param_string(params, "SPINJETAGN:AGN_jet_feedback_model", temp6);
if (strcmp(temp6, "MinimumDistance") == 0)
@@ -1044,6 +1064,47 @@ INLINE static void black_holes_props_init(struct black_holes_props *bp,
"The AGN feedback model must be MinimumDistance, MaximumDistance, "
"SpinAxis or MinimumDensity, not %s",
temp6);
+
+ bp->accretion_efficiency_slim =
+ parser_get_param_float(params, "SPINJETAGN:accretion_efficiency_slim");
+
+ if ((bp->accretion_efficiency_slim < 0) ||
+ (bp->accretion_efficiency_slim > 1)) {
+ error(
+ "The accretion_efficiency_slim parameter must be between 0 and 1, "
+ "(inclusive), not %f",
+ bp->accretion_efficiency_slim);
+ }
+
+ char temp7[60];
+ parser_get_param_string(params, "SPINJETAGN:accretion_efficiency_mode",
+ temp7);
+ if (strcmp(temp7, "Constant") == 0) {
+ bp->accretion_efficiency_mode = BH_accretion_efficiency_constant;
+ bp->accretion_efficiency_thick =
+ parser_get_param_float(params, "SPINJETAGN:accretion_efficiency_thick");
+
+ if ((bp->accretion_efficiency_thick < 0) ||
+ (bp->accretion_efficiency_thick > 1)) {
+ error(
+ "The accretion_efficiency_thick parameter must be between 0 and 1, "
+ "(inclusive), not %f",
+ bp->accretion_efficiency_thick);
+ }
+ } else if (strcmp(temp7, "Variable") == 0) {
+ bp->accretion_efficiency_mode = BH_accretion_efficiency_variable;
+ bp->ADIOS_s = parser_get_param_float(params, "SPINJETAGN:ADIOS_s");
+
+ if ((bp->ADIOS_s < 0) || (bp->ADIOS_s > 1)) {
+ error(
+ "The ADIOS_s parameter must be between 0 and 1, "
+ "(inclusive), not %f",
+ bp->ADIOS_s);
+ }
+ } else {
+ error("The accretion efficiency model must be Constant or Variable, not %s",
+ temp7);
+ }
}
/**
diff --git a/src/black_holes/SPIN_JET/black_holes_spin.h b/src/black_holes/SPIN_JET/black_holes_spin.h
index e9b7ba502745a0551be8667d8cb817547e5cc382..247c66c53d1c31446bfef1b4c979277ca7de701a 100644
--- a/src/black_holes/SPIN_JET/black_holes_spin.h
+++ b/src/black_holes/SPIN_JET/black_holes_spin.h
@@ -30,14 +30,41 @@
#include "inline.h"
#include "physical_constants.h"
+/**
+ * @brief Compute the gravitational radius of a black hole.
+ *
+ * @param a Black hole mass.
+ * @param constants Physical constants (in internal units).
+ */
+__attribute__((always_inline)) INLINE static float
+black_hole_gravitational_radius(float mass,
+ const struct phys_const* constants) {
+
+ const float r_G =
+ mass * constants->const_newton_G /
+ (constants->const_speed_light_c * constants->const_speed_light_c);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (r_G <= 0.f) {
+ error(
+ "Something went wrong with calculation of R_G of black holes. "
+ " R_G is %f instead of R_G > 0.",
+ r_G);
+ }
+#endif
+
+ return r_G;
+}
+
/**
* @brief Compute the radius of the horizon of a BH particle in gravitational
* units.
*
* @param a Black hole spin, -1 < a < 1.
*/
-__attribute__((always_inline)) INLINE static float r_hor(float a) {
- return 1. + sqrtf((1. - a) * (1. + a));
+__attribute__((always_inline)) INLINE static float black_hole_horizon_radius(
+ float a) {
+ return 1.f + sqrtf((1.f - a) * (1.f + a));
}
/**
@@ -49,32 +76,33 @@ __attribute__((always_inline)) INLINE static float r_hor(float a) {
*
* @param a Black hole spin, -1 < a < 1.
*/
-__attribute__((always_inline)) INLINE static float r_isco(float a) {
- const float Z1 = 1. + (cbrtf((1. + fabsf(a)) * (1. - a * a)) +
- cbrtf((1. - fabsf(a)) * (1. - a * a)));
- const float Z2 = sqrtf(3. * a * a + Z1 * Z1);
+__attribute__((always_inline)) INLINE static float black_hole_isco_radius(
+ float a) {
+ const float Z1 = 1.f + (cbrtf((1.f + fabsf(a)) * (1.f - a * a)) +
+ cbrtf((1.f - fabsf(a)) * (1.f - a * a)));
+ const float Z2 = sqrtf(3.f * a * a + Z1 * Z1);
const float R_ISCO =
- 3. + Z2 - a / fabsf(a) * sqrtf((3. - Z1) * (3. + Z1 + 2. * Z2));
+ 3. + Z2 - a / fabsf(a) * sqrtf((3.f - Z1) * (3.f + Z1 + 2.f * Z2));
#ifdef SWIFT_DEBUG_CHECKS
- if (Z1 > 3.) {
+ if (Z1 > 3.f) {
error(
"Something went wrong with calculation of Z1 factor for r_isco of"
" black holes. Z1 is %f instead of Z1 > 3.",
Z1);
}
- if ((3. + Z1 + 2. * Z2) < 0.) {
+ if ((3.f + Z1 + 2.f * Z2) < 0.f) {
error(
"Something went wrong with calculation of (3. + Z1 + 2. * Z2 ) "
"factor for r_isco of black holes. (3. + Z1 + 2. * Z2 ) is %f instead "
"of"
" (3. + Z1 + 2. * Z2 ) > 0.",
- 3. + Z1 + 2. * Z2);
+ 3.f + Z1 + 2.f * Z2);
}
- if (R_ISCO < 1.) {
+ if (R_ISCO < 1.f) {
error(
"Something went wrong with calculation of R_ISCO of black holes. "
"R_ISCO is %f instead >= 1.",
@@ -92,15 +120,16 @@ __attribute__((always_inline)) INLINE static float r_isco(float a) {
* @param a Black hole spin magnitude, 0 < a < 1.
* @param constants Physical constants (in internal units).
*/
-__attribute__((always_inline)) INLINE static float j_BH(
- struct bpart* bp, const struct phys_const* constants) {
+__attribute__((always_inline)) INLINE static float
+black_hole_angular_momentum_magnitude(struct bpart* bp,
+ const struct phys_const* constants) {
const float J_BH =
fabs(bp->subgrid_mass * bp->subgrid_mass * bp->spin *
constants->const_newton_G / constants->const_speed_light_c);
#ifdef SWIFT_DEBUG_CHECKS
- if (J_BH <= 0.) {
+ if (J_BH <= 0.f) {
error(
"Something went wrong with calculation of j_BH of black holes. "
" J_BH is %f instead of J_BH > 0.",
@@ -111,31 +140,6 @@ __attribute__((always_inline)) INLINE static float j_BH(
return J_BH;
}
-/**
- * @brief Compute the gravitational radius of a black hole.
- *
- * @param a Black hole mass.
- * @param constants Physical constants (in internal units).
- */
-__attribute__((always_inline)) INLINE static float R_gravitational(
- float mass, const struct phys_const* constants) {
-
- const float r_G =
- mass * constants->const_newton_G /
- (constants->const_speed_light_c * constants->const_speed_light_c);
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (r_G <= 0.) {
- error(
- "Something went wrong with calculation of R_G of black holes. "
- " R_G is %f instead of R_G > 0.",
- r_G);
- }
-#endif
-
- return r_G;
-}
-
/**
* @brief Compute the warp radius of a black hole particle.
*
@@ -159,27 +163,29 @@ __attribute__((always_inline)) INLINE static float R_gravitational(
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static float r_warp(
+__attribute__((always_inline)) INLINE static float black_hole_warp_radius(
struct bpart* bp, const struct phys_const* constants,
const struct black_holes_props* props) {
- /* Define placeholder variable for the result */
- float Rw = -1.;
+ /* Define placeholder value for the result. We will assign the final result
+ to this variable. */
+ float Rw = -1.f;
/* Gravitational radius */
- const float R_G = R_gravitational(bp->subgrid_mass, constants);
+ const float R_G =
+ black_hole_gravitational_radius(bp->subgrid_mass, constants);
/* Start branching depending on which accretion mode the BH is in */
if (bp->accretion_mode == BH_thick_disc) {
/* Eqn. 22 from Lubow et al. (2002) with H/R=h_0_ADAF (thick disk) */
- const float base = 15.36 * fabsf(bp->spin) / props->h_0_ADAF_2;
- Rw = R_G * powf(base, 0.4);
+ const float base = 15.36f * fabsf(bp->spin) / props->h_0_ADAF_2;
+ Rw = R_G * powf(base, 0.4f);
} else if (bp->accretion_mode == BH_slim_disc) {
/* Eqn. 22 from Lubow et al. (2002) with H/R=1/gamma_SD (slim disk) */
- const float base = 15.36 * fabsf(bp->spin) * props->gamma_SD;
- Rw = R_G * powf(base, 0.4);
+ const float base = 15.36f * fabsf(bp->spin) * props->gamma_SD;
+ Rw = R_G * powf(base, 0.4f);
} else if (bp->accretion_mode == BH_thin_disc) {
/* Start branching depending on which region of the thin disk we wish to
@@ -190,20 +196,20 @@ __attribute__((always_inline)) INLINE static float r_warp(
/* Calculate different factors in eqn. 11 (Griffin et al. 2019) for warp
radius of region b in Shakura & Sunyaev (1973) */
float mass_factor =
- powf(bp->subgrid_mass / (1e8 * constants->const_solar_mass), 0.2);
- float edd_factor = powf(bp->eddington_fraction, 0.4);
+ powf(bp->subgrid_mass / (1e8f * constants->const_solar_mass), 0.2f);
+ float edd_factor = powf(bp->eddington_fraction, 0.4f);
/* Gather the factors and finalize calculation */
const float base = mass_factor * fabsf(bp->spin) /
(props->xi_TD * props->alpha_factor_08 * edd_factor);
- const float rw = 3410. * 2. * R_G * powf(base, 0.625);
+ const float rw = 3410.f * 2.f * R_G * powf(base, 0.625f);
/* Self-gravity radius in region b: eqn. 16 in Griffin et al. */
- mass_factor =
- powf(bp->subgrid_mass / (1e8 * constants->const_solar_mass), -0.961);
- edd_factor = powf(bp->eddington_fraction, -0.353);
+ mass_factor = powf(
+ bp->subgrid_mass / (1e8f * constants->const_solar_mass), -0.961f);
+ edd_factor = powf(bp->eddington_fraction, -0.353f);
- const float rs = 4790. * 2. * R_G * mass_factor *
+ const float rs = 4790.f * 2.f * R_G * mass_factor *
props->alpha_factor_0549 * edd_factor;
/* Take the minimum */
@@ -214,20 +220,20 @@ __attribute__((always_inline)) INLINE static float r_warp(
/* Calculate different factors in eqn. A8 (Fiacconi et al. 2018) */
float mass_factor =
- powf(bp->subgrid_mass / (1e6 * constants->const_solar_mass), 0.2);
- float edd_factor = powf(bp->eddington_fraction, 0.3);
+ powf(bp->subgrid_mass / (1e6f * constants->const_solar_mass), 0.2f);
+ float edd_factor = powf(bp->eddington_fraction, 0.3f);
/* Gather the factors and finalize calculation */
const float base = mass_factor * fabsf(bp->spin) /
(props->xi_TD * props->alpha_factor_02 * edd_factor);
- const float rw = 1553. * 2. * R_G * powf(base, 0.5714);
+ const float rw = 1553.f * 2.f * R_G * powf(base, 0.5714f);
/* Repeat the same for self-gravity radius - eqn. A6 in F2018 */
- mass_factor =
- powf(bp->subgrid_mass / (1e6 * constants->const_solar_mass), -1.1556);
- edd_factor = powf(bp->eddington_fraction, -0.48889);
+ mass_factor = powf(
+ bp->subgrid_mass / (1e6f * constants->const_solar_mass), -1.1556f);
+ edd_factor = powf(bp->eddington_fraction, -0.48889f);
- const float rs = 1.2 * 100000. * 2. * R_G * mass_factor *
+ const float rs = 1.2f * 100000.f * 2.f * R_G * mass_factor *
props->alpha_factor_06222 * edd_factor;
/* Take the minimum */
@@ -236,7 +242,7 @@ __attribute__((always_inline)) INLINE static float r_warp(
}
#ifdef SWIFT_DEBUG_CHECKS
- if (Rw < 0.) {
+ if (Rw < 0.f) {
error(
"Something went wrong with calculation of Rw of black holes. "
" Rw is %f instead of Rw >= 0.",
@@ -272,15 +278,17 @@ __attribute__((always_inline)) INLINE static float r_warp(
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static double m_warp(
+__attribute__((always_inline)) INLINE static double black_hole_warp_mass(
struct bpart* bp, const struct phys_const* constants,
const struct black_holes_props* props) {
- /* Define placeholder variable for the result */
+ /* Define placeholder value for the result. We will assign the final result
+ to this variable. */
double Mw = -1.;
/* Gravitational radius */
- const float R_G = R_gravitational(bp->subgrid_mass, constants);
+ const float R_G =
+ black_hole_gravitational_radius(bp->subgrid_mass, constants);
/* Start branching depending on which accretion mode the BH is in */
if ((bp->accretion_mode == BH_thick_disc) ||
@@ -299,7 +307,7 @@ __attribute__((always_inline)) INLINE static double m_warp(
Mw = 2. * bp->accretion_rate /
(3. * props->alpha_acc * v_0 *
sqrtf(bp->subgrid_mass * constants->const_newton_G)) *
- powf(r_warp(bp, constants, props), 1.5);
+ powf(black_hole_warp_radius(bp, constants, props), 1.5);
} else {
/* Start branching depending on which region of the thin disk we wish to
@@ -312,7 +320,7 @@ __attribute__((always_inline)) INLINE static double m_warp(
powf(bp->subgrid_mass / (1e8 * constants->const_solar_mass), 2.2);
const float edd_factor = powf(bp->eddington_fraction, 0.6);
const float R_factor =
- powf(r_warp(bp, constants, props) / (2. * R_G), 1.4);
+ powf(black_hole_warp_radius(bp, constants, props) / (2. * R_G), 1.4);
/* Gather factors and finalize calculation */
Mw = constants->const_solar_mass * 1.35 * mass_factor *
@@ -325,7 +333,7 @@ __attribute__((always_inline)) INLINE static double m_warp(
powf(bp->subgrid_mass / (1e6 * constants->const_solar_mass), 2.2);
const float edd_factor = powf(bp->eddington_fraction, 0.7);
const float R_factor =
- powf(r_warp(bp, constants, props) / (2. * R_G), 1.25);
+ powf(black_hole_warp_radius(bp, constants, props) / (2. * R_G), 1.25);
Mw = constants->const_solar_mass * 0.01 * mass_factor *
props->alpha_factor_08_inv_10 * edd_factor * R_factor;
@@ -364,11 +372,13 @@ __attribute__((always_inline)) INLINE static double m_warp(
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static double j_warp(
- struct bpart* bp, const struct phys_const* constants,
- const struct black_holes_props* props) {
+__attribute__((always_inline)) INLINE static double
+black_hole_warp_angular_momentum(struct bpart* bp,
+ const struct phys_const* constants,
+ const struct black_holes_props* props) {
- /* Define placeholder variable for the result */
+ /* Define placeholder value for the result. We will assign the final result
+ to this variable. */
double Jw = -1.;
/* Start branching depending on which accretion mode the BH is in */
@@ -389,8 +399,9 @@ __attribute__((always_inline)) INLINE static double j_warp(
}
/* Gather factors for the final result */
+ const double r_warp = black_hole_warp_radius(bp, constants, props);
Jw = 2. * bp->accretion_rate * omega_0 / (2. * props->alpha_acc * v_0) *
- r_warp(bp, constants, props) * r_warp(bp, constants, props);
+ r_warp * r_warp;
} else {
/* Start branching depending on which region of the thin disk we wish to
@@ -402,14 +413,14 @@ __attribute__((always_inline)) INLINE static double j_warp(
For region b, c=-3/5 (see Griffin et al. 2019), and for region c, c=-3/4
(see Fiacconi et al. 2018). */
if (props->TD_region == TD_region_B) {
- Jw = 0.737 * m_warp(bp, constants, props) *
+ Jw = 0.737 * black_hole_warp_mass(bp, constants, props) *
sqrtf(bp->subgrid_mass * constants->const_newton_G *
- r_warp(bp, constants, props));
+ black_hole_warp_radius(bp, constants, props));
}
if (props->TD_region == TD_region_C) {
- Jw = 0.714 * m_warp(bp, constants, props) *
+ Jw = 0.714 * black_hole_warp_mass(bp, constants, props) *
sqrtf(bp->subgrid_mass * constants->const_newton_G *
- r_warp(bp, constants, props));
+ black_hole_warp_radius(bp, constants, props));
}
}
@@ -433,46 +444,64 @@ __attribute__((always_inline)) INLINE static double j_warp(
*
* @param a Black hole spin, -1 < a < 1.
*/
-__attribute__((always_inline)) INLINE static float eps_NT(float a) {
+__attribute__((always_inline)) INLINE static float eps_Novikov_Thorne(float a) {
#ifdef SWIFT_DEBUG_CHECKS
- if (r_isco(a) <= 0.6667) {
+ if (black_hole_isco_radius(a) <= 0.6667f) {
error(
- "Something went wrong with calculation of eps_NT of black holes. "
- " r_isco is %f instead of r_isco > 1.",
- r_isco(a));
+ "Something went wrong with calculation of eps_Novikov_Thorn of. "
+ "black holes. r_isco is %f instead of r_isco > 1.",
+ black_hole_isco_radius(a));
}
#endif
- return 1. - sqrtf(1. - 2. / 3. / r_isco(a));
+ return 1. - sqrtf(1. - 2.f / (3.f * black_hole_isco_radius(a)));
}
/**
* @brief Compute the spin- and accretion rate-dependant radiative efficiency
* of a BH particle in the super-Eddington (slim disk) regime.
*
- * This is eqn. 3 in Madau et al. (2014), which is based on numerical GR
+ * This is eqn. 3-6 in Madau et al. (2014), which is based on numerical GR
* results by Sadowski (2009).
*
* @param a Black hole spin, -1 < a < 1.
* @param m_dot Accretion rate normalized to the Eddington rate.
*/
-__attribute__((always_inline)) INLINE static float eps_SD(float a, float mdot) {
- const float B = powf(4.627 - 4.445 * a, -0.5524);
- const float C = powf(827.3 - 718.1 * a, -0.706);
- const float A = powf(0.9663 - 0.9292 * a, -0.5693);
+__attribute__((always_inline)) INLINE static float eps_slim_disc(float a,
+ float mdot) {
+ const float B = powf(4.627f - 4.445f * a, -0.5524f);
+ const float C = powf(827.3f - 718.1f * a, -0.706f);
+ const float A = powf(0.9663f - 0.9292f * a, -0.5693f);
#ifdef SWIFT_DEBUG_CHECKS
- if (mdot <= 0.) {
+ if (mdot <= 0.f) {
error(
- "The calculation of eps_SD was called even though mdot is %f. "
+ "The calculation of eps_slim_disc was called even though mdot is %f. "
" This function should not have been called if the accretion rate is "
" not > 0.",
mdot);
}
#endif
- return 0.1 / mdot * (0.985 / (B + 1.6 / mdot) + 0.015 / (C + 1.6 / mdot)) * A;
+ /* Since we use a definition of the Eddington ratio (mdot) that includes
+ the varying (Novikov-Thorne) radiative efficiency, we need to rescale this
+ back to a constant one, as the paper provides a formula assuming a
+ constant radiative efficiency. They use a value of 1/16, so we redefine
+ the Eddington ratio using the ratio of efficiencies. */
+ const float constant_rad_efficiency = 1.f / 16.f;
+
+ mdot = mdot * constant_rad_efficiency / eps_Novikov_Thorne(a);
+
+ /* Return radiative efficiency as given by Eqn 3 from Madau et al. (2014).
+ Note that the equation provided in the paper is for L / L_Edd, rather than
+ for L / (f_Edd * M_Edd * c^2). We thus need to multiply their Eqn 3 by
+ L_Edd / (f_Edd * M_Edd * c^2) = eps_rad_constant / mdot. Here we have used
+ M_Edd = L_Edd / (eps_rad_constant * c^2). Also note that mdot = f_Edd in
+ our notation. */
+
+ return (constant_rad_efficiency / mdot) *
+ (0.985f / (B + mdot) + 0.015f / (C + mdot)) * A;
}
/**
@@ -480,18 +509,19 @@ __attribute__((always_inline)) INLINE static float eps_SD(float a, float mdot) {
*
* The possible modes are the thick disk, thin disk and slim disk, in
* order of increasing accretion rate. The transition from thick to thin disk
- * is at 0.4*alpha^2, based on current theory (Yuan & Narayan 2014). The
- * transition from thin to slim disk occurs when the slim disk efficiency
- * becomes sufficiently weak compared to the thin disk one. We parametrize
- * this transition as occuring at eps_SD = props->TD_SD_eps_r_threshold *
- * eps_TD, with props->TD_SD_eps_r_threshold < 1 of order 0.5
+ * is currently governed by a free parameter, props->mdot_crit_ADAF (of order
+ * 0.01. The transition between the thin and slim disc is assumed to take place
+ * at mdot = 1, i.e. for super-Eddington accretion. Note that this assumption
+ * only works if we define mdot by using the spin-dependent radiative
+ * efficiency, which we do.
*
* @param bp Pointer to the b-particle data.
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static void decide_mode(
- struct bpart* bp, const struct black_holes_props* props) {
+__attribute__((always_inline)) INLINE static void
+black_hole_select_accretion_mode(struct bpart* bp,
+ const struct black_holes_props* props) {
/* For deciding the accretion mode, we want to use the Eddington fraction
* calculated using the raw, unsuppressed accretion rate. This means that
@@ -499,22 +529,15 @@ __attribute__((always_inline)) INLINE static void decide_mode(
* already suppressed, needs to be unsuppressed (increased) to retrieve the
* raw Bondi-based Eddington ratio. */
float eddington_fraction_Bondi = bp->eddington_fraction;
- if (bp->accretion_mode == BH_thick_disc)
- eddington_fraction_Bondi *= 1. / props->accretion_efficiency;
+ eddington_fraction_Bondi *= 1.f / bp->accretion_efficiency;
if (eddington_fraction_Bondi < props->mdot_crit_ADAF) {
bp->accretion_mode = BH_thick_disc;
} else {
/* The disc is assumed to be slim (super-Eddington) if the Eddington
- * fraction is above 1. However, the Eddington fraction is typically
- * (also in this code and model) defined using an Eddington luminosity
- * with an assumed radiative efficiency of eps_r = 0.1. For this purpose
- * it is more sensible to define the Eddington luminosity using the spin-
- * dependent eps_NT(a). We thus assume that the disc becomes super-
- * Eddington at a critical f_Edd(a) = f_Edd * (eps_NT(a)/0.1) = 1. */
- if ((eddington_fraction_Bondi * eps_NT(bp->spin) / 0.1 > 1.) &&
- (props->include_slim_disk)) {
+ * fraction is above 1. */
+ if ((eddington_fraction_Bondi > 1.f) && (props->include_slim_disk)) {
bp->accretion_mode = BH_slim_disc;
} else {
bp->accretion_mode = BH_thin_disc;
@@ -522,142 +545,182 @@ __attribute__((always_inline)) INLINE static void decide_mode(
}
/* If we do not include radiative feedback, then we force the disk to be in
- the thick disk mode */
+ the thick disk mode always. */
if (props->turn_off_radiative_feedback) {
bp->accretion_mode = BH_thick_disc;
}
/* similar for if we do not include jets - we force the disk to be thin */
- if (props->include_jets == 0) {
+ if (!props->include_jets) {
bp->accretion_mode = BH_thin_disc;
}
}
/**
- * @brief Compute the aspect ratio of the subgrid accretion disk.
+ * @brief Compute the accretion efficiency of a BH particle.
*
* The result depends on bp->accretion_mode (thick disk, thin disk or
- * slim disk). For the thick disk and slim disk, the aspect ratio is
- * a constant, H/R = h_0.
- *
- * For the thin disk, the result depends on props->TD_region (B - region b from
- * Shakura & Sunyaev 1973, C - region c from Shakura & Sunyaev 1973). In region
- * b, we take H/R as eqn. 8 in Griffin et al. (2019), and in region c H/r is
- * taken directly as eqn. 2.19 from Shakura & Sunyaev (1973).
+ * slim disk). We assume no accretion efficiency (100%) in the thin disk,
+ * and allow for options for a non-zero accretion efficiency in the thick
+ * and slim disc. For both we allow the option of constant values, and for the
+ * thick disc we allow an option for a scaling with Eddington ratio that is
+ * motivated by simulations.
*
* @param bp Pointer to the b-particle data.
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static float aspect_ratio(
- struct bpart* bp, const struct phys_const* constants,
- const struct black_holes_props* props) {
-
- /* Define placeholder variable for the result */
- float h_0 = -1.;
-
- /* Start branching depending on which accretion mode the BH is in.
- * We assume the thick and slim disc to be very similar in geometrical
- * terms. */
- if ((bp->accretion_mode == BH_thick_disc) ||
- (bp->accretion_mode == BH_slim_disc)) {
- h_0 = props->h_0_ADAF;
- } else {
-
- /* Start branching depending on which region of the thin disk we wish to
- base the model upon (TD_region=B: region b from Shakura & Sunyaev 1973,
- or TD_region=C: region c). */
- if (props->TD_region == TD_region_B) {
+__attribute__((always_inline)) INLINE static float
+black_hole_accretion_efficiency(struct bpart* bp,
+ const struct black_holes_props* props,
+ const struct phys_const* constants,
+ const struct cosmology* cosmo) {
+
+ if (bp->accretion_mode == BH_thick_disc ||
+ bp->accretion_mode == BH_slim_disc) {
+
+ if (props->accretion_efficiency_mode == BH_accretion_efficiency_constant) {
+ if (bp->accretion_mode == BH_thick_disc) {
+ return props->accretion_efficiency_thick;
+ } else if (bp->accretion_mode == BH_slim_disc) {
+ return props->accretion_efficiency_slim;
+ } else {
- /* Compute factors for eqn. 8 in Griffin et al. (2019). */
- const float mass_factor =
- powf(bp->subgrid_mass / (1e8 * constants->const_solar_mass), -0.1);
- const float edd_factor = powf(bp->eddington_fraction, 0.2);
- const float R_G = R_gravitational(bp->subgrid_mass, constants);
- const float R_factor =
- powf(r_warp(bp, constants, props) / (2. * R_G), 0.05);
+#ifdef SWIFT_DEBUG_CHECKS
+ error(
+ "This branch of the function accretion_efficiency() should not"
+ " have been reached!");
+#endif
- /* Gather factors and finalize calculation. */
- h_0 = 1.25 * 0.001 * mass_factor * props->alpha_factor_01 * edd_factor *
- R_factor;
- }
- if (props->TD_region == TD_region_C) {
+ return 1.f;
+ }
+ } else if (props->accretion_efficiency_mode ==
+ BH_accretion_efficiency_variable) {
+
+ if (bp->accretion_mode == BH_thick_disc) {
+
+ /* Compute the transition radius between an outer thin disc and an
+ * inner thick disc. This is assumed to happen at 10 R_G at the
+ * critical value of the Eddington ratio between the two regimes.
+ * The transition radius then increases as 1 / f_Edd^2. Note that
+ * we also need to use the raw (unsuppressed) Eddington ratio here,
+ * hence the multiplication by accretion efficiencies. Note that the
+ * units of the transition radius here are in R_G. */
+ float R_tr = 10.f * props->mdot_crit_ADAF * props->mdot_crit_ADAF *
+ bp->accretion_efficiency * bp->accretion_efficiency /
+ (bp->eddington_fraction * bp->eddington_fraction);
+
+ /* We need to also compute the Bondi radius (in units of R_G), which
+ * can be expressed in terms of the ratio between speed of light and
+ * sound speed. */
+ const double c = constants->const_speed_light_c;
+ float gas_c_phys2 = bp->sound_speed_gas * cosmo->a_factor_sound_speed;
+ gas_c_phys2 = gas_c_phys2 * gas_c_phys2;
+ float R_B = c * c / gas_c_phys2;
+
+ /* Limit the transition radius to no larger than R_B and no smaller
+ * than 10 R_G. */
+ R_tr = fminf(R_B, R_tr);
+ R_tr = fmaxf(10.f, R_tr);
+
+ /* Implement the actual scaling of accretion efficiency with transition
+ * radius as found by GRMHD simulations. */
+ float suppr_factor = powf(10.f / R_tr, props->ADIOS_s);
+ return suppr_factor;
+ } else if (bp->accretion_mode == BH_slim_disc) {
+ return props->accretion_efficiency_slim;
+ } else {
- /* Compute factors for eqn. 2.19 in Shakura & Sunyaev (1973). */
- const float mass_factor =
- powf(bp->subgrid_mass / (1e8 * constants->const_solar_mass), -0.1);
- const float edd_factor = powf(bp->eddington_fraction, 0.15);
- const float R_G = R_gravitational(bp->subgrid_mass, constants);
- const float R_factor =
- powf(r_warp(bp, constants, props) / (2. * R_G), 0.125);
+#ifdef SWIFT_DEBUG_CHECKS
+ error(
+ "This branch of the function accretion_efficiency() should not"
+ " have been reached!");
+#endif
- /* Gather factors and finalize calculation. */
- h_0 = 1.15 * 0.001 * mass_factor * props->alpha_factor_01 * edd_factor *
- R_factor;
- }
- }
+ return 1.f;
+ }
+ } else {
#ifdef SWIFT_DEBUG_CHECKS
- if (h_0 <= 0.) {
- error(
- "Something went wrong with calculation of h_0 of black holes. "
- " h_0 is %f instead of h_0 > 0.",
- h_0);
- }
+ error(
+ "This branch of the function accretion_efficiency() should not"
+ " have been reached!");
#endif
- return h_0;
+ return 1.f;
+ }
+ } else {
+ return 1.f;
+ }
}
/**
* @brief Compute the jet efficiency of a BH particle.
*
* The result depends on bp->accretion_mode (thick disk, thin disk or
- * slim disk), through the varying H/R aspect ratios.
+ * slim disk).
*
* The equation implemented is eqn. 9 from Tchekhovskoy et al. (2010), with the
- * dimensionless magnetic flux phi taken as eqn. 9 from Narayan et al. (2021).
- *
- * The dependence on the aspect ratio comes from results in Tchekhovskoy et al.
- * (2014) and the dependence in classical Blandford & Znajek (1979) jet theory.
+ * dimensionless magnetic flux phi taken as eqn. 9 from Narayan et al. (2022),
+ * and an additional modification from Ricarte et al. (2023).
*
* @param bp Pointer to the b-particle data.
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static float jet_efficiency(
+__attribute__((always_inline)) INLINE static float black_hole_jet_efficiency(
struct bpart* bp, const struct black_holes_props* props) {
- float jet_eff = -1.;
+ /* Define placeholder value for the result. We will assign the final result
+ to this variable. */
+ float jet_eff = -1.f;
if (props->fix_jet_efficiency) {
jet_eff = props->jet_efficiency;
} else {
- const float kappa = 0.05;
+
+ /* Numerical prefactor that appears in the jet power formula, related to
+ the geometry of the magnetic field. */
+ const float kappa = 0.05f;
+
+ /* Definition of angular velocity at the BH event horizon */
const float horizon_ang_vel =
- bp->spin / (2. * (1. + sqrtf(1. - bp->spin * bp->spin)));
- const float phi = -20.2 * bp->spin * bp->spin * bp->spin -
- 14.9 * bp->spin * bp->spin + 34. * bp->spin + 52.6;
- jet_eff = kappa * 0.25 * M_1_PI * phi * phi *
- powf(bp->aspect_ratio * 3.333, props->jet_h_r_slope) *
- horizon_ang_vel * horizon_ang_vel *
- (1. + 1.38 * horizon_ang_vel * horizon_ang_vel -
- 9.2 * horizon_ang_vel * horizon_ang_vel * horizon_ang_vel *
+ bp->spin / (2.f * (1.f + sqrtf(1.f - bp->spin * bp->spin)));
+
+ /* Dimensionless magnetic flux as a function of BH spin, using Eqn. (15)
+ from Narayan et al. (2022). */
+ float phi = -20.2f * bp->spin * bp->spin * bp->spin -
+ 14.9f * bp->spin * bp->spin + 34.f * bp->spin + 52.6f;
+
+ float Eddington_ratio = bp->eddington_fraction;
+
+ /* Suppress the magnetic flux if we are in the thin or slim disc,
+ * according to results from Ricarte et al. (2023). */
+ if ((bp->accretion_mode == BH_slim_disc) ||
+ (props->use_jets_in_thin_disc && bp->accretion_mode == BH_thin_disc)) {
+ phi *= powf(Eddington_ratio / 1.88f, 1.29f) /
+ (1.f + powf(Eddington_ratio / 1.88f, 1.29f));
+ }
+
+ /* Full jet efficiency formula as in Tchekhovskoy et al. (2010). */
+ jet_eff = kappa * 0.25f * M_1_PI * phi * phi * horizon_ang_vel *
+ horizon_ang_vel *
+ (1.f + 1.38f * horizon_ang_vel * horizon_ang_vel -
+ 9.2f * horizon_ang_vel * horizon_ang_vel * horizon_ang_vel *
horizon_ang_vel);
}
/* Turn off jet feedback if we want to do that */
- if (props->include_jets == 0) {
- jet_eff = 0.;
+ if (!props->include_jets) {
+ jet_eff = 0.f;
}
/* Turn off jets in thin disk mode if we want to do that */
- if ((props->turn_off_secondary_feedback) &&
- (bp->accretion_mode == BH_thin_disc)) {
- jet_eff = 0.;
+ if ((bp->accretion_mode == BH_thin_disc) && (!props->use_jets_in_thin_disc)) {
+ jet_eff = 0.f;
}
#ifdef SWIFT_DEBUG_CHECKS
- if (jet_eff < 0.) {
+ if (jet_eff < 0.f) {
error(
"Something went wrong with calculation of jet efficiency of black "
"holes. jet_eff is %f instead of jet_eff >= 0.",
@@ -684,14 +747,16 @@ __attribute__((always_inline)) INLINE static float jet_efficiency(
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static float rad_efficiency(
- struct bpart* bp, const struct black_holes_props* props) {
+__attribute__((always_inline)) INLINE static float
+black_hole_radiative_efficiency(struct bpart* bp,
+ const struct black_holes_props* props) {
/* Calculate Novikov-Thorne efficiency, which will be needed twice. */
- const float eps_TD = eps_NT(bp->spin);
+ const float eps_TD = eps_Novikov_Thorne(bp->spin);
- /* Define placeholder variable for the result */
- float rad_eff = -1;
+ /* Define placeholder value for the result. We will assign the final result
+ to this variable. */
+ float rad_eff = -1.f;
if (props->fix_radiative_efficiency) {
rad_eff = props->radiative_efficiency;
@@ -705,11 +770,11 @@ __attribute__((always_inline)) INLINE static float rad_efficiency(
} else if (bp->accretion_mode == BH_slim_disc) {
/* Assign Madau 2014 efficiency to the slim disk. */
- rad_eff = eps_SD(bp->spin, bp->eddington_fraction);
+ rad_eff = eps_slim_disc(bp->spin, bp->eddington_fraction);
} else {
#ifdef SWIFT_DEBUG_CHECKS
- if (props->beta_acc > 1.) {
+ if (props->beta_acc > 1.f) {
error(
"Something went wrong with calculation of radiative efficiency of "
" black holes. beta_acc is %f instead of beta_acc < 1.",
@@ -717,30 +782,36 @@ __attribute__((always_inline)) INLINE static float rad_efficiency(
}
#endif
- /* Assign Mahadevan 1997 efficiency to the thick disk. */
- if (bp->eddington_fraction < props->edd_crit_thick) {
- rad_eff = 4.8 * eps_TD / r_isco(bp->spin) * (1. - props->beta_acc) *
- props->delta_ADAF;
+ /* Assign Mahadevan 1997 efficiency to the thick disk. We implement these
+ using Eqns. (29) and (30) from Griffin et al. (2019). */
+ if (bp->eddington_fraction < props->mdot_crit_ADAF) {
+ rad_eff = 4.8f * eps_TD / black_hole_isco_radius(bp->spin) *
+ (1.f - props->beta_acc) * props->delta_ADAF;
} else {
- rad_eff = 2.4 * eps_TD / r_isco(bp->spin) * props->beta_acc *
- bp->eddington_fraction * props->alpha_acc_2_inv;
+ rad_eff = 2.4f * eps_TD / black_hole_isco_radius(bp->spin) *
+ props->beta_acc * bp->eddington_fraction *
+ props->alpha_acc_2_inv;
+ }
+
+ /* Add contribution of truncated thin disc from larger radii */
+ if (props->accretion_efficiency_mode ==
+ BH_accretion_efficiency_variable) {
+ float R_tr = 10.f * props->mdot_crit_ADAF * props->mdot_crit_ADAF *
+ bp->accretion_efficiency * bp->accretion_efficiency /
+ (bp->eddington_fraction * bp->eddington_fraction);
+ R_tr = fmaxf(10.f, R_tr);
+ rad_eff += 1.f - sqrtf(1. - 2.f / (3.f * R_tr));
}
}
}
/* Turn off radiative feedback if we want to do that */
if (props->turn_off_radiative_feedback) {
- rad_eff = 0.;
- }
-
- /* Turn off radiation in the thick disk mode if we want to do that */
- if ((props->turn_off_secondary_feedback) &&
- (bp->accretion_mode == BH_thick_disc)) {
- rad_eff = 0.;
+ rad_eff = 0.f;
}
#ifdef SWIFT_DEBUG_CHECKS
- if (rad_eff < 0.) {
+ if (rad_eff < 0.f) {
error(
"Something went wrong with calculation of radiative efficiency of "
" black holes. rad_eff is %f instead of rad_eff >= 0.",
@@ -751,6 +822,67 @@ __attribute__((always_inline)) INLINE static float rad_efficiency(
return rad_eff;
}
+/**
+ * @brief Compute the wind efficiency of a BH particle.
+ *
+ * The result depends on bp->accretion_mode (thick disk, thin disk or
+ * slim disk), with no wind assumed for the thin disc (effectively, the
+ * radiation launches its own wind, while in the thick/slim disc, it is gas
+ * pressure/MHD effects that launch the wind. In all cases, the wind is dumped
+ * as thermal energy, alongside radiation.
+ *
+ * For the thick disk, we take the results from Sadowski et al. (2013)
+ * (2013MNRAS.436.3856S), which is applicable to MAD discs. For the slim disc,
+ * we constructed a fitting function by using the total MHD efficiency from
+ * Ricarte et al. (2023) (2023ApJ...954L..22R), which includes both winds and
+ * jets, and subtracting from that the jet efficiency used by our model.
+ *
+ * @param bp Pointer to the b-particle data.
+ * @param constants Physical constants (in internal units).
+ * @param props Properties of the black hole scheme.
+ */
+__attribute__((always_inline)) INLINE static float black_hole_wind_efficiency(
+ struct bpart* bp, const struct black_holes_props* props) {
+
+ /* (Dimensionless) magnetic flux on the BH horizon, as given by the
+ Narayan et al. (2022) fitting function for MAD discs. */
+ float phi = -20.2f * bp->spin * bp->spin * bp->spin -
+ 14.9f * bp->spin * bp->spin + 34.f * bp->spin + 52.6f;
+
+ if (bp->accretion_mode == BH_slim_disc) {
+
+ /* We need to suppress the magnetic flux by an Eddington-ratio-dependent
+ factor (Equation 3 from Ricarte et al. 2023). */
+ float Eddington_ratio = bp->eddington_fraction;
+ phi *= powf(Eddington_ratio / 1.88f, 1.29f) /
+ (1.f + powf(Eddington_ratio / 1.88f, 1.29f));
+ float phi_factor = (1.f + (phi / 50.f) * (phi / 50.f));
+
+ float horizon_ang_vel =
+ bp->spin / (2.f * (1.f + sqrtf(1.f - bp->spin * bp->spin)));
+ float spin_factor =
+ 1.f - 8.f * horizon_ang_vel * horizon_ang_vel + 1.f * horizon_ang_vel;
+ if (bp->spin > 0.f) {
+ spin_factor = max(0.4f, spin_factor);
+ } else {
+ spin_factor = max(0.f, spin_factor);
+ }
+
+ /* Final result for slim disc wind efficiency. (Not published
+ yet anywhere) */
+ return props->slim_disc_wind_factor * 0.0635f * phi_factor * spin_factor;
+ } else if (bp->accretion_mode == BH_thick_disc && props->use_ADIOS_winds) {
+
+ /* Equation (29) from Sadowski et al. (2013). */
+ float horizon_ang_vel =
+ bp->spin / (2.f * (1.f + sqrtf(1.f - bp->spin * bp->spin)));
+ return 0.005f * (1.f + 3.f * phi * phi / 2500.f * horizon_ang_vel *
+ horizon_ang_vel / 0.04f);
+ } else {
+ return 0.f;
+ }
+}
+
/**
* @brief Compute the spec. ang. mom. at the inner radius of a BH particle.
*
@@ -771,31 +903,34 @@ __attribute__((always_inline)) INLINE static float l_acc(
struct bpart* bp, const struct phys_const* constants,
const struct black_holes_props* props) {
- /* Define placeholder variable for the result */
- float L = -1.;
+ /* Define placeholder value for the result. We will assign the final result
+ to this variable. */
+ float L = -1.f;
#ifdef SWIFT_DEBUG_CHECKS
- if (r_isco(bp->spin) <= 0.6667) {
+ if (black_hole_isco_radius(bp->spin) <= 0.6667f) {
error(
"Something went wrong with calculation of l_acc of black holes. "
" r_isco is %f instead of r_isco > 1.",
- r_isco(bp->spin));
+ black_hole_isco_radius(bp->spin));
}
#endif
/* Spec. ang. mom. at ISCO */
- const float L_ISCO = 0.385 * (1. + 2. * sqrtf(3. * r_isco(bp->spin) - 2.));
+ const float L_ISCO =
+ 0.385f *
+ (1.f + 2.f * sqrtf(3.f * black_hole_isco_radius(bp->spin) - 2.f));
/* Branch depending on which accretion mode the BH is in */
if ((bp->accretion_mode == BH_thick_disc) ||
(bp->accretion_mode == BH_slim_disc)) {
- L = 0.45 * L_ISCO;
+ L = 0.45f * L_ISCO;
} else {
L = L_ISCO;
}
#ifdef SWIFT_DEBUG_CHECKS
- if (L <= 0.) {
+ if (L <= 0.f) {
error(
"Something went wrong with calculation of l_acc of black holes. "
" l_acc is %f instead of l_acc > 0.",
@@ -807,51 +942,85 @@ __attribute__((always_inline)) INLINE static float l_acc(
}
/**
- * @brief Compute the evolution of the spin of a BH particle.
+ * @brief Compute the evolution of the spin of a BH particle. This
+ * spinup/spindown rate is equal to da / dln(M_BH)_0, or
+ * da / (d(M_BH,0)/M_BH ), where the subscript '0' means that it is
+ * the mass increment before losses due to jets, radiation or winds
+ * (i.e. without the effect of efficiencies).
*
* The result depends on bp->accretion_mode (thick disk, thin disk or
* slim disk), due to differing spec. ang. momenta as well as jet and
* radiative efficiencies.
*
- * This equation corresponds to eqn. 2 in Benson & Babul (2009), including
- * a jet spindown term.
+ * For the thick disc, we use the jet spindown formula from Narayan et al.
+ * (2022). For the slim and thin disc, we use the formula from Ricarte et al.
+ * (2023).
*
* @param bp Pointer to the b-particle data.
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static float da_dln_mbh_0(
+__attribute__((always_inline)) INLINE static float black_hole_spinup_rate(
struct bpart* bp, const struct phys_const* constants,
const struct black_holes_props* props) {
const float a = bp->spin;
- if ((a == 0.) || (a < -0.9981) || (a > 0.9981)) {
+ if ((a == 0.f) || (a < -0.9981f) || (a > 0.9981f)) {
error(
- "The da_dln_mbh_0 function was called and spin is %f. Spin should "
+ "The spinup function was called and spin is %f. Spin should "
" not be a = 0, a < -0.998 or a > 0.998.",
a);
}
- float spinup_rate = 0.;
-
- if (props->include_GRMHD_spindown) {
- if (bp->accretion_mode == BH_thin_disc) {
- spinup_rate = l_acc(bp, constants, props) -
- 2. * a * (1. - rad_efficiency(bp, props));
- } else {
- spinup_rate = 0.45 - 12.53 * a - 7.8 * a * a + 9.44 * a * a * a +
- 5.71 * a * a * a * a - 4.03 * a * a * a * a * a;
+ if (bp->accretion_mode == BH_thin_disc && !props->use_jets_in_thin_disc) {
+
+ /* If we are in the thin disc and use no jets, we use the simple spinup /
+ * spindown formula, e.g. from Benson & Babul (2009). This accounts for
+ * accretion only. */
+ return l_acc(bp, constants, props) -
+ 2.f * a * (1.f - bp->radiative_efficiency);
+
+ } else if (bp->accretion_mode == BH_thick_disc) {
+
+ /* Fiting function from Narayan et al. (2022) */
+ return 0.45f - 12.53f * a - 7.8f * a * a + 9.44f * a * a * a +
+ 5.71f * a * a * a * a - 4.03f * a * a * a * a * a;
+
+ } else if (bp->accretion_mode == BH_slim_disc ||
+ (bp->accretion_mode == BH_thin_disc &&
+ props->use_jets_in_thin_disc)) {
+
+ /* Fitting function from Ricarte et al. (2023). */
+ float Eddington_ratio =
+ bp->eddington_fraction * eps_Novikov_Thorne(bp->spin) / 0.1f;
+ float xi = Eddington_ratio * 0.017f;
+ float s_min = 0.86f - 1.94f * bp->spin;
+ float L_ISCO =
+ 0.385f *
+ (1.f + 2.f * sqrtf(3.f * black_hole_isco_radius(bp->spin) - 2.f));
+ float s_thin = L_ISCO - 2.f * a * (1.f - eps_Novikov_Thorne(bp->spin));
+ float s_HD = (s_thin + s_min * xi) / (1.f + xi);
+
+ float horizon_ang_vel =
+ fabsf(bp->spin) / (2.f * (1.f + sqrtf(1.f - bp->spin * bp->spin)));
+ float k_EM = 0.23f;
+ if (bp->spin > 0.f) {
+ k_EM = min(0.1f + 0.5f * bp->spin, 0.35f);
}
+
+ float s_EM = -1.f * bp->spin / fabsf(bp->spin) * bp->jet_efficiency *
+ (1.f / (k_EM * horizon_ang_vel) - 2.f * bp->spin);
+ return s_HD + s_EM;
} else {
- spinup_rate =
- l_acc(bp, constants, props) -
- 2. * a * (1. - rad_efficiency(bp, props)) -
- sqrtf(1. - a * a) / a *
- (a * a + (1. + sqrtf(1. - a * a)) * (1. + sqrtf(1. - a * a))) *
- jet_efficiency(bp, props);
- }
- return spinup_rate;
+#ifdef SWIFT_DEBUG_CHECKS
+ error(
+ "We shouldn't have reached this branch of the "
+ "black_hole_spinup_rate() function!");
+#endif
+
+ return 0;
+ }
}
/**
@@ -866,7 +1035,7 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_delta_T(
const struct bpart* bp, const struct black_holes_props* props,
const struct cosmology* cosmo, const struct phys_const* constants) {
- float delta_T = -1.;
+ float delta_T = -1.f;
if (props->AGN_heating_temperature_model ==
AGN_heating_temperature_constant) {
delta_T = props->AGN_delta_T_desired;
@@ -898,7 +1067,7 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_delta_T(
/* Calculate heating temperature from the power, smoothing length (proper,
not comoving), neighbour sound speed and neighbour mass. Apply floor. */
const float delta_T_repl =
- (2.f * 0.6 * constants->const_proton_mass * feedback_power *
+ (2.f * 0.6f * constants->const_proton_mass * feedback_power *
replenishment_time_scale) /
(3.f * constants->const_boltzmann_k * bp->ngb_mass);
@@ -911,7 +1080,7 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_delta_T(
(0.6 * constants->const_proton_mass) / (constants->const_boltzmann_k) *
powf(2.f * bp->h * cosmo->a * feedback_power /
(sqrtf(15.f) * bp->ngb_mass / ((double)bp->num_ngbs)),
- 0.6667);
+ 0.6667f);
/* Calculate minimum temperature from Dalla Vecchia & Schaye (2012) to
prevent numerical overcooling. This is in Kelvin. */
@@ -952,28 +1121,13 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_dv_jet(
float v_jet = -1.;
if (props->AGN_jet_velocity_model == AGN_jet_velocity_BH_mass) {
- /* Assign the halo mass according to an empirical relation given in the
- parameter file */
- const float halo_mass =
- powf(bp->subgrid_mass / props->v_jet_BH_mass_scaling_reference_mass,
+ v_jet =
+ powf((bp->subgrid_mass / props->v_jet_BH_mass_scaling_reference_mass),
props->v_jet_BH_mass_scaling_slope);
- /* Get the critical density and virial overdensity at this redshift */
- const float critical_density = cosmo->critical_density;
- const float overdensity = cosmo->overdensity_BN98;
-
- /* Gather the previous factors and compute the virial radius, virial
- velocity and finally the sound speed in the hot gas */
- const float virial_radius =
- cbrtf(3. * halo_mass / (4. * M_PI * overdensity * critical_density));
- const float virial_velocity =
- sqrtf(halo_mass * constants->const_newton_G / virial_radius);
- const float sound_speed = sqrtf(5. / 3. * 0.5) * virial_velocity;
-
- /* Return the jet velocity as some factor times the sound speed, apply
- * floor and ceiling values. */
- v_jet = fmaxf(props->v_jet_min, props->v_jet_cs_ratio * sound_speed);
- v_jet = fminf(props->v_jet_max, v_jet);
+ /* Apply floor and ceiling values */
+ v_jet = props->v_jet_max * fminf(v_jet, 1.f);
+ v_jet = fmaxf(v_jet, props->v_jet_min);
} else if (props->AGN_jet_velocity_model == AGN_jet_velocity_constant) {
v_jet = props->v_jet;
@@ -1016,7 +1170,7 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_dv_jet(
* power, smoothing length (proper, not comoving), neighbour sound speed
* and (total) neighbour mass. */
const float v_jet_repl =
- sqrtf(jet_power * replenishment_time_scale / (2. * bp->ngb_mass));
+ sqrtf(jet_power * replenishment_time_scale / (2.f * bp->ngb_mass));
/* Calculate jet velocity from the crossing condition, i.e. set the
* velocity such that a new particle pair will be launched roughly when
@@ -1024,7 +1178,7 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_dv_jet(
* power, smoothing length and neighbour mass (per particle, not total). */
const float v_jet_cross =
cbrtf(bp->h * cosmo->a * jet_power /
- (4. * bp->ngb_mass / ((double)bp->num_ngbs)));
+ (4.f * bp->ngb_mass / ((double)bp->num_ngbs)));
/* Find whichever of these two is the minimum, and multiply it by an
* arbitrary scaling factor (whose fiducial value is 1, i.e. no
@@ -1041,7 +1195,7 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_dv_jet(
"supported.");
}
- if (v_jet <= 0.) {
+ if (v_jet <= 0.f) {
error(
"The black_hole_feedback_dv_jet returned a value less than 0. which "
" is v_jet = %f.",
@@ -1069,7 +1223,7 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_dv_jet(
* @param cos_gamma cosine of the angle between the second spin and the initial
* total angular momentu
*/
-__attribute__((always_inline)) INLINE static float l_variable(
+__attribute__((always_inline)) INLINE static float black_hole_l_variable(
const float a1, const float a2, const float q, const float eta,
const float cos_alpha, const float cos_beta, const float cos_gamma) {
@@ -1095,38 +1249,6 @@ __attribute__((always_inline)) INLINE static float l_variable(
return term1 + term2 + term3 + term4 + term5;
}
-/**
- * @brief Auxilliary function used for the calculation of final spin of
- * a BH merger.
- *
- * This implements the fitting formula for the final spin from Barausse &
- * Rezolla (2009), ApJ, 704, Equation 6. It is used in the merger_spin_evolve()
- * function.
- *
- * @param a1 spin of the first (more massive) black hole
- * @param a2 spin of the less massive black hole
- * @param q mass ratio of the two black holes, 0 < q < 1
- * @param cos_alpha cosine of the angle between the two spins
- * @param cos_beta cosine of the angle between the first spin and the initial
- * total angular momentum
- * @param cos_gamma cosine of the angle between the second spin and the initial
- * total angular momentu
- */
-__attribute__((always_inline)) INLINE static float final_spin(
- const float a1, const float a2, const float q, const float cos_alpha,
- const float cos_beta, const float cos_gamma, const float l) {
-
- /* Gather the terms of Eqn. 6 */
- const float term1 = a1 * a1;
- const float term2 = a2 * a2 * q * q * q * q;
- const float term3 = 2.f * a1 * a2 * q * q * cos_alpha;
- const float term4 = 2.f * (a1 * cos_beta + a2 * q * q * cos_gamma) * l * q;
- const float term5 = l * l * q * q;
-
- /* Calculate the final spin */
- return sqrtf(term1 + term2 + term3 + term4 + term5) / ((1.f + q) * (1.f + q));
-}
-
/**
* @brief Auxilliary function used for the calculation of mass lost to GWs.
*
@@ -1159,14 +1281,14 @@ __attribute__((always_inline)) INLINE static float mass_fraction_lost_to_GWs(
* @param constants Physical constants (in internal units).
* @param props Properties of the black hole scheme.
*/
-__attribute__((always_inline)) INLINE static float merger_spin_evolve(
- struct bpart* bpi, const struct bpart* bpj,
- const struct phys_const* constants) {
+__attribute__((always_inline)) INLINE static float
+black_hole_merger_spin_evolve(struct bpart* bpi, const struct bpart* bpj,
+ const struct phys_const* constants) {
/* Check if something is wrong with the masses. This is important and could
possibly happen as a result of jet spindown and mass loss at any time,
so we want to know about it. */
- if ((bpj->subgrid_mass <= 0.) || (bpi->subgrid_mass <= 0.)) {
+ if ((bpj->subgrid_mass <= 0.f) || (bpi->subgrid_mass <= 0.f)) {
error(
"Something went wrong with calculation of spin of a black hole "
" merger remnant. The black hole masses are %f and %f, instead of > "
@@ -1190,7 +1312,7 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve(
/* Check if the BHs have been spun down to 0. This is again an important
potential break point, we want to know about it. */
- if ((spin1 == 0.) || (spin2 == 0.)) {
+ if ((spin1 == 0.f) || (spin2 == 0.f)) {
error(
"Something went wrong with calculation of spin of a black hole "
" merger remnant. The black hole spins are %f and %f, instead of > 0.",
@@ -1207,8 +1329,8 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve(
/* We want to compute the direction of the orbital angular momentum of the
two BHs, which is used in the fits. Start by defining the coordinates in
- the frame of the centre of mass. For this we first need to compute the
- centre of mass coodinates and velocity. */
+ the frame of one of the BHs (it doesn't matter which one, the total
+ angular momentum is the same). */
const float centre_of_mass[3] = {
(m1 * bpi->x[0] + m2 * bpj->x[0]) / (m1 + m2),
(m1 * bpi->x[1] + m2 * bpj->x[1]) / (m1 + m2),
@@ -1250,7 +1372,7 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve(
m2 * (relative_coordinates_2[0] * relative_velocities_2[1] -
relative_coordinates_2[1] * relative_velocities_2[0])};
- /* The total, final orbital angular momentum. */
+ /* Calculate the orbital angular momentum itself. */
const float orbital_angular_momentum[3] = {
angular_momentum_1[0] + angular_momentum_2[0],
angular_momentum_1[1] + angular_momentum_2[1],
@@ -1283,6 +1405,7 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve(
const float j_BH_2 =
fabs(bpj->subgrid_mass * bpj->subgrid_mass * bpj->spin *
constants->const_newton_G / constants->const_speed_light_c);
+
float total_angular_momentum_direction[3] = {
j_BH_1 * spin_vec1[0] + j_BH_2 * spin_vec2[0] +
orbital_angular_momentum[0],
@@ -1327,11 +1450,28 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve(
/* Get the variable l used in the fit, see Eqn. 10 in Barausse & Rezolla
(2009). */
- const float l = l_variable(spin1, spin2, mass_ratio, sym_mass_ratio,
- cos_alpha, cos_beta, cos_gamma);
+ const float l = black_hole_l_variable(
+ spin1, spin2, mass_ratio, sym_mass_ratio, cos_alpha, cos_beta, cos_gamma);
+
+ const float l_vector[3] = {l * orbital_angular_momentum_direction[0],
+ l * orbital_angular_momentum_direction[1],
+ l * orbital_angular_momentum_direction[2]};
+
+ /* Final spin vector, constructed from the two spins and the auxilliary l
+ vector. */
+ const float spin_vector[3] = {
+ (spin_vec1[0] +
+ spin_vec2[0] * mass_ratio * mass_ratio * l_vector[0] * mass_ratio) /
+ ((1.f + mass_ratio) * (1.f + mass_ratio)),
+ (spin_vec1[1] +
+ spin_vec2[1] * mass_ratio * mass_ratio * l_vector[1] * mass_ratio) /
+ ((1.f + mass_ratio) * (1.f + mass_ratio)),
+ (spin_vec1[2] +
+ spin_vec2[2] * mass_ratio * mass_ratio * l_vector[2] * mass_ratio) /
+ ((1.f + mass_ratio) * (1.f + mass_ratio))};
#ifdef SWIFT_DEBUG_CHECKS
- if (l < 0.) {
+ if (l < 0.f) {
error(
"Something went wrong with calculation of spin of a black hole "
" merger remnant. The l factor is %f, instead of >= 0.",
@@ -1339,13 +1479,13 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve(
}
#endif
- /* Calculate the magnitude of final spin from Barausse & Rezolla (2009),
- Eqn. 6. */
+ /* Get magnitude of the final spin simply as the magnitude of the vector. */
const float final_spin_magnitude =
- final_spin(spin1, spin2, mass_ratio, cos_alpha, cos_beta, cos_gamma, l);
+ sqrtf(spin_vector[0] * spin_vector[0] + spin_vector[1] * spin_vector[1] +
+ spin_vector[2] * spin_vector[2]);
#ifdef SWIFT_DEBUG_CHECKS
- if (final_spin_magnitude <= 0.) {
+ if (final_spin_magnitude <= 0.f) {
error(
"Something went wrong with calculation of spin of a black hole "
" merger remnant. The final spin magnitude is %f, instead of > 0.",
@@ -1355,15 +1495,15 @@ __attribute__((always_inline)) INLINE static float merger_spin_evolve(
/* Assign the final spin value to the BH, but also make sure we don't go
above 0.998 nor below 0.001. */
- bpi->spin = min(final_spin_magnitude, 0.998);
- if (fabsf(bpi->spin) < 0.01) {
- bpi->spin = 0.01;
+ bpi->spin = min(final_spin_magnitude, 0.998f);
+ if (fabsf(bpi->spin) < 0.01f) {
+ bpi->spin = 0.01f;
}
/* Assign the directions of the spin to the BH. */
- bpi->angular_momentum_direction[0] = total_angular_momentum_direction[0];
- bpi->angular_momentum_direction[1] = total_angular_momentum_direction[1];
- bpi->angular_momentum_direction[2] = total_angular_momentum_direction[2];
+ bpi->angular_momentum_direction[0] = spin_vector[0] / final_spin_magnitude;
+ bpi->angular_momentum_direction[1] = spin_vector[1] / final_spin_magnitude;
+ bpi->angular_momentum_direction[2] = spin_vector[2] / final_spin_magnitude;
/* Finally we also want to calculate the fraction of total mass-energy
lost during the merger to gravitational waves. We use Eqn. 16 and 18
diff --git a/src/tracers/EAGLE/tracers.h b/src/tracers/EAGLE/tracers.h
index aa74338bd9aaee5264daa0457deaba4b2cfb1f17..e7fbf37636a35b5e3080382cc0e4cbbf0a4fa9d6 100644
--- a/src/tracers/EAGLE/tracers.h
+++ b/src/tracers/EAGLE/tracers.h
@@ -220,6 +220,8 @@ static INLINE void tracers_first_init_xpart(
xp->tracers_data.last_AGN_jet_feedback_scale_factor = 0.f;
xp->tracers_data.last_AGN_jet_feedback_time = 0.f;
xp->tracers_data.last_jet_kick_velocity = 0.f;
+ xp->tracers_data.last_jet_kick_accretion_mode = BH_thick_disc;
+ xp->tracers_data.last_jet_kick_BH_id = 0;
}
/**
@@ -324,7 +326,7 @@ static INLINE void tracers_after_black_holes_feedback(
static INLINE void tracers_after_jet_feedback(
const struct part *p, struct xpart *xp, const int with_cosmology,
const float scale_factor, const double time, const double delta_energy,
- const float vel_kick) {
+ const float vel_kick, BH_accretion_modes accretion_mode, long long id) {
if (with_cosmology)
xp->tracers_data.last_AGN_jet_feedback_scale_factor = scale_factor;
@@ -333,6 +335,8 @@ static INLINE void tracers_after_jet_feedback(
xp->tracers_data.hit_by_jet_feedback++;
xp->tracers_data.jet_feedback_energy += delta_energy;
xp->tracers_data.last_jet_kick_velocity = vel_kick;
+ xp->tracers_data.last_jet_kick_accretion_mode = accretion_mode;
+ xp->tracers_data.last_jet_kick_BH_id = id;
}
/**
diff --git a/src/tracers/EAGLE/tracers_io.h b/src/tracers/EAGLE/tracers_io.h
index 4cb34b461ea9d3de317e0fe5217929940b8693f8..53c02b192f619778228bb8cc6c32820e725368d8 100644
--- a/src/tracers/EAGLE/tracers_io.h
+++ b/src/tracers/EAGLE/tracers_io.h
@@ -228,7 +228,21 @@ __attribute__((always_inline)) INLINE static int tracers_write_particles(
tracers_data.last_jet_kick_velocity,
"Kick velocity at last AGN jet event.");
- return 15;
+ list[15] = io_make_output_field("LastAGNJetKickMode", CHAR, 1,
+ UNIT_CONV_NO_UNITS, 0.f, xparts,
+ tracers_data.last_jet_kick_accretion_mode,
+ "The accretion/feedback mode the BH was "
+ "in when it kicked this particle. 0 "
+ "corresponds to the thick disc, 1 to the "
+ "thin disc and 2 to the slim disc.");
+
+ list[16] =
+ io_make_output_field("LastAGNJetKickBHId", CHAR, 1, UNIT_CONV_NO_UNITS,
+ 0.f, xparts, tracers_data.last_jet_kick_BH_id,
+ "The id of the BH that last kicked this "
+ "particle.");
+
+ return 17;
} else {
return 11;
}
@@ -376,7 +390,21 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles(
tracers_data.last_jet_kick_velocity,
"Kick velocity at last AGN jet event.");
- return 15;
+ list[15] = io_make_output_field("LastAGNJetKickMode", CHAR, 1,
+ UNIT_CONV_NO_UNITS, 0.f, sparts,
+ tracers_data.last_jet_kick_accretion_mode,
+ "The accretion/feedback mode the BH was "
+ "in when it kicked this particle. 0 "
+ "corresponds to the thick disc, 1 to the "
+ "thin disc and 2 to the slim disc.");
+
+ list[16] =
+ io_make_output_field("LastAGNJetKickBHId", CHAR, 1, UNIT_CONV_NO_UNITS,
+ 0.f, sparts, tracers_data.last_jet_kick_BH_id,
+ "The id of the BH that last kicked this "
+ "particle.");
+
+ return 17;
} else {
return 11;
}
diff --git a/src/tracers/EAGLE/tracers_struct.h b/src/tracers/EAGLE/tracers_struct.h
index 5e8c7df73249f4ac0149c7dcfc989a54a894acc9..735025d7cb0c59866420bb3a85b27b3f9285be2a 100644
--- a/src/tracers/EAGLE/tracers_struct.h
+++ b/src/tracers/EAGLE/tracers_struct.h
@@ -95,8 +95,15 @@ struct tracers_xpart_data {
/*! Has this particle been hit by AGN feedback? */
char hit_by_AGN_feedback;
- /* Kick velocity at last AGN jet event */
+ /*! Kick velocity at last AGN jet event */
float last_jet_kick_velocity;
+
+ /*! The accretion/feedback mode of the BH when this particle was last
+ * kicked */
+ enum BH_accretion_modes last_jet_kick_accretion_mode;
+
+ /*! The ID of the BH that did the last kick */
+ long long last_jet_kick_BH_id;
};
/**