Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • dc-oman1/swiftsim
  • swift/swiftsim
  • pdraper/swiftsim
  • tkchan/swiftsim
  • dc-turn5/swiftsim
5 results
Select Git revision
Show changes
Showing
with 1161 additions and 560 deletions
...@@ -46,38 +46,40 @@ Below we give an example of parameter choices applicable for e.g. a 50 Mpc box. ...@@ -46,38 +46,40 @@ Below we give an example of parameter choices applicable for e.g. a 50 Mpc box.
AGN_use_deterministic_feedback: 1 # Deterministic (1) or stochastic (0) AGN feedback model AGN_use_deterministic_feedback: 1 # Deterministic (1) or stochastic (0) AGN feedback model
AGN_feedback_model: Isotropic # AGN feedback model (Isotropic or MinimumDistance) AGN_feedback_model: Isotropic # AGN feedback model (Isotropic or MinimumDistance)
minimum_timestep_yr: 1000.0 # Minimum time-step of black-hole particles minimum_timestep_yr: 1000.0 # Minimum time-step of black-hole particles
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate. 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]. 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. 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.1 # 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. 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.
seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1. 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.
AGN_jet_velocity_model: BlackHoleMass # 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. seed_spin: 0.01 # The (randomly-directed) black hole spin assigned to BHs when they are seeded. Should be strictly between 0 and 1.
v_jet_km_p_s: 10000. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s. 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_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_km_p_s: 3160. # Jet velocity to use if 'AGN_jet_velocity_model' is 'Constant'. Units are km/s.
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_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.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_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_mass_loading: 400. # The constant mass loading to use if 'AGN_jet_velocity_model' is MassLoading. 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_xi: 0.707 # The numerical multiplier by which the jet velocity formula is scaled, if 'AGN_jet_velocity_model' is 'Local' or 'SoundSpeed'. The appropriate values (to exactly obtain the formulas as derived) are 0.63 and 0.707 for the two, respectively. v_jet_max_km_p_s: 1e4 # The maximal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading' or 'Local'.
v_jet_min_km_p_s: 500 # The minimal jet velocity. This is used if 'AGN_jet_velocity_model' is 'BlackHoleMass', 'MassLoading', 'Local' or 'SoundSpeed'.
opening_angle_in_degrees: 7.5 # The half-opening angle of the jet in degrees. Should use values < 15 unless for tests. 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). 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. 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. 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: Variable # How the accretion efficiencies are calculated for the thick accretion disc. If 'Constant', the value of 'accretion_efficiency_thick' will be used. If 'Variable', 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.
ADIOS_s: 0.5 # The exponent of the scaling between accretion efficiency and transition radius of the accretion disc, used if 'accretion_efficiency_mode' is 'Variable'.
ADIOS_R_in: 1e4 # The normalisation (the value) of the transition radius of the accretion disc at the critical Eddington ratio (0.01), used if 'accretion_efficiency_mode' is 'Variable'.
fix_radiative_efficiency: 0 # Global switch whether to fix the radiative efficiency to a particular value [1], or use a spin-dependant formula [0]. 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. 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_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.
ADIOS_R_in: 30. # If include_ADIOS_accr_suppression is set to 1, this parameter controls the inner radius within which winds are not important.
ADIOS_s: 0.4 # Slope of the accretion rate - radius relationship if include_ADIOS_accr_suppression is set to 1.
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. 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.
TD_SD_eps_r_threshold: 0.75 # Parameter controlling the transition from thin to slim disk. Accretion disk will be slim if radiative efficiency satisfies eps_slim < TD_SD_eps_r_threshold * eps_thin. This parameter is only used if include_slim_disk is set to 1. 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. 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, separated by a value of ``mdot_crit_ADAF``), 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.
If you set ``use_var_v_jet: 0``, you will need to change ``v_jet``, the kicking velocity of particles, depending on what system you are simulating. You should typically choose values at least 10 times larger than the sound speed of the hot gas in your most massive haloes (e.g. 1500 km/s for a MW-type galaxy and 10 000 km/s for a :math:`10^{14}` :math:`\mathrm{M}_\odot` halo). If, on the other hand, you set ``use_var_v_jet: 1``, the launching velocities will vary on their own depending on the typical sound speed (virial velocity) of the hot gas in the haloes. You then need to set ``v_jet_cs_ratio`` to values :math:`\gg1` (10-30 works well) in order to have significant shocking. Turning on ``use_jets_in_thin_disc`` or ``use_ADIOS_winds`` will cause jets to also be used in the thin disk and winds (thermal isotropic feedback) in the thick disk. Similarly, ``use_ADIOS_winds`` will lead to winds in the slim disk, with the value of the parameter being used to rescale the formula that is implemented (a value of 1 leading to maximal winds).
...@@ -188,6 +188,39 @@ def L_adv(x, alpha): ...@@ -188,6 +188,39 @@ def L_adv(x, alpha):
) / eta(x, alpha) ) / eta(x, alpha)
def jet_eff(f_Edd, a):
horizon_ang_vel = abs(a) / (2.0 * (1.0 + np.sqrt(1 - a * a)))
phi = -20.2 * a ** 3 - 14.9 * a ** 2 + 34.0 * a + 52.6
phi = phi * (f_Edd / 1.88) ** 1.29 / (1 + (f_Edd / 1.88) ** 1.29)
return (
0.05
/ (4.0 * np.pi)
* phi ** 2
* horizon_ang_vel ** 2
* (1.0 + 1.38 * horizon_ang_vel ** 2 - 9.2 * horizon_ang_vel ** 4)
)
def s_HD(f_Edd, a):
xi = f_Edd * 0.017
s_min = 0.86 - 1.94 * a
L_ISCO = 0.385 * (1.0 + 2.0 * np.sqrt(3.0 * r_isco(a) - 2.0))
s_thin = L_ISCO - 2.0 * a * (1.0 - eps_NT(a))
return (s_thin + s_min * xi) / (1 + xi)
def s(f_Edd, a):
horizon_ang_vel = abs(a) / (2.0 * (1.0 + np.sqrt(1 - a * a)))
k_EM = 0.23 * np.ones(np.size(a))
k_EM[a > 0] = np.minimum(0.1 + 0.5 * a[a > 0], 0.35 * np.ones(np.size(a[a > 0])))
s_EM = (
-1 * a / abs(a) * jet_eff(f_Edd, a) * (1.0 / (k_EM * horizon_ang_vel) - 2.0 * a)
)
return s_HD(f_Edd, a) + s_EM
a = np.arange(-1, 1, 0.0001) a = np.arange(-1, 1, 0.0001)
mdotcrit1 = m_dot_crit1(a, 0.5) mdotcrit1 = m_dot_crit1(a, 0.5)
mdotcrit2 = m_dot_crit2(a, 0.5) mdotcrit2 = m_dot_crit2(a, 0.5)
...@@ -196,6 +229,7 @@ m_a_75 = [find_root(x, 0.75) for x in a] ...@@ -196,6 +229,7 @@ m_a_75 = [find_root(x, 0.75) for x in a]
m_a_50 = [find_root(x, 0.5) for x in a] m_a_50 = [find_root(x, 0.5) for x in a]
import matplotlib import matplotlib
import pylab
matplotlib.use("Agg") matplotlib.use("Agg")
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -204,17 +238,17 @@ import matplotlib.gridspec as gridspec ...@@ -204,17 +238,17 @@ import matplotlib.gridspec as gridspec
fig = plt.figure(figsize=(8, 6)) fig = plt.figure(figsize=(8, 6))
plt.style.use("classic") plt.style.use("classic")
plt.fill_between(a, [0.0001 for x in a], [0.028 for x in a], color="blue", alpha=0.2) plt.fill_between(a, [0.0001 for x in a], [0.01 for x in a], color="blue", alpha=0.2)
plt.fill_between(a, [0.028 for x in a], mdotcrit1, color="red", alpha=0.2) plt.fill_between(a, [0.01 for x in a], [1 for x in a], color="red", alpha=0.2)
plt.fill_between(a, mdotcrit1, [375 for x in a], color="orange", alpha=0.2) plt.fill_between(a, [1 for x in a], [375 for x in a], color="orange", alpha=0.2)
plt.ylabel("$\dot{m}$", fontsize=24, usetex=True) plt.ylabel(r"$f_\mathrm{Edd}$", fontsize=24, usetex=True)
plt.xlabel("$a$", fontsize=24, usetex=True) plt.xlabel(r"$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in") plt.tick_params(axis="y", right=True, direction="in")
plt.yscale("log") plt.yscale("log")
plt.axis([-1, 1, 0.001, 100]) plt.axis([-1, 1, 0.0001, 100])
plt.text(-0.22, 0.004, "Thick disc", fontsize=20) plt.text(-0.22, 0.0008, "Thick disk", fontsize=20)
plt.text(-0.2, 0.33, "Thin disc", fontsize=20) plt.text(-0.2, 0.08, "Thin disk", fontsize=20)
plt.text(-0.2, 18, "Slim disc", fontsize=20) plt.text(-0.2, 8, "Slim disk", fontsize=20)
plt.minorticks_on() plt.minorticks_on()
plt.tick_params( plt.tick_params(
axis="x", axis="x",
...@@ -260,6 +294,7 @@ plt.tick_params( ...@@ -260,6 +294,7 @@ plt.tick_params(
plt.savefig("modes.png", bbox_inches="tight") plt.savefig("modes.png", bbox_inches="tight")
plt.close() plt.close()
a = np.arange(-1, 1, 0.0001)
phi = -20.2 * a ** 3 - 14.9 * a ** 2 + 34.0 * a + 52.6 phi = -20.2 * a ** 3 - 14.9 * a ** 2 + 34.0 * a + 52.6
horizon_ang_vel = a / (2 * (1 + np.sqrt(1 - a ** 2))) horizon_ang_vel = a / (2 * (1 + np.sqrt(1 - a ** 2)))
jet_factor = ( jet_factor = (
...@@ -271,14 +306,14 @@ jet_factor = ( ...@@ -271,14 +306,14 @@ jet_factor = (
* horizon_ang_vel ** 2 * horizon_ang_vel ** 2
* (1.0 + 1.38 * horizon_ang_vel ** 2 - 9.2 * horizon_ang_vel ** 4) * (1.0 + 1.38 * horizon_ang_vel ** 2 - 9.2 * horizon_ang_vel ** 4)
) )
Z1_j = np.array( Z_1 = np.array(
[ [
1 + (1 - x ** 2) ** 0.333 * ((1 + abs(x)) ** 0.333 + (1 - abs(x)) ** 0.333) 1 + (1 - x ** 2) ** 0.333 * ((1 + abs(x)) ** 0.333 + (1 - abs(x)) ** 0.333)
for x in a for x in a
] ]
) )
Z2_j = np.array(np.sqrt(3 * a ** 2 + Z1_j ** 2)) Z_2 = np.array(np.sqrt(3 * a ** 2 + Z_1 ** 2))
r_iso = 3 + Z2_j - np.sign(np.array(a)) * np.sqrt((3 - Z1_j) * (3 + Z1_j + 2 * Z2_j)) r_iso = 3 + Z_2 - np.sign(np.array(a)) * np.sqrt((3 - Z_1) * (3 + Z_1 + 2 * Z_2))
eps_TD = 1 - np.sqrt(1 - 2 / (3 * r_iso)) eps_TD = 1 - np.sqrt(1 - 2 / (3 * r_iso))
eps_ADAF1 = 0.144 * (6 / r_iso) * eps_TD * min(1, 0.028 / 0.0044) eps_ADAF1 = 0.144 * (6 / r_iso) * eps_TD * min(1, 0.028 / 0.0044)
eps_ADAF2 = 0.144 * (6 / r_iso) * eps_TD * min(1, 0.001 / 0.0044) eps_ADAF2 = 0.144 * (6 / r_iso) * eps_TD * min(1, 0.001 / 0.0044)
...@@ -286,6 +321,7 @@ Jet_ADAF = jet_factor * 0.3 ...@@ -286,6 +321,7 @@ Jet_ADAF = jet_factor * 0.3
Jet_SD = 0.22 * jet_factor Jet_SD = 0.22 * jet_factor
Jet_TD1 = 10 ** -3 * 0.1 ** (-0.1) * 100 ** 0.2 * 10 ** (2 * 0.1) * jet_factor Jet_TD1 = 10 ** -3 * 0.1 ** (-0.1) * 100 ** 0.2 * 10 ** (2 * 0.1) * jet_factor
Jet_TD2 = 10 ** -3 * 0.1 ** (-0.1) * 10 ** (-1 * 0.1) * jet_factor Jet_TD2 = 10 ** -3 * 0.1 ** (-0.1) * 10 ** (-1 * 0.1) * jet_factor
eps_SD1 = ( eps_SD1 = (
1 1
/ 10 / 10
...@@ -317,109 +353,82 @@ mdot_bh_TD1 = (1 - Jet_TD1 / 4.447) * (1 - eps_TD - Jet_TD1) ...@@ -317,109 +353,82 @@ mdot_bh_TD1 = (1 - Jet_TD1 / 4.447) * (1 - eps_TD - Jet_TD1)
mdot_bh_TD2 = (1 - Jet_TD2 / 4.447) * (1 - eps_TD - Jet_TD2) mdot_bh_TD2 = (1 - Jet_TD2 / 4.447) * (1 - eps_TD - Jet_TD2)
fig = plt.figure(figsize=(18, 4)) def omega(spin):
fig.subplots_adjust(wspace=0, hspace=0, top=1, bottom=0) return spin / (2 * (1 + np.sqrt(1 - spin ** 2)))
gs = gridspec.GridSpec(1, 3, width_ratios=[1, 1, 1])
fig = plt.figure(figsize=(13, 4))
fig.subplots_adjust(top=1, bottom=0, wspace=0.25)
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 1])
plt.style.use("classic") plt.style.use("classic")
plt.subplot(gs[0]) plt.subplot(gs[0])
plt.plot( plt.plot(
a, a,
eps_ADAF2, 100 * 0.005 * (1 + 3 * (phi / 50) ** 2 * (horizon_ang_vel / 0.2) ** 2),
linewidth=2, linewidth=2,
label="$\epsilon_\mathrm{rad}(\dot{m}<0.0044)$", label=r"$\epsilon_\mathrm{wind,thick}$",
color="red", color="blue",
) )
plt.plot( plt.plot(
a, a,
eps_ADAF1, 100 * 0.1 * (1 - np.sqrt(1 - 2 / (3 * r_iso))),
linewidth=2, linewidth=2,
label="$\epsilon_\mathrm{rad}(\dot{m}=0.028)$", label=r"$\epsilon_\mathrm{f}\epsilon_\mathrm{rad,NT}$ $\mathrm{(thin}$ $\mathrm{disc})$",
color="red", color="red",
linestyle="--",
)
plt.plot(a, 0.97 * Jet_ADAF, linewidth=2, label="$\epsilon_\mathrm{jet}$", color="blue")
plt.fill_between(a, eps_ADAF1, eps_ADAF2, color="red", alpha=0.2)
plt.ylabel("$\epsilon_\mathrm{feedback}$", fontsize=24, usetex=True)
plt.xlabel("$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in")
plt.legend(loc="upper left", prop={"size": 13})
plt.xticks([-1, -0.5, 0, 0.5, 1], [-1, -0.5, 0, 0.5, 1])
plt.yticks(
[0.0001, 0.001, 0.01, 0.1, 1, 10, 100],
["", 10 ** (-3), 10 ** (-2), 10 ** (-1), 10 ** (-0), 10 ** (1), 10 ** (2)],
)
plt.minorticks_on()
plt.axis([-1, 1, 0.0001, 10])
plt.yscale("log")
plt.tick_params(
axis="x",
direction="in",
bottom=True,
top=True,
length=8,
width=1.2,
which="major",
labelsize=16,
) )
plt.tick_params( plt.plot(
axis="y", a,
direction="in", 100
left=True, * 0.0635
right=True, * (1 + ((1 / 1.88) ** 1.29 / (1 + (1 / 1.88) ** 1.29) * phi / 50) ** 2)
length=8, * np.maximum((1 - 8 * omega(a) ** 2 + 1 * omega(a)), np.zeros(np.size(a))),
width=1.2, linestyle=":",
which="major", linewidth=1.5,
labelsize=16, label=r"$\epsilon_\mathrm{wind,slim},$ $f_\mathrm{Edd}=1$",
) color="orange",
plt.tick_params(
axis="x",
direction="in",
bottom=True,
top=True,
length=4,
width=0.9,
which="minor",
labelsize=16,
)
plt.tick_params(
axis="y",
direction="in",
left=True,
right=True,
length=4,
width=0.9,
which="minor",
labelsize=16,
) )
plt.title("Thick disc", fontsize=16)
plt.subplot(gs[1])
plt.plot(a, 0.97 * eps_TD, linewidth=2, label="$\epsilon_\mathrm{rad}$", color="red")
plt.plot( plt.plot(
a, a,
Jet_TD2, 100
linewidth=2, * 0.0635
label="$\epsilon_\mathrm{jet}(\dot{m}=0.028,M_\mathrm{BH}=10^9 \mathrm{M}_\odot)$", * (1 + ((10 / 1.88) ** 1.29 / (1 + (10 / 1.88) ** 1.29) * phi / 50) ** 2)
color="blue", * np.maximum((1 - 8 * omega(a) ** 2 + 1 * omega(a)), np.zeros(np.size(a))),
linestyle="-.",
linewidth=1.5,
label=r"$\epsilon_\mathrm{wind,slim},$ $f_\mathrm{Edd}=10$",
color="orange",
) )
plt.plot( plt.plot(
a, a,
Jet_TD1, 100
linewidth=2, * 0.0635
label="$\epsilon_\mathrm{jet}(\dot{m}=1,M_\mathrm{BH}=10^6 \mathrm{M}_\odot)$", * (1 + ((100 / 1.88) ** 1.29 / (1 + (100 / 1.88) ** 1.29) * phi / 50) ** 2)
color="blue", * np.maximum((1 - 8 * omega(a) ** 2 + 1 * omega(a)), np.zeros(np.size(a))),
linestyle="--", linestyle="--",
linewidth=1.5,
label=r"$\epsilon_\mathrm{wind,slim},$ $f_\mathrm{Edd}=100$",
color="orange",
) )
plt.fill_between(a, Jet_TD1, Jet_TD2, color="blue", alpha=0.2) plt.plot(
plt.xlabel("$a$", fontsize=24, usetex=True) a,
100
* 0.0635
* (1 + ((1000 / 1.88) ** 1.29 / (1 + (1000 / 1.88) ** 1.29) * phi / 50) ** 2)
* np.maximum((1 - 8 * omega(a) ** 2 + 1 * omega(a)), np.zeros(np.size(a))),
linestyle="-",
linewidth=1.5,
label=r"$\epsilon_\mathrm{wind,slim},$ $f_\mathrm{Edd}=1000$",
color="orange",
)
plt.fill_between(a, eps_ADAF1, eps_ADAF2, color="red", alpha=0.2)
plt.ylabel(r"$\epsilon_\mathrm{wind}$ $[\%]$", fontsize=24, usetex=True)
plt.xlabel(r"$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in") plt.tick_params(axis="y", right=True, direction="in")
plt.yscale("log") pylab.legend(loc="upper left", prop={"size": 12}, ncol=2)
plt.legend(loc="upper left", prop={"size": 13})
plt.xticks([-1, -0.5, 0, 0.5, 1], ["", -0.5, 0, 0.5, 1])
plt.axis([-1, 1, 0.0001, 10])
plt.yticks([0.001, 0.01, 0.1, 1, 10], ["", "", "", "", ""])
plt.minorticks_on() plt.minorticks_on()
plt.axis([-1, 1, 0, 25])
plt.tick_params( plt.tick_params(
axis="x", axis="x",
direction="in", direction="in",
...@@ -460,29 +469,58 @@ plt.tick_params( ...@@ -460,29 +469,58 @@ plt.tick_params(
which="minor", which="minor",
labelsize=16, labelsize=16,
) )
plt.title("Thin disc", fontsize=16) plt.title("Wind efficiency", fontsize=16)
plt.subplot(gs[2]) plt.subplot(gs[1])
plt.plot( plt.plot(
a, eps_SD1, linewidth=2, label="$\epsilon_\mathrm{rad}(\dot{m}=1)$", color="red" a, 100 * Jet_ADAF, linewidth=2, label=r"$\epsilon_\mathrm{jet,thick}$", color="blue"
) )
plt.plot( plt.plot(
a, a,
eps_SD2, 100 * Jet_ADAF * ((0.01 / 1.88) ** 1.29 / (1 + (0.01 / 1.88) ** 1.29)) ** 2,
linewidth=2, linewidth=1.5,
label="$\epsilon_\mathrm{rad}(\dot{m}=50)$", linestyle=":",
label=r"$\epsilon_\mathrm{jet,thin},$ $f_\mathrm{Edd}=0.01$",
color="red", color="red",
)
plt.plot(
a,
100 * Jet_ADAF * ((0.1 / 1.88) ** 1.29 / (1 + (0.1 / 1.88) ** 1.29)) ** 2,
linewidth=1.5,
linestyle="-.",
label=r"$\epsilon_\mathrm{jet,thin},$ $f_\mathrm{Edd}=0.1$",
color="red",
)
plt.plot(
a,
100 * Jet_ADAF * ((1 / 1.88) ** 1.29 / (1 + (1 / 1.88) ** 1.29)) ** 2,
linewidth=1.5,
linestyle="--", linestyle="--",
label=r"$\epsilon_\mathrm{jet,thin},$ $f_\mathrm{Edd}=1$",
color="red",
) )
plt.plot(a, Jet_SD, linewidth=2, label="$\epsilon_\mathrm{jet}$", color="blue") plt.plot(
plt.fill_between(a, eps_SD1, eps_SD2, color="red", alpha=0.2) a,
plt.xlabel("$a$", fontsize=24, usetex=True) 100 * Jet_ADAF * ((10 / 1.88) ** 1.29 / (1 + (10 / 1.88) ** 1.29)) ** 2,
linewidth=1.5,
linestyle="--",
label=r"$\epsilon_\mathrm{jet,slim},$ $f_\mathrm{Edd}=10$",
color="orange",
)
plt.plot(
a,
100 * Jet_ADAF * ((100 / 1.88) ** 1.29 / (1 + (100 / 1.88) ** 1.29)) ** 2 - 2,
linewidth=1.5,
linestyle="-",
label=r"$\epsilon_\mathrm{jet,slim},$ $f_\mathrm{Edd}=100$",
color="orange",
)
plt.ylabel(r"$\epsilon_\mathrm{jet}$ $[\%]$", fontsize=24, usetex=True)
plt.xlabel(r"$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in") plt.tick_params(axis="y", right=True, direction="in")
plt.yscale("log") pylab.legend(loc="upper left", prop={"size": 15})
plt.legend(loc="upper left", prop={"size": 13}) plt.axis([-1, 1, 0, 200])
plt.xticks([-1, -0.5, 0, 0.5, 1], ["", -0.5, 0, 0.5, 1])
plt.axis([-1, 1, 0.0001, 10])
plt.yticks([0.001, 0.01, 0.1, 1, 10], ["", "", "", "", ""])
plt.minorticks_on() plt.minorticks_on()
plt.tick_params( plt.tick_params(
axis="x", axis="x",
...@@ -524,7 +562,7 @@ plt.tick_params( ...@@ -524,7 +562,7 @@ plt.tick_params(
which="minor", which="minor",
labelsize=16, labelsize=16,
) )
plt.title("Slim disc", fontsize=16) plt.title("Jet efficiency", fontsize=16)
plt.savefig("efficiencies.png", bbox_inches="tight") plt.savefig("efficiencies.png", bbox_inches="tight")
...@@ -532,7 +570,7 @@ L_isco1 = [2 / 3 * 1 / np.sqrt(3) * (1 + 2 * np.sqrt(3 * r_isco(x) - 2)) for x i ...@@ -532,7 +570,7 @@ L_isco1 = [2 / 3 * 1 / np.sqrt(3) * (1 + 2 * np.sqrt(3 * r_isco(x) - 2)) for x i
plt.style.use("classic") plt.style.use("classic")
fig = plt.figure(figsize=(8, 6), linewidth=4) fig = plt.figure(figsize=(8, 6), linewidth=4)
plt.plot(a, L_isco1, linewidth=2, label="$\ell_\mathrm{ISCO}$", color="red") plt.plot(a, L_isco1, linewidth=2, label=r"$\ell_\mathrm{ISCO}$", color="red")
plt.plot( plt.plot(
a, a,
0.45 * np.array(L_isco1), 0.45 * np.array(L_isco1),
...@@ -562,8 +600,8 @@ plt.plot( ...@@ -562,8 +600,8 @@ plt.plot(
label=r"$\ell_\mathrm{adv},\alpha=0.3$", label=r"$\ell_\mathrm{adv},\alpha=0.3$",
color="purple", color="purple",
) )
plt.ylabel("$\ell_\mathrm{in}$", fontsize=24, usetex=True) plt.ylabel(r"$\ell_\mathrm{in}$", fontsize=24, usetex=True)
plt.xlabel("$a$", fontsize=24, usetex=True) plt.xlabel(r"$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in") plt.tick_params(axis="y", right=True, direction="in")
plt.legend(loc="upper right", prop={"size": 14}) plt.legend(loc="upper right", prop={"size": 14})
plt.minorticks_on() plt.minorticks_on()
...@@ -670,125 +708,75 @@ da_SD_Benson = ( ...@@ -670,125 +708,75 @@ da_SD_Benson = (
) )
fig = plt.figure(figsize=(18, 4)) fig = plt.figure(figsize=(7, 5))
fig.subplots_adjust(wspace=0, hspace=0, top=1, bottom=0)
gs = gridspec.GridSpec(1, 3, width_ratios=[1, 1, 1])
plt.style.use("classic") plt.style.use("classic")
plt.subplot(gs[0]) z1 = np.array(
plt.plot(a, da_TD_acc_only, linewidth=2, label="Accretion only", color="black") [
plt.plot(a, da_TD_Benson, linewidth=1.5, label="Jet spindown included", color="blue") 1 + (1 - x ** 2) ** 0.333 * ((1 + abs(x)) ** 0.333 + (1 - abs(x)) ** 0.333)
plt.plot(a, [0 for x in a], linewidth=1.5, color="black", linestyle="--") for x in a
plt.ylabel("$\mathrm{d}a/\mathrm{d}\ln M_\mathrm{BH,0}$", fontsize=24, usetex=True) ]
plt.xlabel("$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in")
plt.legend(loc="lower left", prop={"size": 15})
plt.minorticks_on()
plt.axis([-1, 1, -4, 7])
plt.tick_params(
axis="x",
direction="in",
bottom=True,
top=True,
length=8,
width=1.2,
which="major",
labelsize=16,
) )
plt.tick_params( z2 = np.array(np.sqrt(3 * a ** 2 + z1 ** 2))
axis="y", r_iso = 3 + z2 - np.sign(np.array(a)) * np.sqrt((3 - z1) * (3 + z1 + 2 * z2))
direction="in", da_TD_acc_only = 2 / 3 * 1 / np.sqrt(3) * (
left=True, 1 + 2 * np.sqrt(3 * r_iso - 2)
right=True, ) - 2 * a * np.sqrt(1 - 2 / (3 * r_iso))
length=8, da_ADAF_Narayan = (
width=1.2, 0.45 - 12.53 * a - 7.8 * a ** 2 + 9.44 * a ** 3 + 5.71 * a ** 4 - 4.03 * a ** 5
which="major",
labelsize=16,
) )
plt.tick_params(
axis="x", plt.plot(a, da_ADAF_Narayan, linewidth=2, label="Thick disk", color="blue")
direction="in", plt.plot(
bottom=True, a,
top=True, s(0.01, a),
length=4, linewidth=2,
width=0.9, label=r"Thin disk, $f_\mathrm{Edd}=0.01$",
which="minor", linestyle=":",
labelsize=16, color="red",
) )
plt.tick_params( plt.plot(
axis="y", a,
direction="in", s(0.1, a),
left=True, linewidth=2,
right=True, label=r"Thin disk, $f_\mathrm{Edd}=0.1$",
length=4, linestyle="-.",
width=0.9, color="red",
which="minor",
labelsize=16,
) )
plt.title("Thin disc", fontsize=16) plt.plot(
a,
plt.subplot(gs[1]) s(1, a),
plt.plot(a, da_ADAF_acc_only, linewidth=2, label="Accretion only", color="black") linewidth=2,
plt.plot(a, da_ADAF_Benson, linewidth=1.5, label="Jet spindown included", color="blue") label=r"Thin disk, $f_\mathrm{Edd}=1$",
plt.plot(a, [0 for x in a], linewidth=1.5, color="black", linestyle="--") linestyle="--",
plt.xlabel("$a$", fontsize=24, usetex=True) color="red",
plt.tick_params(axis="y", right=True, direction="in")
plt.xticks([-1.0, -0.5, 0.0, 0.5, 1.0], ["", -0.5, 0.0, 0.5, 1.0])
plt.yticks([-8, -6, -4, -2, 0, 2, 4, 6, 8], ["", "", "", "", "", "", "", "", ""])
plt.minorticks_on()
plt.axis([-1, 1, -4, 7])
plt.tick_params(
axis="x",
direction="in",
bottom=True,
top=True,
length=8,
width=1.2,
which="major",
labelsize=16,
) )
plt.tick_params( plt.plot(
axis="y", a,
direction="in", s(10, a),
left=True, linewidth=2,
right=True, label=r"Slim disk, $f_\mathrm{Edd}=10$",
length=8, linestyle="--",
width=1.2, color="orange",
which="major",
labelsize=16,
) )
plt.tick_params( plt.plot(
axis="x", a,
direction="in", s(100, a),
bottom=True, linewidth=2,
top=True, label=r"Slim disk, $f_\mathrm{Edd}=100$",
length=4, linestyle="-",
width=0.9, color="orange",
which="minor",
labelsize=16,
) )
plt.tick_params( plt.plot(a, [0 for x in a], linewidth=1.0, color="black", linestyle="--")
axis="y", plt.plot([-0.0001, 0.0001], [-200, 200], linewidth=1.0, color="black", linestyle="--")
direction="in", plt.ylabel(
left=True, r"$\mathrm{d}a/(\mathrm{d} M_\mathrm{BH,0}/M_\mathrm{BH})$", fontsize=24, usetex=True
right=True,
length=4,
width=0.9,
which="minor",
labelsize=16,
) )
plt.title("Thick disc", fontsize=16)
plt.subplot(gs[2])
plt.plot(a, da_SD_acc_only, linewidth=2, label="Accretion only", color="black")
plt.plot(a, da_SD_Benson, linewidth=1.5, label="Jet spindown included", color="blue")
plt.plot(a, [0 for x in a], linewidth=1.5, color="black", linestyle="--")
plt.xlabel("$a$", fontsize=24, usetex=True) plt.xlabel("$a$", fontsize=24, usetex=True)
plt.tick_params(axis="y", right=True, direction="in") plt.tick_params(axis="y", right=True, direction="in")
plt.xticks([-1.0, -0.5, 0.0, 0.5, 1.0], ["", -0.5, 0.0, 0.5, 1.0]) pylab.legend(loc="lower left", prop={"size": 13})
plt.yticks([-8, -6, -4, -2, 0, 2, 4, 6, 8], ["", "", "", "", "", "", "", "", ""])
plt.minorticks_on() plt.minorticks_on()
plt.axis([-1, 1, -4, 7]) plt.axis([-1, 1, -10, 10])
plt.tick_params( plt.tick_params(
axis="x", axis="x",
direction="in", direction="in",
...@@ -829,6 +817,5 @@ plt.tick_params( ...@@ -829,6 +817,5 @@ plt.tick_params(
which="minor", which="minor",
labelsize=16, labelsize=16,
) )
plt.title("Slim disc", fontsize=16)
plt.savefig("spinup.png", bbox_inches="tight") plt.savefig("spinup.png", bbox_inches="tight")
...@@ -11,7 +11,7 @@ Here we provide a comprehensive summary of the model. In order to more easily vi ...@@ -11,7 +11,7 @@ Here we provide a comprehensive summary of the model. In order to more easily vi
Any model for realistic AGN jets must include black hole spin since jet powers depend steeply on spin, and because including spin provides a well-defined direction for the jets to be launched in. The spin (angular momentum) of BHs is best represented through the dimensionlesss spin parameter :math:`a=J_\mathrm{BH}c/M_\mathrm{BH}^2 G`, where :math:`J_\mathrm{BH}` is its angular momentum. For theoretical reasons, its magnitude cannot grow above 1. It can be positive, representing prograde accretion, or negative, representing retrograde accretion. Any model for realistic AGN jets must include black hole spin since jet powers depend steeply on spin, and because including spin provides a well-defined direction for the jets to be launched in. The spin (angular momentum) of BHs is best represented through the dimensionlesss spin parameter :math:`a=J_\mathrm{BH}c/M_\mathrm{BH}^2 G`, where :math:`J_\mathrm{BH}` is its angular momentum. For theoretical reasons, its magnitude cannot grow above 1. It can be positive, representing prograde accretion, or negative, representing retrograde accretion.
Jet powers, in addition to depending on spin, also depend on which accretion state the black hole is in. We refer to these states by the shape of the accretion disk that surrounds the BH. We include three accretion states: the thick (or advection-dominated accretion flow; ADAF), thin (standard) and slim disk. Our main reference points for these disks are the following papers: `Shakura & Sunyaev (1973) <https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract>`_, `Narayan & Yi (1994) <https://ui.adsabs.harvard.edu/abs/1994ApJ...428L..13N/abstract>`_ and `Wang & Zhou. (1999) <https://ui.adsabs.harvard.edu/abs/1999ApJ...516..420W/abstract>`_, respectively. Jet powers, in addition to depending on spin, also depend on which accretion state the black hole is in. We refer to these states by the shape of the accretion disk that surrounds the BH. We include three accretion states: the thick (or advection-dominated accretion flow; ADAF), thin (standard) and slim disk (super-Eddington accretion). Our main reference points for these disks are the following papers: `Shakura & Sunyaev (1973) <https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract>`_, `Narayan & Yi (1994) <https://ui.adsabs.harvard.edu/abs/1994ApJ...428L..13N/abstract>`_ and `Wang & Zhou. (1999) <https://ui.adsabs.harvard.edu/abs/1999ApJ...516..420W/abstract>`_, respectively.
The thick disk appears at low accretion rates, has very strong jets and is inefficient at spinning up the black hole. The thin disk, appearing at intermediate accretion rates, typically has weak jets, strong radiation and efficiently spins up the black hole. The slim disk, corresponding to super-Eddington accretion, has features of both: in terms of geometry, orbits and angular momentum, it is similar to the thick disk. It is optically thin, leading to strong radiation. However, it also has strong jets. We assume that each subgrid accretion disk launches jets and radiates at the same time, regardless of the type it is. However, we use expressions for the jet and radiative efficiencies that depend on the type of the disk, and which are physically motivated. The thick disk appears at low accretion rates, has very strong jets and is inefficient at spinning up the black hole. The thin disk, appearing at intermediate accretion rates, typically has weak jets, strong radiation and efficiently spins up the black hole. The slim disk, corresponding to super-Eddington accretion, has features of both: in terms of geometry, orbits and angular momentum, it is similar to the thick disk. It is optically thin, leading to strong radiation. However, it also has strong jets. We assume that each subgrid accretion disk launches jets and radiates at the same time, regardless of the type it is. However, we use expressions for the jet and radiative efficiencies that depend on the type of the disk, and which are physically motivated.
...@@ -24,42 +24,63 @@ Transitions from one accretion state to another ...@@ -24,42 +24,63 @@ Transitions from one accretion state to another
:figclass: align-center :figclass: align-center
:alt: Accretion regimes :alt: Accretion regimes
The type of accretion disk surrounding the BHs depending on their accretion rates and spins. The transition between the thick and thin disk is calculated assuming the viscosity parameter :math:`\alpha=0.2`, while the transition from thin to slim disk is assumed to occur when the latter is :math:`F=0.5` times as radiatively efficienct as the former. The type of accretion disk surrounding the BHs depending on their accretion rates and spins.
The state of the subgrid accretion disk depends mostly on the Eddington fraction, i.e. the (dimensionless) accretion rate of the BH in units of the Eddington accretion rate, which we denote as :math:`\dot{m}`. We assume that the subgrid accretion disk is thick for :math:`\dot{m}<0.03`, based on observations (`Russell et al. 2013 <https://ui.adsabs.harvard.edu/abs/2013MNRAS.432..530R/abstract>`_). This also allows us to constrain the value of one of the main parameters in any accretion model: the viscosity parameter :math:`\alpha` (which is related to the kinematic viscosity :math:`\nu`and sound speed :math:`c_\mathrm{s}` through :math:`\nu=\alpha c_\mathrm{s}H`, with :math:`c_\mathrm{s}` the sound speed and :math:`H` the disk half-thickness). Numerical calculations suggest that thick disks are present for :math:`\dot{m}<0.4\alpha^2` (`Yuan & Narayan 2014 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_), and this agrees with observations if :math:`\alpha=0.25-0.3`. These values agree very well with more direct observational estimates, which suggest :math:`\alpha=0.2-0.3` (`Martin et al. 2019 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_). The state of the subgrid accretion disk depends mostly on the Eddington fraction, i.e. the (dimensionless) accretion rate of the BH in units of the Eddington accretion rate, which we denote as :math:`f_\mathrm{Edd}`. We assume that the subgrid accretion disk is thick for :math:`f_\mathrm{Edd}<f_\mathrm{Edd,crit,thick}`, where :math:`f_\mathrm{Edd,crit,thick}\approx0.01-0.03` (based on observations; `Russell et al. 2013 <https://ui.adsabs.harvard.edu/abs/2013MNRAS.432..530R/abstract>`_) is a free parameter in the model. The accretion disc is assumed to be slim for :math:`f_\mathrm{Edd}>1`.
The transition from the thin to the slim disk should occur around :math:`\dot{m}\approx 1`. However, the exact physics of this transition is not well understood. There is likely some spin dependence of the critical accretion rate, due to different radiative physics depending on spin. One of the main properties of slim disks is that they are less radiatively efficient than thin disks (`Sadowski et al. 2014 <https://ui.adsabs.harvard.edu/abs/2014MNRAS.439..503S/abstract>`_). We thus assume that the transition occurs when the radiative efficiency of a slim disk, :math:`\epsilon_\mathrm{r,SD}`, falls below some fraction of the radiative efficiency of a thin disk, :math:`\epsilon_\mathrm{r,TD}`. We quantify this as :math:`\epsilon_\mathrm{r,SD}<F\epsilon_\mathrm{r,SD}`, with :math:`F\approx 0.5` a free parameter. We give the expressions for both of the efficiencies below. Accretion efficiencies
-----------------------------------------------
Jet efficiencies Our model requires the usage of accretion efficiencies, minimally in the thick disk regime. These accretion efficiencies arise due to winds that take away most of the accreting mass as it falls towards the BH. We supress the large-scale accretion rate (e.g. the Bondi rate), :math:`\dot{M}_\mathrm{BH,0}:` such that the net accretion rate is equal to
----------------
The jet efficiency is related to the jet power through :math:`\epsilon_\mathrm{j}=P_\mathrm{j}/\dot{M}_\mathrm{BH,0}c^2`, where :math:`\dot{M}_\mathrm{BH,0}` is the accretion rate measured in the simulation, e.g. the Bondi rate). We use the formula for the jet efficiency based on general-relativistic, magneto-hydrodynamical (GRMHD) simulations by `Tchekhovskoy et al. (2010) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_: .. math::
\dot{M}_\mathrm{BH} = (1 - \epsilon_\mathrm{rad} - \epsilon_\mathrm{wind} - \epsilon_\mathrm{jet})\epsilon_\mathrm{acc}\dot{M}_\mathrm{BH,0},
where the terms in parenthesis are feedback efficiencies (discussed below), defined as :math:`\epsilon_i=P_i/\dot{M}_\mathrm{BH}c^2`, while :math:`\epsilon_\mathrm{acc}` is the accretion efficiency.
In the thick and slim disk, we allow options where the accretion efficiencies are free parameters in the model, to be tuned somehow (in practice, for the thick disc efficiency, to the local AGN bolometric luminosity function). We also allow a more complex scaling for the thick disk, based on GRMHD simulations (e.g. `Cho et al. 2024 <https://arxiv.org/abs/2405.13887>`_). These simulations show that the accretion efficiency in thick disks can be calculated as
.. math:: .. math::
\epsilon_\mathrm{j}=\frac{\kappa}{4\pi}\bigg(\frac{H/R}{0.3}\bigg)^\eta \phi_\mathrm{BH}^2\Omega_\mathrm{BH}^2\big(1+1.38\Omega_\mathrm{BH}^2-9.2\Omega_\mathrm{BH}^4\big), \epsilon_\mathrm{acc} = \bigg(\frac{R_0}{R_\mathrm{thick}}\bigg)^s,
where :math:`R_0\approx5-10` (in units of :math:`R_\mathrm{G}`), :math:`R_\mathrm{thick}` is the radius of the hot accretion flow and :math:`s=0.5` in recent such simulations. The radius :math:`R_\mathrm{thick}` is the same as the Bondi radius if the accretion rate is very low. However, at high accretion rates (but still below :math:`f_\mathrm{Edd}\approx0.01`), the accretion disk may be thin at large distances and thick at smaller ones, with the transition occuring at some radius :math:`R_\mathrm{tr}`. In this case, the winds (and mass loading) operate only at smaller radii. Given these considerations, we write
where :math:`\kappa\approx0.05` is a numerical factor which depends on the initial geometry of the magnetic field, :math:`\phi_\mathrm{BH}` is the dimensionless magnetic flux threading the horizon (see original paper for precise definition), and :math:`\Omega_\mathrm{BH}=a/2r_\mathrm{H}` is the (dimensionless) angular velocity of the black hole event horizon. Here, :math:`r_\mathrm{H}=1+\sqrt{1-a^2}` is the radius of the horizon in units of the gravitational radius :math:`R_\mathrm{G}=M_\mathrm{BH}G/c^2`. The formula above, for the jet efficiency, agrees very well with the results from higher-resolution simulations performed by `Narayan et al. (2021) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_, who provide the following fit for the magnetic flux as a function of spin: .. math::
R_\mathrm{thick} = \min(R_\mathrm{B},R_\mathrm{tr}),
and we choose to parametrize :math:`R_\mathrm{tr}` based on original calculations by `Narayan & Yi (1995) <https://ui.adsabs.harvard.edu/abs/1995ApJ...452..710N/abstract>`_, who found the transition radius as the radius where half of the energy gets radiated away, and half advected inwards. Their formula can be written in the form
.. math:: .. math::
\phi_\mathrm{BH}(a)=-20.2a^3-14.9a^2+34a+52.6. R_\mathrm{tr} = R_1 \bigg(\frac{0.01}{f_\mathrm{Edd}}\bigg)^2,
The `Tchekhovskoy et al. (2010) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_ jet efficiency depends very steeply on spin (:math:`\epsilon_\mathrm{j}\propto a^2` for small spin and :math:`\epsilon_\mathrm{j}\propto a^6` near :math:`a=1`). It can reach values above 100 per cent for large spins, and is also different (weaker) for negative spins. where :math:`R_1` is some normalisation radius that depends strongly on the assumed value of the accretion disk viscosity parameter :math:`\alpha`. We instead leave :math:`R_1` as a free parameter in our formulation. Note that at moderate Eddington ratios, where the Bondi radius is not the limiting radius (i.e. where a thin disk component exists outside the thick disk), we may write the accretion efficiency as:
The dependence of the jet efficiency on the type of accretion disk comes through the factor that depends on the aspect ratio :math:`H/R`, since accretion disks differ in this quantity. Theoretical, self-similar models of thick disks suggest :math:`H/R=0.5` (`Narayan & Yi 1995b <https://ui.adsabs.harvard.edu/abs/1995ApJ...452..710N/abstract>`_), but we instead take :math:`H/R=0.3`, more in line with simulations. For slim disks, which have received less attention in simulations, we assume the value :math:`H/R=1/(2\sqrt{5})\approx 0.2` (based on the self-similar model by `Wang & Zhou. 1999 <https://ui.adsabs.harvard.edu/abs/1999ApJ...516..420W/abstract>`_). .. math::
\epsilon_\mathrm{acc} = \sqrt{\frac{R_0}{R_1}}\bigg(\frac{f_\mathrm{Edd}}{0.01}\bigg),
where we have assumed :math:`s=0.5`.
Jet efficiencies
----------------
Thin disks are, not surprisingly, much thinner. The value of :math:`H/R` in this regime is not a constant, but rather depends on the BH mass and accretion rate, slightly on radius and also on the viscosity parameter :math:`\alpha`. Thin disks have three different regions in the `Shakura & Sunyaev (1973) <https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract>`_ model. For simplicity, we model the whole disk as being represented with only one region. In region a), the innermost one, radiation dominates over gas pressure. It is typically very small or doesn't exist at all, so we disregard it as a possibility. In regions b) and c), gas pressure dominates over radiation pressure. In b), electrons dominate in the opacity, while in c), free-free absorption dominates. We leave both regions as a possibility, and leave the choice as a free parameter in the model (not likely to lead to large differences in galaxy/BH evolution). The expressions for the aspect ratio in these regions are The jet efficiency is related to the jet power through :math:`\epsilon_\mathrm{j}=P_\mathrm{j}/\dot{M}_\mathrm{BH,0}c^2`, where :math:`\dot{M}_\mathrm{BH,0}` is the accretion rate measured in the simulation, e.g. the Bondi rate). We use the formula for the jet efficiency based on general-relativistic, magneto-hydrodynamical (GRMHD) simulations by `Tchekhovskoy et al. (2010) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_:
.. math:: .. math::
\bigg(\frac{H}{R}\bigg)_\mathrm{TD,b} = 1.25\times10^{-3} \alpha^{-1/10}\dot{m}^{1/5}\bigg(\frac{M_\mathrm{BH}}{10^8\hspace{0.5mm}\mathrm{M}_\odot}\bigg)^{-1/10}\bigg(\frac{R}{2R_\mathrm{G}}\bigg)^{1/20} \epsilon_\mathrm{j}=\frac{\kappa}{4\pi} \phi_\mathrm{BH}^2\Omega_\mathrm{BH}^2\big(1+1.38\Omega_\mathrm{BH}^2-9.2\Omega_\mathrm{BH}^4\big),
in region b) and where :math:`\kappa\approx0.05` is a numerical factor which depends on the initial geometry of the magnetic field, :math:`\phi_\mathrm{BH}` is the dimensionless magnetic flux threading the horizon (see original paper for precise definition), and :math:`\Omega_\mathrm{BH}=a/2r_\mathrm{H}` is the (dimensionless) angular velocity of the black hole event horizon. Here, :math:`r_\mathrm{H}=1+\sqrt{1-a^2}` is the radius of the horizon in units of the gravitational radius :math:`R_\mathrm{G}=M_\mathrm{BH}G/c^2`. The formula above, for the jet efficiency, agrees very well with the results from higher-resolution simulations performed by `Narayan et al. (2021) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_, who provide the following fit for the magnetic flux as a function of spin:
.. math:: .. math::
\bigg(\frac{H}{R}\bigg)_\mathrm{TD,c} = 1.15\times10^{-3} \alpha^{-1/10}\dot{m}^{3/20}\bigg(\frac{M_\mathrm{BH}}{10^8\hspace{0.5mm}\mathrm{M}_\odot}\bigg)^{-1/10}\bigg(\frac{R}{2R_\mathrm{G}}\bigg)^{1/8} \phi_\mathrm{BH,MAD}(a)=-20.2a^3-14.9a^2+34a+52.6.
The `Tchekhovskoy et al. (2010) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_ jet efficiency depends very steeply on spin (:math:`\epsilon_\mathrm{j}\propto a^2` for small spin and :math:`\epsilon_\mathrm{j}\propto a^6` near :math:`a=1`). It can reach values above 100 per cent for large spins, and is also different (weaker) for negative spins.
in region c). The dependence of the jet efficiency on the type of accretion disk is encoded in the fact that thick disks are thought to be in a magnetically-arred state (so-called MAD. see `Narayan et al. 2003 <https://ui.adsabs.harvard.edu/abs/2003PASJ...55L..69N/abstract>`_), while thin disks are likely not, because they do not feature strong advection. The slim disk, on the other hand, is thought to be similar to the thick disk in terms of advection, and thus probably in terms of jet powers. Recent simulations by `Ricarte et al. (2023) <https://ui.adsabs.harvard.edu/abs/2023ApJ...954L..22R/abstract>`_ have found an increase of :math:`\phi_\mathrm{BH}` in the thin and slim disk regime as the Eddington ratio increases, and they parametrise this increase as
The jet efficiency also depends on the slope :math:`\eta`. Classical jet theory (`Meier 2001 <https://ui.adsabs.harvard.edu/abs/2001Sci...291...84M/abstract>`_) suggests that jet powers depend on the aspect ratio linearly, so :math:`\eta=1`. This is also in line with some simulations finding a reduction in jet efficiencies with the aspect ratio (e.g. `Tchekhovskoy et al. 2014 <https://ui.adsabs.harvard.edu/abs/2014MNRAS.437.2744T/abstract>`_). In this scenario, jets launched from thin disks are of order :math:`\approx100` times less powerful than those launched from thick disks. On the other hand, some simulations of thin disks have found jet efficiencies similar to thick disk ones (e.g. `Liska et al. 2019 <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487..550L/abstract>`_), which is supported by observations of blazars (`Ghisellini et al. 2014 <https://ui.adsabs.harvard.edu/abs/2014Natur.515..376G/abstract>`_). In this picture, thin jets are approximately as efficient as thick disk ones, which can in our case be implemented as :math:`\eta=0`. The reality is likely to be somwhere in between. Note that the choice of :math:`\eta` likely has a strong impact on the evolution of galaxies and BHs; our default choice is the classical picture in which :math:`\eta=1`. .. math::
\phi_\mathrm{BH,thin,slim} = \frac{(f_\mathrm{Edd}/1.88)^{1.29}}{1+(f_\mathrm{Edd}/1.88)^{1.29}}\phi_\mathrm{BH,MAD}.
The magnetic flux eventually saturates (at very high :math:`f_\mathrm{Edd}`) at the same value as that reached in the thick disc; :math:`\phi_\mathrm{BH,MAD}`.
.. figure:: efficiencies.png .. figure:: efficiencies.png
:width: 1200px :width: 1200px
...@@ -69,8 +90,8 @@ The jet efficiency also depends on the slope :math:`\eta`. Classical jet theory ...@@ -69,8 +90,8 @@ The jet efficiency also depends on the slope :math:`\eta`. Classical jet theory
Feedback efficiencies (jet - blue, radiation - red) for all three accretion disk types. Shaded regions represent likely ranges of efficiencies (where the efficiencies depend on mass and/or accretion rate). The thin disk jet efficiencies were computed assuming the slope of the efficiency vs. aspect ratio relation is :math:`\eta=1`, and the aspect ratios were computed for region b) of the Shakura & Sunyaev solution. Radiative efficiencies in the thick disk were computed assuming the electron heating parameter :math:`\delta=0.2`. Feedback efficiencies (jet - blue, radiation - red) for all three accretion disk types. Shaded regions represent likely ranges of efficiencies (where the efficiencies depend on mass and/or accretion rate). The thin disk jet efficiencies were computed assuming the slope of the efficiency vs. aspect ratio relation is :math:`\eta=1`, and the aspect ratios were computed for region b) of the Shakura & Sunyaev solution. Radiative efficiencies in the thick disk were computed assuming the electron heating parameter :math:`\delta=0.2`.
Radiative efficiencies Radiative/wind efficiencies
---------------------- ---------------------------
In the EAGLE and COLIBRE models, all subgrid accretion disks are effectively thin, and the BH is always assumed to be in this regime. In our model, the radiative efficiency (defined in an analagous way to the jet efficiency, but using the luminosity) is no longer fixed at a value of order :math:`10` per cent. Instead, we use spin-dependant formulas that vary with the type of disk. In the thin disk, the radiative efficiency :math:`\epsilon_\mathrm{r,TD}` is related to the binding energy at the innermost stable circular orbit (ISCO) and is given by In the EAGLE and COLIBRE models, all subgrid accretion disks are effectively thin, and the BH is always assumed to be in this regime. In our model, the radiative efficiency (defined in an analagous way to the jet efficiency, but using the luminosity) is no longer fixed at a value of order :math:`10` per cent. Instead, we use spin-dependant formulas that vary with the type of disk. In the thin disk, the radiative efficiency :math:`\epsilon_\mathrm{r,TD}` is related to the binding energy at the innermost stable circular orbit (ISCO) and is given by
...@@ -79,7 +100,7 @@ In the EAGLE and COLIBRE models, all subgrid accretion disks are effectively thi ...@@ -79,7 +100,7 @@ In the EAGLE and COLIBRE models, all subgrid accretion disks are effectively thi
Here, :math:`r_\mathrm{ISCO}` is the radius of the ISCO in gravitational radii (see e.g. appendix B of `Fiacconi et al. 2018 <https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.3807F/abstract>`_ for an expression giving the spin dependence). The radiative efficiency of the thin disk grows slowly from its minimum value of :math:`\approx4` per cent for :math:`a=-1` to :math:`\approx5.5` per cent for :math:`a=0`. For positive spins it grows more steeply; it is :math:`10` per cent by :math:`a=0.65`. Beyond that the dependence steepens even further, with values of :math:`20`, :math:`30` and :math:`40` per cent reached at :math:`a=0.95`, :math:`a=0.997` and :math:`a=1`, respectively. Here, :math:`r_\mathrm{ISCO}` is the radius of the ISCO in gravitational radii (see e.g. appendix B of `Fiacconi et al. 2018 <https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.3807F/abstract>`_ for an expression giving the spin dependence). The radiative efficiency of the thin disk grows slowly from its minimum value of :math:`\approx4` per cent for :math:`a=-1` to :math:`\approx5.5` per cent for :math:`a=0`. For positive spins it grows more steeply; it is :math:`10` per cent by :math:`a=0.65`. Beyond that the dependence steepens even further, with values of :math:`20`, :math:`30` and :math:`40` per cent reached at :math:`a=0.95`, :math:`a=0.997` and :math:`a=1`, respectively.
In the thick disk regime, radiative efficiencies are lower by a factor :math:`\approx100` than jet efficiencies. The formulas we use are based on results by `Mahadevan (1997) <https://ui.adsabs.harvard.edu/abs/1997ApJ...477..585M/abstract>`_, who studied cooling processes of electrons (which dominate in the radiation) in the context of the original thick disc solution. They found two different regimes: for :math:`\dot{m}<\dot{m}_\mathrm{crit,visc}`, viscous heating dominates the heating of electrons, whereas for :math:`\dot{m}_\mathrm{crit,visc}<\dot{m}<\dot{m}_\mathrm{crit,ADAF}`, it is dominated by ion-electron heating. Here, :math:`\dot{m}_\mathrm{crit,visc}` is the transitional value between the two thick disc (ADAF) regimes, and :math:`\dot{m}_\mathrm{crit,ADAF}=0.4\alpha^2` is the transitional accretion rate which separates thin and thick discs. The radiative efficiency in the viscous heating regime is given by In the thick disk regime, radiative efficiencies are lower by a factor :math:`\approx100` than jet efficiencies. The formulas we use are based on results by `Mahadevan (1997) <https://ui.adsabs.harvard.edu/abs/1997ApJ...477..585M/abstract>`_, who studied cooling processes of electrons (which dominate in the radiation) in the context of the original thick disc solution. They found two different regimes: for :math:`f_\mathrm{Edd}<f_\mathrm{Edd,crit,visc}`, viscous heating dominates the heating of electrons, whereas for :math:`f_\mathrm{Edd,crit,visc}<f_\mathrm{Edd}<f_\mathrm{Edd,crit,ADAF}`, it is dominated by ion-electron heating. Here, :math:`f_\mathrm{Edd,crit,visc}` is the transitional value between the two thick disc (ADAF) regimes, and :math:`f_\mathrm{Edd,crit,ADAF}=0.4\alpha^2` is the transitional accretion rate which separates thin and thick discs. The radiative efficiency in the viscous heating regime is given by
.. math:: .. math::
\epsilon_\mathrm{r,ADAF}=0.0002\epsilon_\mathrm{r,TD}\bigg(\frac{\delta_\mathrm{ADAF}}{0.0005}\bigg)\bigg(\frac{1-\beta}{0.5}\bigg)\bigg(\frac{6}{r_\mathrm{ISCO}}\bigg), \epsilon_\mathrm{r,ADAF}=0.0002\epsilon_\mathrm{r,TD}\bigg(\frac{\delta_\mathrm{ADAF}}{0.0005}\bigg)\bigg(\frac{1-\beta}{0.5}\bigg)\bigg(\frac{6}{r_\mathrm{ISCO}}\bigg),
...@@ -87,30 +108,32 @@ In the thick disk regime, radiative efficiencies are lower by a factor :math:`\a ...@@ -87,30 +108,32 @@ In the thick disk regime, radiative efficiencies are lower by a factor :math:`\a
while in the ion-heating regime it is given by while in the ion-heating regime it is given by
.. math:: .. math::
\epsilon_\mathrm{r,ADAF}=0.2\epsilon_\mathrm{r,TD}\bigg(\frac{\dot{m}}{\alpha^2}\bigg)\bigg(\frac{\beta}{0.5}\bigg)\bigg(\frac{6}{r_\mathrm{ISCO}}\bigg). \epsilon_\mathrm{r,ADAF}=0.2\epsilon_\mathrm{r,TD}\bigg(\frac{f_\mathrm{Edd}}{\alpha^2}\bigg)\bigg(\frac{\beta}{0.5}\bigg)\bigg(\frac{6}{r_\mathrm{ISCO}}\bigg).
Here, :math:`\beta` is the ratio of gas pressure and total pressure (which includes the magnetic pressure). `Yuan & Narayan (2014) <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_ define a somewhat different parameter, :math:`\beta_\mathrm{ADAF}`, as the ratio of gas pressure and magnetic pressure. The two parameters are related by :math:`\beta=\beta_\mathrm{ADAF}/(1+\beta_\mathrm{ADAF})`. :math:`\beta_\mathrm{ADAF}` is not an independent parameter; many simulations have found that :math:`\alpha\beta_\mathrm{ADAF}\approx0.5` (e.g. `Begelman et al. 2021 <https://ui.adsabs.harvard.edu/abs/2022MNRAS.511.2040B/abstract>`_, see also `Yuan & Narayan 2014 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_ for a review), which we adopt. :math:`\delta_\mathrm{ADAF}` represents the fraction of viscous energy transferred to the electrons, and is constrained in theoretical studies between 0.1 and 0.5 (`Yuan & Narayan 2014 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_, `Sharma et al. 2007 <https://ui.adsabs.harvard.edu/abs/2007ApJ...667..714S/abstract>`_). Observations imply a value close to 0.2 (`Yuan et al. 2003 <https://ui.adsabs.harvard.edu/abs/2003ApJ...598..301Y/abstract>`_, `Liu & Wu 2013 <https://ui.adsabs.harvard.edu/abs/2013ApJ...764...17L/abstract>`_). The critical accretion rate between the two thick disc regimes can be found by ensuring that both formulas presented above yield the same radiative efficiency (at that accretion rate). This gives an accretion rate equal to Here, :math:`\beta` is the ratio of gas pressure and total pressure (which includes the magnetic pressure). `Yuan & Narayan (2014) <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_ define a somewhat different parameter, :math:`\beta_\mathrm{ADAF}`, as the ratio of gas pressure and magnetic pressure. The two parameters are related by :math:`\beta=\beta_\mathrm{ADAF}/(1+\beta_\mathrm{ADAF})`. :math:`\beta_\mathrm{ADAF}` is not an independent parameter; many simulations have found that :math:`\alpha\beta_\mathrm{ADAF}\approx0.5` (e.g. `Begelman et al. 2021 <https://ui.adsabs.harvard.edu/abs/2022MNRAS.511.2040B/abstract>`_, see also `Yuan & Narayan 2014 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_ for a review), which we adopt. :math:`\delta_\mathrm{ADAF}` represents the fraction of viscous energy transferred to the electrons, and is constrained in theoretical studies between 0.1 and 0.5 (`Yuan & Narayan 2014 <https://ui.adsabs.harvard.edu/abs/2014ARA%26A..52..529Y/abstract>`_, `Sharma et al. 2007 <https://ui.adsabs.harvard.edu/abs/2007ApJ...667..714S/abstract>`_). Observations imply a value close to 0.2 (`Yuan et al. 2003 <https://ui.adsabs.harvard.edu/abs/2003ApJ...598..301Y/abstract>`_, `Liu & Wu 2013 <https://ui.adsabs.harvard.edu/abs/2013ApJ...764...17L/abstract>`_). The critical accretion rate between the two thick disc regimes can be found by ensuring that both formulas presented above yield the same radiative efficiency (at that accretion rate). This gives an accretion rate equal to
.. math:: .. math::
\dot{m}_\mathrm{crit,visc}=0.0002\bigg(\frac{\delta_\mathrm{ADAF}}{0.0005}\bigg)\bigg(\frac{1-\beta}{\beta}\bigg)\alpha^2. f_\mathrm{Edd,crit,visc}=0.0002\bigg(\frac{\delta_\mathrm{ADAF}}{0.0005}\bigg)\bigg(\frac{1-\beta}{\beta}\bigg)\alpha^2.
For slim disks we take the radiative efficiency based on GRMHD simulations of super-Eddington accretion (for various BH spins) performed by `Sadowski et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014MNRAS.439..503S/abstract>`_. `Madau et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014ApJ...784L..38M/abstract>`_ found the following fitting function which represents the `Sadowski et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014MNRAS.439..503S/abstract>`_ results: For slim disks we take the radiative efficiency based on GRMHD simulations of super-Eddington accretion (for various BH spins) performed by `Sadowski et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014MNRAS.439..503S/abstract>`_. `Madau et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014ApJ...784L..38M/abstract>`_ found the following fitting function which represents the `Sadowski et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014MNRAS.439..503S/abstract>`_ results:
.. math:: .. math::
\epsilon_\mathrm{r,SD}=\frac{0.1}{\dot{m}}A(a)\bigg( \frac{0.985}{1.6/\dot{m}+B(a)}+\frac{0.015}{1.6/\dot{m}+C(a)}\bigg), \epsilon_\mathrm{r,SD}=\frac{0.1}{f_\mathrm{Edd}}A(a)\bigg( \frac{0.985}{1.6/f_\mathrm{Edd}+B(a)}+\frac{0.015}{1.6/f_\mathrm{Edd}+C(a)}\bigg),
where the three spin-dependant functions are given by :math:`A(a)=(0.9663-0.9292a)^{-0.5639}`, :math:`B(a)=(4.627-4.445a)^{-0.5524}` and :math:`C(a)=(827.3-718.1a)^{-0.7060}`. The radiative efficiency of slim disks, based on this formula, matches the thin disk radiative efficiency (given at the beginning of the section) at low accretion rates. At high accretion rates (:math:`\dot{m}\gtrapprox1`, but depending on spin), the radiative efficiency drops. These two formulas are used to decide when a disk transitions from thin to slim. where the three spin-dependant functions are given by :math:`A(a)=(0.9663-0.9292a)^{-0.5639}`, :math:`B(a)=(4.627-4.445a)^{-0.5524}` and :math:`C(a)=(827.3-718.1a)^{-0.7060}`. The radiative efficiency of slim disks, based on this formula, matches the thin disk radiative efficiency (given at the beginning of the section) at low accretion rates. At high accretion rates (:math:`f_\mathrm{Edd}\gtrapprox1`, but depending on spin), the radiative efficiency drops.
Evolution of the black hole spin magnitude The thin disc radiative efficiency is used to source feedback in the simulations. In the thin disk regime, a fraction :math:`\epsilon_\mathrm{f}\approx0.1` of all of the radiation released by black holes couples to the gas in the form of thermal energy. In the thick and slim disk, we do not use radiation to source feedback. We do, however, assume that winds launched from the accretion disk are present in these two states. In the thick disk, winds are thought to be launched on account of a combination of gas pressure and MHD effects. We use the formula from `Sadowski et al. (2013) <https://ui.adsabs.harvard.edu/abs/2013MNRAS.436.3856S/abstract>`_:
------------------------------------------
.. figure:: spec_ang_mom.png .. math::
:width: 600px \epsilon_\mathrm{wind,thick} = 0.005\bigg[1+0.3\bigg(\frac{\phi_\mathrm{BH,MAD}}{50}\bigg)\bigg(\frac{\Omega_\mathrm{H}}{0.2}\bigg) \bigg].
:align: center
:figclass: align-center For the slim disk, we again use results from `Ricarte et al. (2023) <https://ui.adsabs.harvard.edu/abs/2023ApJ...954L..22R/abstract>`_, as we did for the jet efficiency. We use their total MHD efficiency and subtract from that the analytical jet efficiency as given by the formula we use as a function of spin and magnetic flux. We then found a simple fitting function to the remaining efficiency, representing the wind:
:alt: Angular momenta
Dimensionless pecific angular momentum of the thin disk at the innermost stable circular orbit (ISCO, solid red line), compared with the specific angular momentum at the inner radius (the event horizon) for advection-dominated flows (the thick and slim disk) for a few values of the viscosity parameter :math:`\alpha`. The dashed red line shows that the latter can be approximated as :math:`45` per cent of the former. .. math::
\epsilon_\mathrm{wind,slim} = 0.065\bigg[1+\bigg(\frac{\phi_\mathrm{BH,thin,slim}}{50}\bigg)^2\bigg] \big(1+\Omega_\mathrm{H}-8\Omega_\mathrm{H}^2\big).
Evolution of the black hole spin magnitude
------------------------------------------
The BH spin (or angular momentum) is, naturally, a vector. However, due to Lense-thirring torques (we discuss these in more detail below), the accretion disk is always either aligned or counteraligned with the rotational axis of the black hole. This means that almost all relevant quantities, such as the efficiencies discussed above, can be expressed as depending only on the magnitude of spin, but also allowing for a negative sign to account for counteraligned disks (retrograde accretion). This is also true for the evolution of the magnitude of spin. The BH spin (or angular momentum) is, naturally, a vector. However, due to Lense-thirring torques (we discuss these in more detail below), the accretion disk is always either aligned or counteraligned with the rotational axis of the black hole. This means that almost all relevant quantities, such as the efficiencies discussed above, can be expressed as depending only on the magnitude of spin, but also allowing for a negative sign to account for counteraligned disks (retrograde accretion). This is also true for the evolution of the magnitude of spin.
...@@ -119,16 +142,34 @@ In the absence of jet spindown, the evolution of angular momentum is given simpl ...@@ -119,16 +142,34 @@ In the absence of jet spindown, the evolution of angular momentum is given simpl
.. math:: .. math::
\frac{\mathrm{d}a}{\mathrm{d}\ln M_\mathrm{BH,0}}=\ell_\mathrm{in}-2a e_\mathrm{in}, \frac{\mathrm{d}a}{\mathrm{d}\ln M_\mathrm{BH,0}}=\ell_\mathrm{in}-2a e_\mathrm{in},
where :math:`\ell_\mathrm{in}` is the specific angular momentum in units where :math:`G` and :math:`c` are equal to unity, and :math:`\mathrm{d}\ln M_\mathrm{BH,0}=\mathrm{d}M_\mathrm{BH,0}/M_\mathrm{BH}` is the logarithmic change in mass, not including losses due to radiation (`Fanidakis et al. 2011 <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410...53F/abstract>`_). The specific binding energy can be related to the radiative efficiency through :math:`e_\mathrm{in}=1-\epsilon_\mathrm{r}` for all three accretion states (for the thick disc, the radiative efficiency is negligible for this application). where :math:`\ell_\mathrm{in}` is the specific angular momentum in units where :math:`G` and :math:`c` are equal to unity, and :math:`\mathrm{d}\ln M_\mathrm{BH,0}=\mathrm{d}M_\mathrm{BH,0}/M_\mathrm{BH}` is the logarithmic change in mass, not including losses due to radiation (`Fanidakis et al. 2011 <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410...53F/abstract>`_). The specific binding energy can be related to the radiative efficiency through :math:`e_\mathrm{in}=1-\epsilon_\mathrm{r}` for all three accretion states (for the thick disc, the radiative efficiency is negligible for this application). All of the above quantities are evaluated at some inner radius beyond which gas orbits are unstable.
For the thin disc, the inner radius :math:`R_\mathrm{in}` can be taken to be the radius of the ISCO, since orbits should quickly degrade within it. We thus take :math:`\ell_\mathrm{in}` as the specific angular momentum at the ISCO for the thin disc (the expression for which is given in e.g. Appendix B of `Fiacconi et al. 2018 <https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.3807F/abstract>`_). For the thin disk, the spinup function (the equation shown above) is always positive, meaning that the BH will always be spun up to positive values. This means that the BH will be spun down if spin is negative, or spun up to an equilibrium value of :math:`a_\mathrm{eq}=1` if spin is positive. For advection-dominated disks (the thick and slim disk), we assume that :math:`\ell_\mathrm{in}` is :math:`45` per cent of the ISCO value, based on numerical GR calculations by `Popham & Gammie (1998) <https://ui.adsabs.harvard.edu/abs/1998ApJ...504..419P/abstract>`_. We base this assumption on fits of the `Popham & Gammie (1998) <https://ui.adsabs.harvard.edu/abs/1998ApJ...504..419P/abstract>`_ results done by `Benson & Babul (2009) <https://ui.adsabs.harvard.edu/abs/2009MNRAS.397.1302B/abstract>`_. We independently compared these fits to the ISCO values and found :math:`\ell_\mathrm{in}\approx0.45\ell_\mathrm{ISCO}` with no more than :math:`10` per cent error for all values of spin and relevant values of :math:`\alpha=0.1-0.3`. To be consistent with what we assumed for feedback efficiencies, we take results for the spinup/spindown function directly from GRMHD simulations. For the thick disc, we use the formula from `Narayan et al. (2021) <https://ui.adsabs.harvard.edu/abs/2010ApJ...711...50T/abstract>`_:
For the thick and slim disk, these lower specific angular momenta lead to a non-zero equilibrium spin value :math:`a_\mathrm{eq}<1`. If :math:`a>a_\mathrm{eq}`, the BH will be spun down due to frame-dragging and viscosity; the frame-dragging rotationally accelerates any accreting gas (on account of the BH angular momentum), while viscosity carries away some of that angular momentum. Including jets into the model leads to further spindown. The jet spindown term (to be added to the spinup equation above) can be derived as .. math::
\bigg(\frac{\mathrm{d}a}{\mathrm{d}M_\mathrm{BH,0}/M_\mathrm{BH}}\bigg)_\mathrm{thick}=0.45 - 12.53a - 7.8a^2 +9.44a^3 + 5.71a^4 -4.03a^5.
For the slim and thin disc, we use results from `Ricarte et al. (2023) <https://ui.adsabs.harvard.edu/abs/2023ApJ...954L..22R/abstract>`_, who find a fitting formula that smoothly interpolates between the thin disc regime without significant jet feedback (for :math:`f_\mathrm{Edd}` not close to super-Eddington values), and that where jet feedback essentially matches the thick disc (and so jet spindown should also be similar). Their formula takes the form
.. math:: .. math::
\bigg(\frac{\mathrm{d}a}{\mathrm{d}\ln M_\mathrm{BH,0}}\bigg)_\mathrm{j}=-\epsilon_\mathrm{j}(a)\frac{\sqrt{1-a^2}}{a}\bigg[\Big(\sqrt{1-a^2}+1 \Big)^2+a^2 \bigg] \bigg(\frac{\mathrm{d}a}{\mathrm{d}M_\mathrm{BH,0}/M_\mathrm{BH}}\bigg)_\mathrm{thin/slim}=s_\mathrm{HD} - s_\mathrm{EM},
(see `Benson & Babul 2009 <https://ui.adsabs.harvard.edu/abs/2009MNRAS.397.1302B/abstract>`_ for a derivation, which we have independently verified). Including jet spindown leads to even lower equilibrium spin values; e.g. for the thick disk this is only :math:`a_\mathrm{eq}\approx0.25`. where the first term is a pure hydrodynamical term, while the second is an electromagnetic term. The first term is given by
.. math::
s_\mathrm{HD}=\frac{s_\mathrm{thin}+s_\mathrm{min}\xi}{1+\xi},
where :math:`\xi=0.017f_\mathrm{Edd}`, :math:`s_\mathrm{min}=0.86-1.94a` and :math:`s_\mathrm{thin}=\ell_\mathrm{ISCO}-2a e_\mathrm{ISCO}` is the spinup/spindown function of the 'pure' thin disc (with no outflows and outside the MAD regime), in which :math:`\ell_\mathrm{ISCO}` and :math:`e_\mathrm{ISCO}` are the (dimensionless) specific angular momentum and binding energy, respectively, at the ISCO. The EM term is given by
.. math::
s_\mathrm{EM}=\mathrm{sgn}(a)\epsilon_\mathrm{EM}\bigg(\frac{1}{k\Omega_\mathrm{H}}-2a\bigg),
where :math:`\epsilon_\mathrm{EM}` is the total (jet+wind) EM efficiency, and :math:`k` is given by
.. math::
k=\min(0.35,0.1+0.5a)
for positive spins :math:`a>0` and by :math:`k=0.23` for negative spins :math:`a<0`.
.. figure:: spinup.png .. figure:: spinup.png
:width: 1200px :width: 1200px
...@@ -136,7 +177,7 @@ For the thick and slim disk, these lower specific angular momenta lead to a non- ...@@ -136,7 +177,7 @@ For the thick and slim disk, these lower specific angular momenta lead to a non-
:figclass: align-center :figclass: align-center
:alt: Spinup/spindown function :alt: Spinup/spindown function
Spinup/spindown function (the rate of black hole spin evolution) as a function of spin for all three accretion disk types. Black lines show evolution with only accretion included, while blue lines show the total including jet spindown. These plots show that the thin disk is always spun up to to :math:`a_\mathrm{eq}=1`, even with jets (due to low jet efficiencies). The advection-dominated disks (thick and slim disk) are spun up to positive equilibrium values :math:`a_\mathrm{eq}<1`, or spun down to such an equilibrium value if :math:`a>a_\mathrm{eq}`. This is due to extraction of rotational energy from the BH by frame dragging and transport of the angular momentum to large distances through viscous forces. Including jet spindown pushes these equilibrium spins to even smaller values. Spinup/spindown function (the dimensionless rate of black hole spin evolution) as a function of spin for all three accretion disk types. For the thin and slim disk, we show several curves for different choices of the Eddington ratio.
Evolution of the black hole spin direction Evolution of the black hole spin direction
------------------------------------------ ------------------------------------------
...@@ -149,7 +190,7 @@ In all cases, Lense-Thirring torques are effective only within some radius :math ...@@ -149,7 +190,7 @@ In all cases, Lense-Thirring torques are effective only within some radius :math
In terms of the evolution of the spin direction, the main assumption of our model is as follows (see `King et al. 2005 <https://ui.adsabs.harvard.edu/abs/2005MNRAS.363...49K/abstract>`_ for the original argument, and additional discussions in e.g. `Fanidakis et al. 2011 <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410...53F/abstract>`_, `Fiacconi et al. 2018 <https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.3807F/abstract>`_ and `Griffin et al. 2019a <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487..198G/abstract>`_). All matter that flows through an accretion disk is aligned or counteraligned with the BH spin vector in the accretion process. Due to conservation of angular momentum, the spin vector itself also has to adjust to keep the total angular momentum conserved. In the process of consuming one warp mass :math:`M_\mathrm{warp}`, the direction of the BH spin vector is aligned to match the direction of the total angular momentum of the system comprising the BH and the disk out to the warp radius. The direction of the BH spin vector can then be determined from :math:`\vec{J}_\mathrm{warp}=\vec{J}_\mathrm{BH}+J_\mathrm{warp}\hat{J}_\mathrm{d}`, where :math:`\vec{J}_\mathrm{BH}` is the old BH angular momentum vector, and :math:`\hat{J}_\mathrm{d}` is the direction of the large-scale accretion disk (which we assume matches the direction of the angular momentum of the gas in the BH smoothing kernel). In terms of the evolution of the spin direction, the main assumption of our model is as follows (see `King et al. 2005 <https://ui.adsabs.harvard.edu/abs/2005MNRAS.363...49K/abstract>`_ for the original argument, and additional discussions in e.g. `Fanidakis et al. 2011 <https://ui.adsabs.harvard.edu/abs/2011MNRAS.410...53F/abstract>`_, `Fiacconi et al. 2018 <https://ui.adsabs.harvard.edu/abs/2018MNRAS.477.3807F/abstract>`_ and `Griffin et al. 2019a <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487..198G/abstract>`_). All matter that flows through an accretion disk is aligned or counteraligned with the BH spin vector in the accretion process. Due to conservation of angular momentum, the spin vector itself also has to adjust to keep the total angular momentum conserved. In the process of consuming one warp mass :math:`M_\mathrm{warp}`, the direction of the BH spin vector is aligned to match the direction of the total angular momentum of the system comprising the BH and the disk out to the warp radius. The direction of the BH spin vector can then be determined from :math:`\vec{J}_\mathrm{warp}=\vec{J}_\mathrm{BH}+J_\mathrm{warp}\hat{J}_\mathrm{d}`, where :math:`\vec{J}_\mathrm{BH}` is the old BH angular momentum vector, and :math:`\hat{J}_\mathrm{d}` is the direction of the large-scale accretion disk (which we assume matches the direction of the angular momentum of the gas in the BH smoothing kernel).
In practice, the BH will consume parcels of mass that differ from :math:`M_\mathrm{warp}`. We assume that any such parcel of mass :math:`\Delta M` (e.g. the mass to be consumed within a single time step) can be split up onto :math:`n=\Delta M / M_\mathrm{warp}` individual increments of accretion, so the total angular momentum of the system within that time step is :math:`\vec{J}_\mathrm{warp}=\vec{J}_\mathrm{BH}+n J_\mathrm{warp}\hat{J}_\mathrm{d}`, i.e. :math:`n` warp angular momenta are consumed, with an angular momentum of :math:`\Delta \vec{J}=n J_\mathrm{warp}\hat{J}_\mathrm{d}=(J_\mathrm{warp}/M_\mathrm{warp})\Delta M `. This can also be viewed as the BH consuming material with a specific angular momentum of :math:`L_\mathrm{warp}=J_\mathrm{warp}/M_\mathrm{warp}`. Note that this picture is only valid if the BH spin vector does not change much during this process (in both magnitude and direction), which can be ensured with wisely chosen time steps. In practice, the BH will consume parcels of mass that differ from :math:`M_\mathrm{warp}`. We assume that any such parcel of mass :math:`\Delta M` (e.g. the mass to be consumed within a single time step) can be split up onto :math:`n=\Delta M / M_\mathrm{warp}` individual increments of accretion, so the total angular momentum of the system within that time step is :math:`\vec{J}_\mathrm{warp}=\vec{J}_\mathrm{BH}+n J_\mathrm{warp}\hat{J}_\mathrm{d}`, i.e. :math:`n` warp angular momenta are consumed, with an angular momentum of :math:`\Delta \vec{J}=n J_\mathrm{warp}\hat{J}_\mathrm{d}=(J_\mathrm{warp}/M_\mathrm{warp})\Delta M`. This can also be viewed as the BH consuming material with a specific angular momentum of :math:`L_\mathrm{warp}=J_\mathrm{warp}/M_\mathrm{warp}`. Note that this picture is only valid if the BH spin vector does not change much during this process (in both magnitude and direction), which can be ensured with wisely chosen time steps.
Deciding whether accretion is prograde or retrograde Deciding whether accretion is prograde or retrograde
---------------------------------------------------- ----------------------------------------------------
...@@ -172,34 +213,34 @@ As mentioned already, Lense-Thirring torques have different effects depending on ...@@ -172,34 +213,34 @@ As mentioned already, Lense-Thirring torques have different effects depending on
(`Ogilvie 1999 <https://ui.adsabs.harvard.edu/abs/1999MNRAS.304..557O/abstract>`_, see also `Lodato et al. 2010 <https://ui.adsabs.harvard.edu/abs/2010MNRAS.405.1212L/abstract>`_ for a detailed discussion). We use the relation :math:`\dot{M}=3\pi\nu_1 \Sigma` to calculate :math:`\nu_1`, and therefore :math:`\nu_2`. The warp radius will depend on which region of the thin disc we assume, with each having its own expression for :math:`\Sigma`. In region b) of the `Shakura & Sunyaev (1973) <https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract>`_ thin disk, the surface density can be expressed as (`Ogilvie 1999 <https://ui.adsabs.harvard.edu/abs/1999MNRAS.304..557O/abstract>`_, see also `Lodato et al. 2010 <https://ui.adsabs.harvard.edu/abs/2010MNRAS.405.1212L/abstract>`_ for a detailed discussion). We use the relation :math:`\dot{M}=3\pi\nu_1 \Sigma` to calculate :math:`\nu_1`, and therefore :math:`\nu_2`. The warp radius will depend on which region of the thin disc we assume, with each having its own expression for :math:`\Sigma`. In region b) of the `Shakura & Sunyaev (1973) <https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract>`_ thin disk, the surface density can be expressed as
.. math:: .. math::
\Sigma_\mathrm{TD,b}=6.84 \times 10^{5} \mathrm{~g} \mathrm{~cm}^{-2} \alpha^{-4 / 5} \dot{m}^{3 / 5}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 8}\left(\frac{R}{R_{\mathrm{S}}}\right)^{-3 / 5}, \Sigma_\mathrm{TD,b}=6.84 \times 10^{5} \mathrm{~g} \mathrm{~cm}^{-2} \alpha^{-4 / 5} f_\mathrm{Edd}^{3 / 5}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 8}\left(\frac{R}{R_{\mathrm{S}}}\right)^{-3 / 5},
while in region c) we have while in region c) we have
.. math:: .. math::
\Sigma_\mathrm{TD,c}=3.41 \times 10^{4} \mathrm{~g} \mathrm{~cm}^{-2} \alpha^{-4 / 5} \dot{m}^{7/10}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 20}\left(\frac{R}{R_{\mathrm{S}}}\right)^{-3 / 4}. \Sigma_\mathrm{TD,c}=3.41 \times 10^{4} \mathrm{~g} \mathrm{~cm}^{-2} \alpha^{-4 / 5} f_\mathrm{Edd}^{7/10}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 20}\left(\frac{R}{R_{\mathrm{S}}}\right)^{-3 / 4}.
These relations lead to the following expressions for :math:`R_\mathrm{warp}`: These relations lead to the following expressions for :math:`R_\mathrm{warp}`:
.. math:: .. math::
R_{\text {warp,TD,b}}=3410 R_{S} a^{5 / 8} \xi^{-5/8}\alpha^{-1 / 2} \dot{m}^{-1 / 4}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 8} R_{\text {warp,TD,b}}=3410 R_{S} a^{5 / 8} \xi^{-5/8}\alpha^{-1 / 2} f_\mathrm{Edd}^{-1 / 4}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{1 / 8}
(in region b) and (in region b) and
.. math:: .. math::
R_\mathrm{warp,TD,c}=2629R_\mathrm{S}a^{4/7}\xi^{-4/7}\alpha^{-16/35}\dot{m}^{-6/35}\bigg(\frac{M_\mathrm{BH}}{10^8\hspace{0.5mm}\mathrm{M}_\odot} \bigg)^{4/35}, R_\mathrm{warp,TD,c}=2629R_\mathrm{S}a^{4/7}\xi^{-4/7}\alpha^{-16/35}f_\mathrm{Edd}^{-6/35}\bigg(\frac{M_\mathrm{BH}}{10^8\hspace{0.5mm}\mathrm{M}_\odot} \bigg)^{4/35},
(in region c), with :math:`R_\mathrm{S}=2R_\mathrm{G}` the Schwarzschild radius. These warp radii are generally of order :math:`\approx1000R_\mathrm{G}`, which can lead to fairly quick alignment of the thin disk with the large-scale angular momentum direction (quicker than any significant evolution in mass or spin magnitude, illustrating why the inclusion of the effects of Lense-Thirring torques is important). (in region c), with :math:`R_\mathrm{S}=2R_\mathrm{G}` the Schwarzschild radius. These warp radii are generally of order :math:`\approx1000R_\mathrm{G}`, which can lead to fairly quick alignment of the thin disk with the large-scale angular momentum direction (quicker than any significant evolution in mass or spin magnitude, illustrating why the inclusion of the effects of Lense-Thirring torques is important).
In the context of thin disks, there is a futher complication. The self-gravity of the disk may become important at large radii (see `Lodato 2007 <https://www.sif.it/riviste/sif/ncr/econtents/2007/030/07/article/0>`_ for a review). The disk will fragment in the region where the Toomre parameter is :math:`Q(R)>1`. We thus assume that the disk extends out to where :math:`Q(R_\mathrm{sg})=1`. The self-gravity radius :math:`R_\mathrm{sg}` can be calculated from this condition and the definition of the Toomre parameter :math:`Q=\Omega c_{\mathrm{s}} /(\pi G \Sigma)`, yielding In the context of thin disks, there is a futher complication. The self-gravity of the disk may become important at large radii (see `Lodato 2007 <https://www.sif.it/riviste/sif/ncr/econtents/2007/030/07/article/0>`_ for a review). The disk will fragment in the region where the Toomre parameter is :math:`Q(R)>1`. We thus assume that the disk extends out to where :math:`Q(R_\mathrm{sg})=1`. The self-gravity radius :math:`R_\mathrm{sg}` can be calculated from this condition and the definition of the Toomre parameter :math:`Q=\Omega c_{\mathrm{s}} /(\pi G \Sigma)`, yielding
.. math:: .. math::
R_{\text {sg,TD,b}}=6460 R_{S} \alpha^{28/51} \dot{m}^{-18/51}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{-49/51} R_{\text {sg,TD,b}}=6460 R_{S} \alpha^{28/51} f_\mathrm{Edd}^{-18/51}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{-49/51}
in region b) and in region b) and
.. math:: .. math::
R_\mathrm{sg,TD,c}=2456 R_{S} \alpha^{28/45} \dot{m}^{-22/45}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{-52/45} R_\mathrm{sg,TD,c}=2456 R_{S} \alpha^{28/45} f_\mathrm{Edd}^{-22/45}\left(\frac{M_{\mathrm{BH}}}{10^{8} M_{\odot}}\right)^{-52/45}
in region c). In all our calculations involving :math:`R_\mathrm{warp}` (for deciding the sign of spin and evolving the direction of angular momentum, as described in the preceeding sections), we always take the minimum of :math:`R_\mathrm{warp}` and :math:`R_\mathrm{sg}`. This is because if :math:`R_\mathrm{sg}<R_\mathrm{warp}`, the entire disk of extent :math:`R_\mathrm{sg}` will be warped. in region c). In all our calculations involving :math:`R_\mathrm{warp}` (for deciding the sign of spin and evolving the direction of angular momentum, as described in the preceeding sections), we always take the minimum of :math:`R_\mathrm{warp}` and :math:`R_\mathrm{sg}`. This is because if :math:`R_\mathrm{sg}<R_\mathrm{warp}`, the entire disk of extent :math:`R_\mathrm{sg}` will be warped.
...@@ -212,7 +253,10 @@ The exact behaviour of the thick and slim disk (which we will collectively call ...@@ -212,7 +253,10 @@ The exact behaviour of the thick and slim disk (which we will collectively call
.. math:: .. math::
R_\mathrm{warp,adv}=R_\mathrm{G}\bigg(\frac{384a}{25(H/R)^2}\bigg)^{2/5}. R_\mathrm{warp,adv}=R_\mathrm{G}\bigg(\frac{384a}{25(H/R)^2}\bigg)^{2/5}.
In our model, we assume that the inner regions of the disks are on average aligned or counteraligned with the spin vector (one can think of this as averaging over the precession, which has periods of :math:`\approx`days, over long enough time scales). For simplicity, we also refer to the radii within which this is true as the warp radii. For both of the advection-dominated disks, these radii are only of order several :math:`R_\mathrm{G}`. Note that similar values are found if one assumes that the Bardeen-Peterson effect operates in these disks. While there are some uncertainties in the assumptions we have made, we point out that using any of these values is much more physically motivated than using thin disk equations (the warp radii of order thousands of :math:`R_\mathrm{G}`), which is what is often done (e.g. `Griffin et al. 2019a <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487..198G/abstract>`_, `Dubois et al. 2012 <https://ui.adsabs.harvard.edu/abs/2014MNRAS.440.1590D/abstract>`_). In our model, we assume that the inner regions of the disks are on
average aligned or counteraligned with the spin vector (one can think
of this as averaging over the precession, which has periods of
:math:`\approx` days, over long enough time scales). For simplicity, we also refer to the radii within which this is true as the warp radii. For both of the advection-dominated disks, these radii are only of order several :math:`R_\mathrm{G}`. Note that similar values are found if one assumes that the Bardeen-Peterson effect operates in these disks. While there are some uncertainties in the assumptions we have made, we point out that using any of these values is much more physically motivated than using thin disk equations (the warp radii of order thousands of :math:`R_\mathrm{G}`), which is what is often done (e.g. `Griffin et al. 2019a <https://ui.adsabs.harvard.edu/abs/2019MNRAS.487..198G/abstract>`_, `Dubois et al. 2012 <https://ui.adsabs.harvard.edu/abs/2014MNRAS.440.1590D/abstract>`_).
In order to determine the sign of spin and evolve the angular momentum direction, expressions for the warp mass :math:`M_\mathrm{warp}` and warp angular momentum :math:`J_\mathrm{warp}` are also needed. We calculate this using surface integrals as In order to determine the sign of spin and evolve the angular momentum direction, expressions for the warp mass :math:`M_\mathrm{warp}` and warp angular momentum :math:`J_\mathrm{warp}` are also needed. We calculate this using surface integrals as
...@@ -234,9 +278,9 @@ where :math:`v_\mathrm{r}=-\alpha v_0 v_\mathrm{K}` is the radial velocity. Here ...@@ -234,9 +278,9 @@ where :math:`v_\mathrm{r}=-\alpha v_0 v_\mathrm{K}` is the radial velocity. Here
Black hole mergers Black hole mergers
------------------ ------------------
In the process of merging, BHs interact in a very complicated manner. Their final spin is not trivial to predict, and it can depend on a very large parameter space (including the mass ratio of the black holes and the relative orientation and magnitude of the spins). Orbital angular momentum plays a role in the merger as well. We use the fitting function found by `Rezzolla et al. (2007) <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.78.044002>`_, whose results have been found to be very accurate in newer and more sophisticated studies that sweep the huge parameter space of possible merger configurations. The only flaw in these formulas is that they do not include the effects of gravitational radiation. However, the effects of this radiation is confined to a :math:`\approx10\%` level, and only if either of the spin vectors is aligned or counteraligned with the direction of the orbital angular momentum (if it is not, the fits are even more accurate). In the process of merging, BHs interact in a very complicated manner. Their final spin is not trivial to predict, and it can depend on a very large parameter space (including the mass ratio of the black holes and the relative orientation and magnitude of the spins). Orbital angular momentum plays a role in the merger as well. We use the fitting function found by `Rezzolla et al. (2009) <https://ui.adsabs.harvard.edu/abs/2009CQGra..26i4023R/abstract>`_, whose results have been found to be very accurate in newer and more sophisticated studies that sweep the huge parameter space of possible merger configurations. These formulas are also applicable to cosmological simulations, since they cover the scenario of inspiral from very large distances.
The final spin, according to `Rezzolla et al. (2007) <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.78.044002>`_ can be calculated as The final spin, according to `Rezzolla et al. (2009) <https://ui.adsabs.harvard.edu/abs/2009CQGra..26i4023R/abstract>`_ can be calculated as
.. math:: .. math::
\mathbf{a}_\mathrm{fin} = \frac{1}{(1+q)^2}(\mathbf{a}_1+\mathbf{a}_2q^2+\mathbf{l}q), \mathbf{a}_\mathrm{fin} = \frac{1}{(1+q)^2}(\mathbf{a}_1+\mathbf{a}_2q^2+\mathbf{l}q),
...@@ -248,4 +292,14 @@ where :math:`q=M_2/M_1` is the mass ratio (such that :math:`M_2<M_1`), :math:`\m ...@@ -248,4 +292,14 @@ where :math:`q=M_2/M_1` is the mass ratio (such that :math:`M_2<M_1`), :math:`\m
\left(\frac{s_{5} \mu+t_{0}+2}{1+q^{2}}\right)\left(\left|\mathbf{a}_{1}\right| \cos \theta+\left|\mathbf{a}_{2}\right| q^{2} \cos \xi\right)+ \\ \left(\frac{s_{5} \mu+t_{0}+2}{1+q^{2}}\right)\left(\left|\mathbf{a}_{1}\right| \cos \theta+\left|\mathbf{a}_{2}\right| q^{2} \cos \xi\right)+ \\
2 \sqrt{3}+t_{2} \mu+t_{3} \mu^{2}. 2 \sqrt{3}+t_{2} \mu+t_{3} \mu^{2}.
Here, :math:`\mu=q/(1+q)^2` is the symmetric mass ratio, and :math:`s_4 = -0.129`, :math:`s_5 = -0.384`, :math:`t_0 = -2.686`, :math:`t_2 = -3.454`, :math:`t_3 = 2.353`. The three cosines depend on the angles between the different vectors which play a role in the merger: :math:`\cos \phi=\hat{\mathbf{a}_{1}} \cdot \hat{\mathbf{a}_{\mathbf{2}}}`, :math:`\cos \theta=\hat{\mathbf{a}_{1}} \cdot \hat{\mathbf{l}}`, :math:`\cos \xi=\hat{\mathbf{a}_{2}} \cdot \hat{\mathbf{l}}`. Here, :math:`\mu=q/(1+q)^2` is the symmetric mass ratio, and :math:`s_4 = -0.1229`, :math:`s_5 = -0.4537`, :math:`t_0 = -2.8904`, :math:`t_2 = -3.5171`, :math:`t_3 = 2.5763`. The three cosines depend on the angles between the different vectors which play a role in the merger: :math:`\cos \phi=\hat{\mathbf{a}_{1}} \cdot \hat{\mathbf{a}_{\mathbf{2}}}`, :math:`\cos \theta=\hat{\mathbf{a}_{1}} \cdot \hat{\mathbf{l}}`, :math:`\cos \xi=\hat{\mathbf{a}_{2}} \cdot \hat{\mathbf{l}}`.
Given the information available within the model, we could in principle calculate the recoil velocity of the remnant, as well as the total mass fraction lost to gravitational waves. We do not implement the former at this stage since we cannot reliably track the movement of black holes in their host galaxies. However, we do implement the latter. We use results from the same series of numerical relativity simulations as above (`Barausse et al. 2012 <https://ui.adsabs.harvard.edu/abs/2012ApJ...758...63B/abstract>`_) and write the final mass of the remnant as:
.. math::
M_\mathrm{BH,fin} = (M_\mathrm{BH,1}+M_\mathrm{BH,2})\Big\{1 - [1 - e_\mathrm{ISCO}(\tilde{a})]\mu - 4\mu^2[4p_0+16p_1\tilde{a}(\tilde{a}+1)+e_\mathrm{ISCO}(\tilde{a})-1]\Big\},
where :math:`p_0=0.04827`, :math:`p_1=0.01707` and :math:`e_\mathrm{ISCO}(\tilde{a})` is the dimensionless specific binding energy at the innermost stable circular orbit calculated using an effective spin variable defined as
.. math::
\tilde{a} = \frac{|\mathbf{a_1}|\cos\theta+|\mathbf{a_2}|\cos\xi}{(1+q)^2}.
.. AGN spin and jet model
Filip Husko, 26 September 2023
.. AGN_spin_jet:
Variable AGN heating temperature model
--------------------------------------
As part of this AGN model we also introduce an option of using variable AGN heating temperatures, as opposed to a constant value. This can be turned on by setting ``AGN_heating_temperature_model`` to ``AGN_heating_temperature_local``. In this scheme, we use three different criteria, all of which are applied simultaneously, and all of which are numerically motivated. These are: 1) the `Dalla Vecchia & Schaye (2012) <https://ui.adsabs.harvard.edu/abs/2012MNRAS.426..140D/abstract>`_ to prevent numerical overcooling, 2) a replenishment condition to make sure that the BH smoothing kernel can be replenished at a rate similar to the rate of evacuation due to heating and 3) a crossing condition that ensures a sufficient separation of heating events in time such that the kernel is never expected to contain more than a single heated particle.
We implement the `Dalla Vecchia & Schaye (2012) <https://ui.adsabs.harvard.edu/abs/2012MNRAS.426..140D/abstract>`_ condition using Eqn. (18) from the paper. We invert it to obtain the minimum heating temperature:
.. math::
\Delta T_\mathrm{DV}= 1.8\times10^6\hspace{0.5mm}\mathrm{K}\hspace{0.5mm}\times\bigg(\frac{N_\mathrm{ngb}m_\mathrm{g}}{60\times10^6\hspace{0.5mm}\mathrm{M}_\odot} \bigg)^{1/3}\bigg(\frac{n_\mathrm{H}}{0.1\hspace{0.5mm}\mathrm{cm}^{-3}} \bigg)^{2/3},
where :math:`N_\mathrm{ngb}` is the number of neighbours in the BH kernel, :math:`m_\mathrm{g}` the average neighbour mass and :math:`n_\mathrm{H}` the number density of gas in the BH kernel.
We derive the second condition by comparing the time-scales of replenishment onto the kernel and evacuation out of it on account of gas heating. The evacuation time-scale is the time required to evacuate the entire kernel of gas, and it is given by :math:`t_\mathrm{evac}=m_\mathrm{g}N_\mathrm{ngb}/\dot{M}`, where :math:`\dot{M}` is the mass flux associated with feedback. It can be expressed using the feedback power :math:`P` and internal energy per unit mass :math:`e` as :math:`\dot{M}=P/e`, where :math:`e=3k_\mathrm{B}\Delta T/2\mu m_\mathrm{p}`, :math:`\Delta T` is the heating temperature, :math:`m_\mathrm{p}` the proton mass and :math:`\mu\approx0.6` the mean molecular weight for ionized gas.
The replenishment time-scale is given by :math:`t_\mathrm{repl}=h/v_\mathrm{repl}`, where :math:`v_\mathrm{repl}` is an effective velocity with which gas can flow inwards to replenish the kernel. Given this formulation, we can set the two time-scales (of evacuation and replenishment) equal to each other and solve for a heating temperature that ensures timely replenishment:
.. math::
\Delta T_\mathrm{repl} = \frac{2\mu m_\mathrm{p}}{3 k_\mathrm{B}}\frac{hP}{N_\mathrm{ngb}m_\mathrm{g}v_\mathrm{repl}}.
Finally, we calculate the replenishment velocity :math:`v_\mathrm{repl}` as follows. We assume that gas can replenish the kernel under the effects of either gas turbulence or pressure. We thus express the replenishment velocity as :math:`v_\mathrm{repl} = \max(\sigma,\Tilde{c}_\mathrm{s})`, where :math:`\sigma` is the velocity dispersion of gas and :math:`\Tilde{c}_\mathrm{s}` an effective sound speed, which we calculate by assuming that there is always some ISM exherting its pressure on gas near the BH, even if all of the gas in the kernel is cold. We choose the form :math:`\Tilde{c}_\mathrm{s}=\max(c_\mathrm{s,hot},\hspace{0.5mm}10\hspace{0.5mm}\mathrm{km}\mathrm{s}^{-1})`, where :math:`c_\mathrm{s,hot}` is the kernel-weighted average sound speed of all particles that have a temperature :math:`T>10^4` K and :math:`10` km :math:`\mathrm{s}^{-1}` is the sound speed of the ISM assuming a typical temperature of :math:`T=10^4` K.
The third and final condition we use is based on the time it takes a single heated particle to cross and exit the kernel before the next one is heated. The time-interval between two heating events is :math:`\Delta = m_\mathrm{g}/\dot{M}`, while the time required for a heated particle to cross the kernel is :math:`\Delta t_\mathrm{cross}= h/c_\mathrm{s,\Delta T}`, where :math:`c_\mathrm{s,\Delta T} = \sqrt{\gamma(\gamma-1)e} = \sqrt{5k_\mathrm{B}\Delta T/3\mu m_\mathrm{p}}` is the sound speed of the heated gas. Equating these two time-scales, we obtain the final heating temperature:
.. math::
\Delta T_\mathrm{cross} = \frac{\mu m_\mathrm{p}}{k_\mathrm{B}}\bigg(\frac{2hP}{\sqrt{15}m_\mathrm{g}}\bigg)^{2/3}.
The final heating temperature scheme we use can be written as:
.. math::
\Delta T = \max[\xi \max(\Delta T_\mathrm{DV}, \Delta T_\mathrm{repl}, \Delta T_\mathrm{cross}), \Delta T_\mathrm{min}],
where :math:`\Delta T_\mathrm{min}` is an additional temperature floor meant to prevent very-low temperature heating events, and :math:`\xi` is an additional free parameter that can be used to rescale heating temperatures from all three conditions to higher values, which may be necessary to correctly calibrate the simulations. Specifically, one may use the hot gas fractions to choose a value of :math:`\xi` if :math:`\Delta T_\mathrm{min}` is set to a low value of :math:`\approx10^{7.5}` K, or one may set :math:`\xi=1`, thus using the numerical conditions exactly as derived, and instead calibrate the simulations by varying :math:`\Delta T_\mathrm{min}`.
...@@ -36,7 +36,7 @@ Star formation ...@@ -36,7 +36,7 @@ Star formation
The AGORA model uses the :ref:`GEAR model <gear_star_formation>` scheme, however with the The AGORA model uses the :ref:`GEAR model <gear_star_formation>` scheme, however with the
``GEARStarFormation:star_formation_mode`` parameter set to ``agora``. Instead of requiring the gas ``GEARStarFormation:star_formation_mode`` parameter set to ``agora``. Instead of requiring the gas
density to reach the pressure floor, we simply require it to be denser than a density density to reach the pressure floor, we simply require it to be denser than a density
threshold defined by ``GEARStarFormation:density_threshold``. threshold defined by ``GEARStarFormation:density_threshold_Hpcm3``.
Recommended parameters for the AGORA model should be: Recommended parameters for the AGORA model should be:
...@@ -47,10 +47,10 @@ Recommended parameters for the AGORA model should be: ...@@ -47,10 +47,10 @@ Recommended parameters for the AGORA model should be:
GEARStarFormation: GEARStarFormation:
star_formation_mode: agora star_formation_mode: agora
star_formation_efficiency: 0.01 star_formation_efficiency: 0.01
maximal_temperature: 1e10 maximal_temperature_K: 1e10
n_stars_per_particle: 1 n_stars_per_particle: 1
min_mass_frac: 0.5 min_mass_frac: 0.5
density_threshold: 1.67e-23 density_threshold_Hpcm3: 10
...@@ -127,9 +127,6 @@ Recommended parameters for the AGORA model should be: ...@@ -127,9 +127,6 @@ Recommended parameters for the AGORA model should be:
.. _agora_pressure_floor: .. _agora_pressure_floor:
Pressure Floor Pressure Floor
...@@ -138,4 +135,12 @@ Pressure Floor ...@@ -138,4 +135,12 @@ Pressure Floor
The AGORA model uses precisely the same pressure floor than the :ref:`GEAR model <gear_pressure_floor>`. The AGORA model uses precisely the same pressure floor than the :ref:`GEAR model <gear_pressure_floor>`.
.. _agora_initial_conditions:
Initial Conditions
~~~~~~~~~~~~~~~~~~
Note that if in the initial conditions, the time of formation of a stellar particle is given (``BirthTime``)
and set to a negative value, the stellar particle will provide no feedback.
A similar behavior will be obtained if the parameter ``overwrite_birth_time`` is set to 1 and
``birth_time`` to -1.
...@@ -5,6 +5,25 @@ ...@@ -5,6 +5,25 @@
Basic model (others) Basic model (others)
==================== ====================
Sinks: Simple Bondi-Hoyle accretion
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The ``Basic`` sink model provides a foundation on which new sink implementations could be built. It includes a prescription for Bondi-Hoyle gas accretion, and a method for sink-sink mergers that is a slightly simplified version of the implementation used in GEAR.
No other physics is implemented for this model. Sinks cannot form - to use this model, sink particles must already be present in the initial conditions. They also cannot spawn stars from the gas they accrete.
Bondi-Hoyle accretion can be done in one of two ways:
* Gas particles within the sink's kernel are stochastically swallowed entirely, with a probability set by the Bondi-Hoyle rate. Specifically, the probability is set by the current difference between the sink's subgrid mass (determined by the accretion rate) and its dynamical mass (which tracks the number of particles/sinks actually swallowed). This mode is equivalent to the EAGLE black hole accretion model.
* Gas particles within the sink's kernel are "nibbled" down to some minimal mass, which can be specified by the user. This method is equivalent to the black hole accretion model of Bahe et al. 2022.
This model has only two parameters that must be specified in your parameter ``yml`` file:
* ``BasicSink:use_nibbling``: determines whether accretion is done by "nibbling" or by swallowing outright.
* ``BasicSink:min_gas_mass_for_nibbling_Msun``: if using "nibbling", the minimum mass to which gas particles can be nibbled. A good default is half the original particle mass.
For an even more bare-bones starting point, the ``Default`` sink model contains no physics at all, and is a totally blank canvas on which to build your sink model.
Cooling: Analytic models Cooling: Analytic models
~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~
......
...@@ -537,6 +537,15 @@ Note that the star formation rates are expressed in internal units and not in ...@@ -537,6 +537,15 @@ Note that the star formation rates are expressed in internal units and not in
solar masses per year as is the case in many other codes. This choice ensures solar masses per year as is the case in many other codes. This choice ensures
consistency between all the fields written to the snapshots. consistency between all the fields written to the snapshots.
Finally, the star formation model can also create more than one star particle
per gas particle in a star formation event. The number of particles generated is
controlled by a runtime parameter. All the stars generated share the same
properties. They are slightly displaced from each others using a random vector
of magnitude :math:`0.1h` added to the position of the gas particle. If only
one star is formed (as is the default), no displacement is added. If N stars are
formed per gas particle, each star is born with a mass of 1/N of the gas
particle's mass.
For a normal EAGLE run, that section of the parameter file reads: For a normal EAGLE run, that section of the parameter file reads:
.. code:: YAML .. code:: YAML
...@@ -559,7 +568,8 @@ For a normal EAGLE run, that section of the parameter file reads: ...@@ -559,7 +568,8 @@ For a normal EAGLE run, that section of the parameter file reads:
threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3. threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
min_over_density: 57.7 # Over-density above which star-formation is allowed. min_over_density: 57.7 # Over-density above which star-formation is allowed.
EOS_entropy_margin_dex: 0.5 # (Optional) Logarithm base 10 of the maximal entropy above the EOS at which stars can form. EOS_entropy_margin_dex: 0.5 # (Optional) Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
num_of_stars_per_gas_particle: 1 # (Optional) The number star particles to form per gas particle converted to stars. (Defaults to 1. Must be > 0)
Alternatively, the code can also use a simple Schmidt law for the SF rate Alternatively, the code can also use a simple Schmidt law for the SF rate
:math:`\dot{m}_* = \epsilon_{ff} \times m_g \times \frac{3 \pi} {32 G :math:`\dot{m}_* = \epsilon_{ff} \times m_g \times \frac{3 \pi} {32 G
......
.. GEAR sub-grid model chemistry
Darwin Roduit, 30th March 2025
.. gear_chemistry:
.. _gear_chemistry:
Chemistry
=========
Feedback mechanisms such as supernova feedback transfer metals to the nearby gas particles. There is no mass exchange (mass advection) in SPH methods, as in grid-based or moving-mesh hydrodynamics solvers. As a result, the metals received by the gas particles are locked in these particles. Grid or moving-mesh codes have mass advection and often implement passive scalar advection of metals to model metal mixing. GEAR implements different methods to model metal mixing.
.. _gear_smoothed_metallicity:
Smoothed metallicity
--------------------
The smoothed metallicity scheme consists in using the SPH to smooth the metallicity of each particle over the neighbors. It is worth to point the fact that we are *not exchanging* any metals but only smoothing it. The parameter ``GEARChemistry:initial_metallicity`` set the (non smoothed) initial mass fraction of each element for all the particles and ``GEARChemistry:scale_initial_metallicity`` use the feedback table to scale the initial metallicity of each element according the Sun's composition. If ``GEARChemistry:initial_metallicity`` is negative, then the metallicities are read from the initial conditions.
For this chemistry scheme the parameters are:
.. code:: YAML
GEARChemistry:
initial_metallicity: 1 # Initial metallicity of the gas (mass fraction)
scale_initial_metallicity: 1 # Should we scale the initial metallicity with the solar one?
.. GEAR sub-grid model feedback and stellar evolution
Loic Hausammann, 17th April 2020
Darwin Roduit, 30th March 2025
.. _gear_stellar_evolution_and_feedback:
Stellar evolution and feedback
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The feedback is composed of a few different models:
- The initial mass function (IMF) defines the quantity of each type of stars,
- The lifetime of a star defines when a star will explode (or simply die),
- The supernovae of type II (SNII) defines the rates and yields,
- The supernovae of type Ia (SNIa) defines the rates and yields,
- The energy injection that defines how to inject the energy / metals into the particles.
Most of the parameters are defined inside a table (``GEARFeedback:yields_table``) but can be override with some parameters in the YAML file.
I will not describe theses parameters more than providing them at the end of this section.
Two different models exist for the supernovae (``GEARFeedback:discrete_yields``).
In the continuous mode, we integrate the quantities over the IMF and then explodes a floating point number of stars (can be below 1 in some cases).
In the discrete mode, we avoid the problem of floating points by rounding the number of supernovae (using a floor and randomly adding a supernovae depending on the fractional part) and then compute the properties for a single star at a time.
Initial mass function
^^^^^^^^^^^^^^^^^^^^^
GEAR is using the IMF model from `Kroupa (2001) <https://ui.adsabs.harvard.edu/abs/2001MNRAS.322..231K/abstract>`_.
We have a difference of 1 in the exponent due to the usage of IMF in mass and not in number.
We also restrict the mass of the stars to be inside :math:`[0.05, 50] M_\odot`.
Here is the default model used, but it can be easily adapted through the initial mass function parameters:
.. math::
\xi(m) \propto m^{-\alpha_i}\, \textrm{where}\,
\begin{cases}
\alpha_0 = 0.3,\, & 0.01 \leq m / M_\odot < 0.08, \\
\alpha_1 = 1.3,\, & 0.08 \leq m / M_\odot < 0.50, \\
\alpha_2 = 2.3,\, & 0.50 \leq m / M_\odot < 1.00, \\
\alpha_3 = 2.3,\, & 1.00 \leq m / M_\odot,
\end{cases}
Lifetime
^^^^^^^^
The lifetime of a star in GEAR depends only on two parameters: first its mass and then its metallicity.
.. math::
\log(\tau(m)) = a(Z) \log^2(m) + b(Z) \log(m) + c(Z) \\ \\
a(Z) = -40.110 Z^2 + 5.509 Z + 0.7824 \\
b(Z) = 141.929 Z^2 - 15.889 Z - 3.2557 \\
c(Z) = -261.365 Z^2 + 17.073 Z + 9.8661
where :math:`\tau` is the lifetime in years, :math:`m` is the mass of the star (in solar mass) and Z the metallicity of the star.
The parameters previously given are the default ones, they can be modified in the parameters file.
Supernovae II
^^^^^^^^^^^^^
The supernovae rate is simply given by the number of stars massive enough that end their life at the required time.
.. math::
\dot{N}_\textrm{SNII}(t) = \int_{M_l}^{M_u} \delta(t - \tau(m)) \frac{\phi(m)}{m} \mathrm{d}m
where :math:`M_l` and :math:`M_u` are the lower and upper mass limits for a star exploding in SNII, :math:`\delta` is the Dirac function and :math:`\phi` is the initial mass function (in mass).
The yields for SNII cannot be written in an analytical form, they depend on a few different tables that are based on the work of `Kobayashi et al. (2000) <https://ui.adsabs.harvard.edu/abs/2000ApJ...539...26K/abstract>`_ and `Tsujimoto et al. (1995) <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..945T/abstract>`_.
Supernovae Ia
^^^^^^^^^^^^^
The supernovae Ia are a bit more complicated as they involve two different stars.
.. math::
\dot{N}_\textrm{SNIa}(t) = \left( \int_{M_{p,l}}^{M_{p,u}} \frac{\phi(m)}{m} \mathrm{d}m \right) \sum_i b_i \int_{M_{d,l,i}}^{M_{d,u,i}}
\delta(t-\tau(m)) \frac{\phi_d(m)}{m}\mathrm{d}m
.. math::
\phi_d(m) \propto m^{-0.35}
where :math:`M_{p,l}` and :math:`M_{p,u}` are the mass limits for a progenitor of a white dwarf, :math:`b_i` is the probability to have a companion and
:math:`M_{d,l,i}` and :math:`M_{d,u,i}` are the mass limits for each type of companion.
The first parenthesis represents the number of white dwarfs and the second one the probability to form a binary.
+------------------+--------------------+-------------------+------------------+
| Companion | :math:`M_{d,l,i}` | :math:`M_{d,u,i}` | :math:`b_i` |
+==================+====================+===================+==================+
| Red giant | 0.9 | 1.5 | 0.02 |
+------------------+--------------------+-------------------+------------------+
| Main sequence | 1.8 | 2.5 | 0.05 |
+------------------+--------------------+-------------------+------------------+
The yields are based on the same papers than the SNII.
Supernova energy injection
^^^^^^^^^^^^^^^^^^^^^^^^^^
When a star goes into a supernova (type II and Ia), the stellar evolution determines how much energy (``GEARFeedback:supernovae_energy_erg`` TO BE CHECKED), mass and metals are released during the explosion. The energy can be ditributed as internal/thermal energy or as momentum. Thus, we need to distribute internal energy, momentum, mass and metals to the gas particles. We will group all these in the “fluxes” term.
We have two models for the distribution of these fluxes and the subgrid modelling of the supernovae : GEAR model and GEAR mechanical model. We describe the two schemes in the :ref:`gear_sn_feedback_models` page.
Generating a new table
^^^^^^^^^^^^^^^^^^^^^^
The feedback table is an HDF5 file with the following structure:
.. graphviz:: feedback_table.dot
where the solid (dashed) squares represent a group (a dataset) with the name of the object underlined and the attributes written below. Everything is in solar mass or without units (e.g. mass fraction or unitless constant).
In ``Data``, the attribute ``elts`` is an array of string with the element names (the last should be ``Metals``, it corresponds to the sum of all the elements), ``MeanWDMass`` is the mass of the white dwarfs
and ``SolarMassAbundances`` is an array of float containing the mass fraction of the different element in the sun.
In ``IMF``, ``n + 1`` is the number of part in the IMF, ``as`` are the exponent (``n+1`` elements), ``ms`` are the mass limits between each part (``n`` elements) and
``Mmin`` (``Mmax``) is the minimal (maximal) mass of a star.
In ``LifeTimes``, the coefficient are given in the form of a single table (``coeff_z`` with a 3x3 shape).
In ``SNIa``, ``a`` is the exponent of the distribution of binaries, ``bb1`` and ``bb2`` are the coefficient :math:`b_i` and the other attributes follow the same names than in the SNIa formulas.
The ``Metals`` group from the ``SNIa`` contains the name of each elements (``elts``) and the metal mass fraction ejected by each supernovae (``data``) in the same order. They must contain the same elements than in ``Data``.
Finally for the ``SNII``, the mass limits are given by ``Mmin`` and ``Mmax``. For the yields, the datasets required are ``Ej`` (mass fraction ejected [processed]), ``Ejnp`` (mass fraction ejected [non processed]) and one dataset for each element present in ``elts``. The datasets should all have the same size, be uniformly sampled in log and contains the attributes ``min`` (mass in log for the first element) and ``step`` (difference of mass in log between two elements).
.. code:: YAML
GEARFeedback:
supernovae_energy_erg: 0.1e51 # Energy released by a single supernovae.
yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 # Table containing the yields.
discrete_yields: 0 # Should we use discrete yields or the IMF integrated one?
GEARInitialMassFunction:
number_function_part: 4 # Number of different part in the IMF
exponents: [0.7, -0.8, -1.7, -1.3] # Exponents of each part of the IMF
mass_limits_msun: [0.05, 0.08, 0.5, 1, 50] # Limits in mass between each part of the IMF
GEARLifetime:
quadratic: [-40.1107, 5.50992, 0.782432] # Quadratic terms in the fit
linear: [141.93, -15.8895, -3.25578] # Linear terms in the fit
constant: [-261.366, 17.0735, 9.86606] # Constant terms in the fit
GEARSupernovaeIa:
exponent: -0.35 # Exponent for the distribution of companions
min_mass_white_dwarf_progenitor: 3 # Minimal mass of a progenitor of white dwarf
max_mass_white_dwarf_progenitor: 8 # Maximal mass of a progenitor of white dwarf
max_mass_red_giant: 1.5 # Maximal mass for a red giant
min_mass_red_giant: 0.9 # Minimal mass for a red giant
coef_red_giant: 0.02 # Coefficient for the distribution of red giants companions
max_mass_main_sequence: 2.6 # Maximal mass for a main sequence star
min_mass_main_sequence: 1.8 # Minimal mass for a main sequence star
coef_main_sequence: 0.05 # Coefficient for the distribution of main sequence companions
white_dwarf_mass: 1.38 # Mass of a white dwarf
GEARSupernovaeII:
interpolation_size: 200 # Number of elements for the interpolation of the data
.. GEAR sub-grid model
Loic Hausammann, 17th April 2020
Darwin Roduit, 30th March 2025
.. _gear_pressure_floor:
Pressure Floor
~~~~~~~~~~~~~~
In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles.
This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by:
.. math::
P_\textrm{Jeans} = \frac{\rho}{\gamma} \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3}
where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
:math:`h` the kernel support and :math:`N_\textrm{Jeans}` (``GEARPressureFloor:jeans_factor`` in the parameter file) is the number of particle required in order to resolve a clump.
This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2, SPHENIX and Pressure-Energy) have the floor available.
In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
.. math::
m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
and simply replace the :math:`P_i, P_j` by the pressure with the floor (when the pressure is below the floor).
Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
.. code:: YAML
GEARPressureFloor:
jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
.. _gear_star_formation:
Star formation
~~~~~~~~~~~~~~
The star formation is done in two steps: first we check if a particle is in the star forming regime and then we use a stochastic approach to transform the gas particles into stars.
A particle is in the star forming regime if:
- The velocity divergence is negative (:math:`\nabla\cdot v < 0`),
- The temperature is lower than a threshold (:math:`T < T_t` where :math:`T_t` is defined with ``GEARStarFormation:maximal_temperature_K``),
- The gas density is higher than a threshold (:math:`\rho > \rho_t` where :math:`\rho_t` is defined with ``GEARStarFormation:density_threshold_Hpcm3``)
- The particle reaches the pressure floor (:math:`\rho > \frac{\pi}{4 G N_\textrm{Jeans}^{2/3} h^2}\frac{\gamma k_B T}{\mu m_p}` where :math:`N_\textrm{Jeans}` is defined in the pressure floor).
If ``GEARStarFormation:star_formation_mode`` is set to ``agora``, the condition on the pressure floor is ignored. Its default value is ``default``.
A star will be able to form if a randomly drawn number is below :math:`\frac{m_g}{m_\star}\left(1 - \exp\left(-c_\star \Delta t / t_\textrm{ff}\right)\right)` where :math:`t_\textrm{ff}` is the free fall time, :math:`\Delta t` is the time step of the particle and :math:`c_\star` is the star formation coefficient (``GEARStarFormation:star_formation_efficiency``), :math:`m_g` the mass of the gas particle and :math:`m_\star` the mass of the possible future star. The mass of the star is computed from the average gas mass in the initial conditions divided by the number of possible stars formed per gas particle (``GEARStarFormation:n_stars_per_particle``). When we cannot have enough mass to form a second star (defined with the fraction of mass ``GEARStarFormation:min_mass_frac``), we fully convert the gas particle into a stellar particle. Once the star is formed, we move it a bit in a random direction and fraction of the smoothing length in order to avoid any division by 0.
Currently, only the following hydro schemes are compatible: SPHENIX, Gadget2, minimal SPH, Gasoline-2 and Pressure-Energy.
Implementing the other hydro schemes is not complicated but requires some careful thinking about the cosmological terms in the definition of the velocity divergence (comoving vs non comoving coordinates and if the Hubble flow is included or not).
.. code:: YAML
GEARStarFormation:
star_formation_efficiency: 0.01 # star formation efficiency (c_*)
maximal_temperature_K: 3e4 # Upper limit to the temperature of a star forming particle
density_threshold_Hpcm3: 10 # Density threshold (Hydrogen atoms/cm^3) for star formation
n_stars_per_particle: 4 # Number of stars that an hydro particle can generate
min_mass_frac: 0.5 # Minimal mass for a stellar particle as a fraction of the average mass for the stellar particles.
Initial Conditions
++++++++++++++++++
Note that if in the initial conditions, the time of formation of a stellar particle is given (``BirthTime``)
and set to a negative value, the stellar particle will provide no feedback.
A similar behavior will be obtained if the parameter ``Stars:overwrite_birth_time`` is set to 1 and
``Stars:birth_time`` to -1.
.. _gear_grackle_cooling:
Cooling: Grackle
~~~~~~~~~~~~~~~~
Grackle is a chemistry and cooling library presented in `B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_
(do not forget to cite if used). Four different modes are available:
equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
H\\(_2^+\\)) and 12 species (adds D, D\\(^+\\) and HD). Following the same
order, the swift cooling options are ``grackle_0``, ``grackle_1``, ``grackle_2``
and ``grackle_3`` (the numbers correspond to the value of
``primordial_chemistry`` in Grackle). It also includes some metal cooling (on/off with ``GrackleCooling:with_metal_cooling``), self-shielding
methods and UV background (on/off with ``GrackleCooling:with_UV_background``). In order to use the Grackle cooling, you will need
to provide a HDF5 table computed by Cloudy (``GrackleCooling:cloudy_table``).
Configuring and compiling SWIFT with Grackle
++++++++++++++++++++++++++++++++++++++++++++
In order to compile SWIFT with Grackle, you need to provide the options ``with-chemistry=GEAR`` and ``with-grackle=$GRACKLE_ROOT``
where ``$GRACKLE_ROOT`` is the root of the install directory (not the ``lib``).
.. warning::
(State 2023) Grackle is experiencing current development, and the API is subject
to changes in the future. For convenience, a frozen version is hosted as a fork
on github here: https://github.com/mladenivkovic/grackle-swift .
The version available there will be tried and tested and ensured to work with
SWIFT.
Additionally, that repository hosts files necessary to install that specific
version of grackle with spack.
To compile it, run
the following commands from the root directory of Grackle:
``./configure; cd src/clib``.
Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
the file ``src/clib/Make.mach.linux-gnu``.
Finish with ``make machine-linux-gnu; make && make install``.
Note that we require the 64 bit float version of Grackle, which should be the default setting.
(The precision can be set while compiling grackle with ``make precision-64``).
If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
Parameters
++++++++++
When starting a simulation without providing the different element fractions in the non equilibrium mode, the code supposes an equilibrium and computes them automatically.
The code uses an iterative method in order to find the correct initial composition and this method can be tuned with two parameters. ``GrackleCooling:max_steps`` defines the maximal number of steps to reach the convergence and ``GrackleCooling:convergence_limit`` defines the tolerance in the relative error.
In the parameters file, a few different parameters are available.
- ``GrackleCooling:redshift`` defines the redshift to use for the UV background (for cosmological simulation, it must be set to -1 in order to use the simulation's redshift).
- ``GrackleCooling:provide_*_heating_rates`` can enable the computation of user provided heating rates (such as with the radiative transfer) in either volumetric or specific units.
- Feedback can be made more efficient by turning off the cooling during a few Myr (``GrackleCooling:thermal_time_myr``) for the particles touched by a supernovae.
- The self shielding method is defined by ``GrackleCooling:self_shielding_method`` where 0 means no self shielding, > 0 means a method defined in Grackle (see Grackle documentation for more information) and -1 means GEAR's self shielding that simply turn off the UV background when reaching a given density (``GrackleCooling:self_shielding_threshold_atom_per_cm3``).
- The initial elemental abundances can be specified with ``initial_nX_to_nY_ratio``, e.g. if you want to specify an initial HII to H abundance. A negative value ignores the parameter. The complete list can be found below.
A maximal (physical) density must be set with the ``GrackleCooling:maximal_density_Hpcm3 parameter``. The density passed to Grackle is *the minimum of this density and the gas particle (physical) density*. A negative value (:math:`< 0`) deactivates the maximal density, i.e. there is no maximal density limit.
The purpose of this parameter is the following. The Cloudy tables provided by Grackle are limited in density (typically to :math:`10^4 \; \mathrm{hydrogen \; atoms/cm}^3`). In high-resolution simulations, particles can have densities higher than :math:`10^4 \; \mathrm{hydrogen \; atoms/cm}^3`. This maximal density ensures that we pass a density within the interpolation ranges of the table, should the density exceed it.
It can be a solution to some of the following errors (with a translation of what the values mean):
.. code:: text
inside if statement solve rate cool: 0 0
MULTI_COOL iter > 10000 at j,k = 1 1
FATAL error (2) in MULTI_COOL
dt = 1.092E-08 ttmin = 7.493E-12
2.8E-19 // sub-cycling timestep
7.5E-12 // time elapsed (in the sub-cycle)
2.2E+25 // derivative of the internal energy
T
.. note::
This problem is particularly relevant with metal cooling enabled. Another solution is to modify the tables. But one is not exempted from exceeding the table maximal density value, since Grackle does not check if the particle density is larger than the table maximal density.
Complete parameter list
-----------------------
Here is the complete section in the parameter file:
.. code:: YAML
GrackleCooling:
cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
with_UV_background: 1 # Enable or not the UV background
redshift: 0 # Redshift to use (-1 means time based redshift)
with_metal_cooling: 1 # Enable or not the metal cooling
provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
max_steps: 10000 # (optional) Max number of step when computing the initial composition
convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
maximal_density_Hpcm3: 1e4 # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
HydrogenFractionByMass : 1. # Hydrogen fraction by mass (default is 0.76)
use_radiative_transfer : 1 # Arrays of ionization and heating rates are provided
RT_heating_rate_cgs : 0 # heating rate in units of / nHI_cgs
RT_HI_ionization_rate_cgs : 0 # HI ionization rate in cgs [1/s]
RT_HeI_ionization_rate_cgs : 0 # HeI ionization rate in cgs [1/s]
RT_HeII_ionization_rate_cgs: 0 # HeII ionization rate in cgs [1/s]
RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s]
volumetric_heating_rates_cgs: 0 # Volumetric heating rate in cgs [erg/s/cm3]
specific_heating_rates_cgs: 0 # Specific heating rate in cgs [erg/s/g]
H2_three_body_rate : 1 # Specific the H2 formation three body rate (0->5,see Grackle documentation)
H2_cie_cooling : 0 # Enable/disable H2 collision-induced emission cooling from Ripamonti & Abel (2004)
H2_on_dust: 0 # Flag to enable H2 formation on dust grains
local_dust_to_gas_ratio : -1 # The ratio of total dust mass to gas mass in the local Universe (-1 to use the Grackle default value).
cmb_temperature_floor : 1 # Enable/disable an effective CMB temperature floor
initial_nHII_to_nH_ratio: -1 # initial nHII to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHeI_to_nH_ratio: -1 # initial nHeI to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHeII_to_nH_ratio: -1 # initial nHeII to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHeIII_to_nH_ratio: -1 # initial nHeIII to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nDI_to_nH_ratio: -1 # initial nDI to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nDII_to_nH_ratio: -1 # initial nDII to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHM_to_nH_ratio: -1 # initial nHM to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nH2I_to_nH_ratio: -1 # initial nH2I to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nH2II_to_nH_ratio: -1 # initial nH2II to nH ratio (number density ratio). Value is ignored if set to -1.
initial_nHDI_to_nH_ratio: -1 # initial nHDI to nH ratio (number density ratio). Value is ignored if set to -1.
.. note::
A simple example running SWIFT with Grackle can be find in ``examples/Cooling/CoolingBox``. A more advanced example combining heating and cooling (with heating and ionization sources) is given in ``examples/Cooling/CoolingHeatingBox``. ``examples/Cooling/CoolingWithPrimordialElements/`` runs a uniform cosmological box with imposed abundances and let them evolve down to redshift 0.
...@@ -6,269 +6,16 @@ GEAR model ...@@ -6,269 +6,16 @@ GEAR model
=========== ===========
GEAR's model are mainly described in `Revaz \& Jablonka <https://ui.adsabs.harvard.edu/abs/2018A%26A...616A..96R/abstract>`_. GEAR's model are mainly described in `Revaz \& Jablonka <https://ui.adsabs.harvard.edu/abs/2018A%26A...616A..96R/abstract>`_.
This model can be selected with the configuration option ``--with-subgrid=GEAR`` and run with the option ``--gear``. A few examples exist and can be found in ``examples/GEAR``. This model can be selected with the configuration option ``--with-subgrid=GEAR`` and run with the option ``--gear``. A few examples exist and can be found in ``examples/GEAR``, ``examples/SinkParticles`` or ``IsolatedGalaxy/IsolatedGalaxy_multi_component``.
.. _gear_pressure_floor:
Pressure Floor .. toctree::
~~~~~~~~~~~~~~ :caption: Table of Contents
In order to avoid the artificial collapse of unresolved clumps, a minimum in pressure is applied to the particles. gear_model
This additional pressure can be seen as the pressure due to unresolved hydrodynamics turbulence and is given by: chemistry
feedback
supernova_feedback
sinks/index
output
.. math::
P_\textrm{Jeans} = \frac{\rho}{\gamma} \frac{4}{\pi} G h^2 \rho N_\textrm{Jeans}^{2/3}
where :math:`\rho` is the density, :math:`\gamma` the adiabatic index, :math:`G` is the gravitational constant,
:math:`h` the kernel support and :math:`N_\textrm{Jeans}` (``GEARPressureFloor:jeans_factor`` in the parameter file) is the number of particle required in order to resolve a clump.
This must be directly implemented into the hydro schemes, therefore only a subset of schemes (Gadget-2, SPHENIX and Pressure-Energy) have the floor available.
In order to implement it, you need equation 12 in `Hopkins 2013 <https://arxiv.org/abs/1206.5006>`_:
.. math::
m_i \frac{\mathrm{d}v_i}{\mathrm{d}t} = - \sum_j x_i x_j \left[ \frac{P_i}{y_i^2} f_{ij} \nabla_i W_{ij}(h_i) + \frac{P_j}{y_j^2} f_{ji} \nabla_j W_{ji}(h_j) \right]
and simply replace the :math:`P_i, P_j` by the pressure with the floor (when the pressure is below the floor).
Here the :math:`x, y` are simple weights that should never have the pressure floor included even if they are related to the pressure (e.g. pressure-entropy).
.. code:: YAML
GEARPressureFloor:
jeans_factor: 10. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
.. _gear_grackle_cooling:
Cooling: Grackle
~~~~~~~~~~~~~~~~
Grackle is a chemistry and cooling library presented in `B. Smith et al. 2017 <https://ui.adsabs.harvard.edu/abs/2017MNRAS.466.2217S>`_
(do not forget to cite if used). Four different modes are available:
equilibrium, 6 species network (H, H\\( ^+ \\), e\\( ^- \\), He, He\\( ^+ \\)
and He\\( ^{++} \\)), 9 species network (adds H\\(^-\\), H\\(_2\\) and
H\\(_2^+\\)) and 12 species (adds D, D\\(^+\\) and HD). Following the same
order, the swift cooling options are ``grackle_0``, ``grackle_1``, ``grackle_2``
and ``grackle_3`` (the numbers correspond to the value of
``primordial_chemistry`` in Grackle). It also includes some metal cooling (on/off with ``GrackleCooling:with_metal_cooling``), self-shielding
methods and UV background (on/off with ``GrackleCooling:with_UV_background``). In order to use the Grackle cooling, you will need
to provide a HDF5 table computed by Cloudy (``GrackleCooling:cloudy_table``).
When starting a simulation without providing the different element fractions in the non equilibrium mode, the code supposes an equilibrium and computes them automatically.
The code uses an iterative method in order to find the correct initial composition and this method can be tuned with two parameters. ``GrackleCooling:max_steps`` defines the maximal number of steps to reach the convergence and ``GrackleCooling:convergence_limit`` defines the tolerance in the relative error.
In order to compile SWIFT with Grackle, you need to provide the options ``with-chemistry=GEAR`` and ``with-grackle=$GRACKLE_ROOT``
where ``$GRACKLE_ROOT`` is the root of the install directory (not the ``lib``).
You will need a Grackle version later than 3.1.1. To compile it, run
the following commands from the root directory of Grackle:
``./configure; cd src/clib``.
Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
the file ``src/clib/Make.mach.linux-gnu``.
Finish with ``make machine-linux-gnu; make && make install``.
Note that we require the 64 bit float version of Grackle, which should be the default setting.
(The precision can be set while compiling grackle with ``make precision-64``).
If you encounter any problem, you can look at the `Grackle documentation <https://grackle.readthedocs.io/en/latest/>`_
You can now provide the path given for ``MACH_INSTALL_PREFIX`` to ``with-grackle``.
In the parameters file, a few different parameters are available.
``GrackleCooling:redshift`` defines the redshift to use for the UV background (for cosmological simulation, it must be set to -1 in order to use the simulation's redshift) and ``GrackleCooling:provide_*_heating_rates`` can enable the computation of user provided heating rates (such as with the radiative transfer) in either volumetric or specific units.
For the feedback, it can be made more efficient by turning off the cooling during a few Myr (``GrackleCooling:thermal_time_myr``) for the particles touched by a supernovae.
The self shielding method is defined by ``GrackleCooling:self_shielding_method`` where 0 means no self shielding, > 0 means a method defined in Grackle (see Grackle documentation for more information) and -1 means GEAR's self shielding that simply turn off the UV background when reaching a given density (``GrackleCooling:self_shielding_threshold_atom_per_cm3``).
.. code:: YAML
GrackleCooling:
cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
with_UV_background: 1 # Enable or not the UV background
redshift: 0 # Redshift to use (-1 means time based redshift)
with_metal_cooling: 1 # Enable or not the metal cooling
provide_volumetric_heating_rates: 0 # (optional) User provide volumetric heating rates
provide_specific_heating_rates: 0 # (optional) User provide specific heating rates
max_steps: 10000 # (optional) Max number of step when computing the initial composition
convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
.. _gear_star_formation:
Star formation
~~~~~~~~~~~~~~
The star formation is done in two steps: first we check if a particle is in the star forming regime and then we use a stochastic approach to transform the gas particles into stars.
A particle is in the star forming regime if:
- The velocity divergence is negative (:math:`\nabla\cdot v < 0`),
- The temperature is lower than a threshold (:math:`T < T_t` where :math:`T_t` is defined with ``GEARStarFormation:maximal_temperature``),
- The gas density is higher than a threshold (:math:`\rho > \rho_t` where :math:`\rho_t` is defined with ``GEARStarFormation:density_threshold``)
- The particle reaches the pressure floor (:math:`\rho > \frac{\pi}{4 G N_\textrm{Jeans}^{2/3} h^2}\frac{\gamma k_B T}{\mu m_p}` where :math:`N_\textrm{Jeans}` is defined in the pressure floor).
If ``GEARStarFormation:star_formation_mode`` is set to ``agora``, the condition on the pressure floor is ignored. Its default value is ``default``.
A star will be able to form if a randomly drawn number is below :math:`\frac{m_g}{m_\star}\left(1 - \exp\left(-c_\star \Delta t / t_\textrm{ff}\right)\right)` where :math:`t_\textrm{ff}` is the free fall time, :math:`\Delta t` is the time step of the particle and :math:`c_\star` is the star formation coefficient (``GEARStarFormation:star_formation_efficiency``), :math:`m_g` the mass of the gas particle and :math:`m_\star` the mass of the possible future star. The mass of the star is computed from the average gas mass in the initial conditions divided by the number of possible stars formed per gas particle (``GEARStarFormation:n_stars_per_particle``). When we cannot have enough mass to form a second star (defined with the fraction of mass ``GEARStarFormation:min_mass_frac``), we fully convert the gas particle into a stellar particle. Once the star is formed, we move it a bit in a random direction and fraction of the smoothing length in order to avoid any division by 0.
Currently, only the following hydro schemes are compatible: SPHENIX and Gadget2.
Implementing the other hydro schemes is not complicated but requires some careful thinking about the cosmological terms in the definition of the velocity divergence (comoving vs non comoving coordinates and if the Hubble flow is included or not).
.. code:: YAML
GEARStarFormation:
star_formation_efficiency: 0.01 # star formation efficiency (c_*)
maximal_temperature: 3e4 # Upper limit to the temperature of a star forming particle
n_stars_per_particle: 4 # Number of stars that an hydro particle can generate
min_mass_frac: 0.5 # Minimal mass for a stellar particle as a fraction of the average mass for the stellar particles.
Chemistry
~~~~~~~~~
In the chemistry, we are using the smoothed metallicity scheme that consists in using the SPH to smooth the metallicity of each particle over the neighbors. It is worth to point the fact that we are not exchanging any metals but only smoothing it. The parameter ``GEARChemistry:initial_metallicity`` set the (non smoothed) initial mass fraction of each element for all the particles and ``GEARChemistry:scale_initial_metallicity`` use the feedback table to scale the initial metallicity of each element according the Sun's composition.
.. code:: YAML
GEARChemistry:
initial_metallicity: 1 # Initial metallicity of the gas (mass fraction)
scale_initial_metallicity: 1 # Should we scale the initial metallicity with the solar one?
Feedback
~~~~~~~~
The feedback is composed of a few different models:
- The initial mass function (IMF) defines the quantity of each type of stars,
- The lifetime of a star defines when a star will explode (or simply die),
- The supernovae of type II (SNII) defines the rates and yields,
- The supernovae of type Ia (SNIa) defines the rates and yields,
- The energy injection that defines how to inject the energy / metals into the particles.
Most of the parameters are defined inside a table (``GEARFeedback:yields_table``) but can be override with some parameters in the YAML file.
I will not describe theses parameters more than providing them at the end of this section.
Two different models exist for the supernovae (``GEARFeedback:discrete_yields``).
In the continuous mode, we integrate the quantities over the IMF and then explodes a floating point number of stars (can be below 1 in some cases).
In the discrete mode, we avoid the problem of floating points by rounding the number of supernovae (using a floor and randomly adding a supernovae depending on the fractional part) and then compute the properties for a single star at a time.
Initial mass function
^^^^^^^^^^^^^^^^^^^^^
GEAR is using the IMF model from `Kroupa (2001) <https://ui.adsabs.harvard.edu/abs/2001MNRAS.322..231K/abstract>`_.
We have a difference of 1 in the exponent due to the usage of IMF in mass and not in number.
We also restrict the mass of the stars to be inside :math:`[0.05, 50] M_\odot`.
Here is the default model used, but it can be easily adapted through the initial mass function parameters:
.. math::
\xi(m) \propto m^{-\alpha_i}\, \textrm{where}\,
\begin{cases}
\alpha_0 = 0.3,\, & 0.01 \leq m / M_\odot < 0.08, \\
\alpha_1 = 1.3,\, & 0.08 \leq m / M_\odot < 0.50, \\
\alpha_2 = 2.3,\, & 0.50 \leq m / M_\odot < 1.00, \\
\alpha_3 = 2.3,\, & 1.00 \leq m / M_\odot,
\end{cases}
Lifetime
^^^^^^^^
The lifetime of a star in GEAR depends only on two parameters: first its mass and then its metallicity.
.. math::
\log(\tau(m)) = a(Z) \log^2(m) + b(Z) \log(m) + c(Z) \\ \\
a(Z) = -40.110 Z^2 + 5.509 Z + 0.7824 \\
b(Z) = 141.929 Z^2 - 15.889 Z - 3.2557 \\
c(Z) = -261.365 Z^2 + 17.073 Z + 9.8661
where :math:`\tau` is the lifetime in years, :math:`m` is the mass of the star (in solar mass) and Z the metallicity of the star.
The parameters previously given are the default ones, they can be modified in the parameters file.
Supernovae II
^^^^^^^^^^^^^
The supernovae rate is simply given by the number of stars massive enough that end their life at the required time.
.. math::
\dot{N}_\textrm{SNII}(t) = \int_{M_l}^{M_u} \delta(t - \tau(m)) \frac{\phi(m)}{m} \mathrm{d}m
where :math:`M_l` and :math:`M_u` are the lower and upper mass limits for a star exploding in SNII, :math:`\delta` is the Dirac function and :math:`\phi` is the initial mass function (in mass).
The yields for SNII cannot be written in an analytical form, they depend on a few different tables that are based on the work of `Kobayashi et al. (2000) <https://ui.adsabs.harvard.edu/abs/2000ApJ...539...26K/abstract>`_ and `Tsujimoto et al. (1995) <https://ui.adsabs.harvard.edu/abs/1995MNRAS.277..945T/abstract>`_.
Supernovae Ia
^^^^^^^^^^^^^
The supernovae Ia are a bit more complicated as they involve two different stars.
.. math::
\dot{N}_\textrm{SNIa}(t) = \left( \int_{M_{p,l}}^{M_{p,u}} \frac{\phi(m)}{m} \mathrm{d}m \right) \sum_i b_i \int_{M_{d,l,i}}^{M_{d,u,i}}
\delta(t-\tau(m)) \frac{\phi_d(m)}{m}\mathrm{d}m
.. math::
\phi_d(m) \propto m^{-0.35}
where :math:`M_{p,l}` and :math:`M_{p,u}` are the mass limits for a progenitor of a white dwarf, :math:`b_i` is the probability to have a companion and
:math:`M_{d,l,i}` and :math:`M_{d,u,i}` are the mass limits for each type of companion.
The first parenthesis represents the number of white dwarfs and the second one the probability to form a binary.
+------------------+--------------------+-------------------+------------------+
| Companion | :math:`M_{d,l,i}` | :math:`M_{d,u,i}` | :math:`b_i` |
+==================+====================+===================+==================+
| Red giant | 0.9 | 1.5 | 0.02 |
+------------------+--------------------+-------------------+------------------+
| Main sequence | 1.8 | 2.5 | 0.05 |
+------------------+--------------------+-------------------+------------------+
The yields are based on the same papers than the SNII.
Energy injection
^^^^^^^^^^^^^^^^
All the supernovae (type II and Ia) inject the same amount of energy into the surrounding gas (``GEARFeedback:supernovae_energy_erg``) and distribute it according to the hydro kernel.
The same is done with the metals and the mass.
Generating a new table
^^^^^^^^^^^^^^^^^^^^^^
The feedback table is an HDF5 file with the following structure:
.. graphviz:: feedback_table.dot
where the solid (dashed) squares represent a group (a dataset) with the name of the object underlined and the attributes written below. Everything is in solar mass or without units (e.g. mass fraction or unitless constant).
In ``Data``, the attribute ``elts`` is an array of string with the element names (the last should be ``Metals``, it corresponds to the sum of all the elements), ``MeanWDMass`` is the mass of the white dwarfs
and ``SolarMassAbundances`` is an array of float containing the mass fraction of the different element in the sun.
In ``IMF``, ``n + 1`` is the number of part in the IMF, ``as`` are the exponent (``n+1`` elements), ``ms`` are the mass limits between each part (``n`` elements) and
``Mmin`` (``Mmax``) is the minimal (maximal) mass of a star.
In ``LifeTimes``, the coefficient are given in the form of a single table (``coeff_z`` with a 3x3 shape).
In ``SNIa``, ``a`` is the exponent of the distribution of binaries, ``bb1`` and ``bb2`` are the coefficient :math:`b_i` and the other attributes follow the same names than in the SNIa formulas.
The ``Metals`` group from the ``SNIa`` contains the name of each elements (``elts``) and the metal mass fraction ejected by each supernovae (``data``) in the same order. They must contain the same elements than in ``Data``.
Finally for the ``SNII``, the mass limits are given by ``Mmin`` and ``Mmax``. For the yields, the datasets required are ``Ej`` (mass fraction ejected [processed]), ``Ejnp`` (mass fraction ejected [non processed]) and one dataset for each element present in ``elts``. The datasets should all have the same size, be uniformly sampled in log and contains the attributes ``min`` (mass in log for the first element) and ``step`` (difference of mass in log between two elements).
.. code:: YAML
GEARFeedback:
supernovae_energy_erg: 0.1e51 # Energy released by a single supernovae.
yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5 # Table containing the yields.
discrete_yields: 0 # Should we use discrete yields or the IMF integrated one?
GEARInitialMassFunction:
number_function_part: 4 # Number of different part in the IMF
exponents: [0.7, -0.8, -1.7, -1.3] # Exponents of each part of the IMF
mass_limits_msun: [0.05, 0.08, 0.5, 1, 50] # Limits in mass between each part of the IMF
GEARLifetime:
quadratic: [-40.1107, 5.50992, 0.782432] # Quadratic terms in the fit
linear: [141.93, -15.8895, -3.25578] # Linear terms in the fit
constant: [-261.366, 17.0735, 9.86606] # Constant terms in the fit
GEARSupernovaeIa:
exponent: -0.35 # Exponent for the distribution of companions
min_mass_white_dwarf_progenitor: 3 # Minimal mass of a progenitor of white dwarf
max_mass_white_dwarf_progenitor: 8 # Maximal mass of a progenitor of white dwarf
max_mass_red_giant: 1.5 # Maximal mass for a red giant
min_mass_red_giant: 0.9 # Minimal mass for a red giant
coef_red_giant: 0.02 # Coefficient for the distribution of red giants companions
max_mass_main_sequence: 2.6 # Maximal mass for a main sequence star
min_mass_main_sequence: 1.8 # Minimal mass for a main sequence star
coef_main_sequence: 0.05 # Coefficient for the distribution of main sequence companions
white_dwarf_mass: 1.38 # Mass of a white dwarf
GEARSupernovaeII:
interpolation_size: 200 # Number of elements for the interpolation of the data
.. Sink particles in GEAR model
Darwin Roduit, 14 July 2024
.. sink_GEAR_model:
Snapshots ouputs
----------------
Here, we provide a summary of the quantities written in the snapshots, in addition to positions, velocities, masses, smoothing lengths and particle IDs.
Sink particles
~~~~~~~~~~~~~~
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| Name | Description | Units | Comments |
+=======================================+=============================================+========================+===================================================+
| ``NumberOfSinkSwallows`` | | Number of merger events with other sinks | [-] | |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``NumberOfGasSwallows`` | | Number of gas particles accreted | [-] | |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``TargetMass`` | | Target mass required to spawn the next | [U_M] | | You can use it to determine if the target mass |
| | | star particle | | | is so huge that the sink's mass cannot spawn |
| | | | | such a star. Such rare behaviour may bias the |
| | | | | IMF towards high masses. |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``Nstars`` | | Number of stars created by this sink | [-] | |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``SwallowedAngularMomentum`` | | Total angular momentum of accreted | [U_M U_L^2 U_T^{-1}] | |
| | | material | | |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``MetalMassFractions`` | | Mass fraction of each tracked metal | [-] | | *Only in GEAR chemistry module.* |
| | | element | | | Array of length ``N`` (number of elements), |
| | | | | set at compile time via |
| | | | | ``--with-chemistry=GEAR_N``. |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``BirthScaleFactors`` | | Scale factor at the time of sink creation | [-] | | Only used in *cosmological* runs. |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``BirthTimes`` | | Time when the sink was created | [U_T] | | Only used in *non-cosmological* runs. |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
Stars
~~~~~
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| Name | Description | Units | Comments |
+=======================================+=============================================+========================+===================================================+
| ``BirthScaleFactors`` | | Scale-factors when the stars were born | [-] | | Only used in cosmological runs. |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``BirthTimes`` | | Time when the stars were born | [U_T] | | Only used in non-cosmological runs. |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``BirthMasses`` | | Masses of the stars at birth time | [U_M] | | SF and sinks modules |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``ProgenitorIDs`` | | ID of the progenitor sinks or gas | [-] | | SF and sinks modules |
| | | particles | | |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``BirthDensities`` | | Gas density at star formation | [U_M U_L^{-3}] | | *Only in SF module* |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``BirthTemperatures`` | | Gas temperature at star formation | [K] | | *Only in SF module* |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``Potentials`` | | Gravitational potential of the star | [U_L^2 U_T^{-2}] | |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``StellarParticleType`` | | Type of the stellar particle: | [-] | | 0: (discrete) single star |
| | | | | | 1: continuous IMF part star |
| | | | | | 2: single population star |
| | | | | | The last type corresponds to legacy IMF stars. |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
| ``MetalMassFractions`` | | Mass fraction of each metal element | [-] | | *Only in GEAR chemistry module*. |
| | | | | | Array of length ``N`` (number of elements), |
| | | | | | set at compile time by |
| | | | | | ``--with-chemistry=GEAR_N``. |
+---------------------------------------+---------------------------------------------+------------------------+---------------------------------------------------+
Gas particles
~~~~~~~~~~~~~
Since hydro scheme writes its own set of outputs, we only provide the outputs that ``GEAR`` writes for gas particles.
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| Name | Description | Units | Comments |
+==========================================+=============================================================+========================+====================================================+
| ``SmoothedMetalMassFractions`` | | Mass fraction of each metal element | [-] | | *Only in GEAR chemistry module.* |
| | | smoothed over the SPH kernel | | | Array of length ``N``, set at compile time by |
| | | | | ``--with-chemistry=GEAR_N`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``MetalMassFractions`` | | Raw (non-smoothed) mass fraction of | [-] | | *Only in GEAR chemistry module.* |
| | | each metal element | | | Same layout as above. |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``HI`` | | Mass fraction of neutral H (:math:`\mathrm{H}`) | [-] | | *Only if* ``GRACKLE_1 to 3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``HII`` | | Mass fraction of ionized H (:math:`\mathrm{H}^+`) | [-] | | *Only if* ``GRACKLE_1 to 3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``HeI`` | | Mass fraction of neutral He (:math:`\mathrm{He}`) | [-] | | *Only if* ``GRACKLE_1 to 3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``HeII`` | | Mass fraction of singly ionized He (:math:`\mathrm{He}^+`)| [-] | | *Only if* ``GRACKLE_1 to 3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``HeIII`` | | Mass fraction of doubly ionized He | [-] | | *Only if* ``GRACKLE_1 to 3`` |
| | | (:math:`\mathrm{He}^{++}`) | | |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``e`` | | Free electron mass fraction (:math:`\mathrm{e}^-`) | [-] | | *Only if* ``GRACKLE_1 to 3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``HM`` | | Mass fraction of :math:`\mathrm{H}^-` | [-] | | *Only if* ``GRACKLE_2 or 3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``H2I`` | | Mass fraction of neutral :math:`\mathrm{H}_2` | [-] | | *Only if* ``GRACKLE_2 or 3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``H2II`` | | Mass fraction of ionized :math:`\mathrm{H}_2^+` | [-] | | *Only if* ``GRACKLE_2 or 3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``DI`` | | Mass fraction of neutral D (:math:`\mathrm{D}`) | [-] | | *Only if* ``GRACKLE_3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``DII`` | | Mass fraction of ionized D (:math:`\mathrm{D}^+`) | [-] | | *Only if* ``GRACKLE_3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
| ``HDI`` | | Mass fraction of :math:`\mathrm{HD}` | [-] | | *Only if* ``GRACKLE_3`` |
+------------------------------------------+-------------------------------------------------------------+------------------------+----------------------------------------------------+
.. Sink particles in GEAR model
Darwin Roduit, 17 April 2024
.. sink_GEAR_model:
.. _sink_GEAR_model:
Sink particles and star formation in GEAR model
===============================================
.. toctree::
:caption: Table of Contents
introduction
theory
sink_timesteps
params
.. Sink particles in GEAR model
Darwin Roduit, 15 March 2024
.. sink_GEAR_model:
Introduction
------------
GEAR sink particles provide an alternative model to star formation. Instead of stochastically transforming gas particles into stars as is done in GEAR star formation scheme under some conditions, we transform a gas into a sink particle. The main property of the sink particle is its accretion radius. When gas particles within this accretion radius are eligible to be swallowed by the sink, we remove them and transfer their mass, momentum, angular momentum, chemistry properties, etc to the sink particle.
With the sink particles, the IMF splits into two parts: the continuous part and the discrete part. Those parts will correspond to two kinds of stars. Particles in the discrete part of the IMF represent individual stars. It means that discrete IMF-sampled stars have different masses. Particles in the continuous part represent a population of stars, all with the same mass.
The sink particle will randomly choose a target mass, accrete gas until it reaches this target mass and finally spawn a star. Then, the sink chooses a new target mass and repeats the same procedure. When stars are spawned, they are given a new position and velocity as well as chemical properties.
In ``theory.rst`` we outline all of the theory which is implemented as part of the model. This includes how sink particles are formed, how the gas is accreted, how sinks are merged and how stars are spawned. In ``params.rst`` we list and discuss all parameters used by the model. Below we outline how to configure and run the model.
Compiling and running the model
-------------------------------
You can configure the model with ``--with-sink=GEAR`` in combination with other configure options of the GEAR model. The model will then be used when the ``--sinks`` flag is among the runtime options. In particular, the sink particles require ``--feedback`` runtime option.
Notice that you also need to compile with ``--with-star-formation=GEAR``. The star formation module is required to collect and write star formation data. However, you do not need to run Swift with the option ``--star-formation``.
Then, you do not need to do anything special. Sink particles will be created during your runs. If you want, you can have sink particles in your ICs. At the moment, sink particles do not have any special fields to be set.
A full list of all relevant parameters of the model is in :ref:`sink_GEAR_parameters`. We also briefly describe the most important parameters which need to be set to run the model, as well as how to run it in different configurations.
.. warning::
Currently, MPI is not implemented for the sink particles. If you try to run with MPI, you will encounter an error. We thus recommend configuring Swift with ``--disable-mpi`` to avoid any surprises.
.. Sink particles in GEAR model
Darwin Roduit, 15 March 2024
.. sink_GEAR_model:
.. _sink_GEAR_parameters:
Model parameters
----------------
The parameters of the GEAR sink model are grouped into the ``GEARSink`` section of the parameter file.
The first three parameters are:
* Whether we are using a fixed cut-off radius for gas and sink accretion: ``use_fixed_cut_off_radius``,
* The sink cut-off radius: ``cut_off_radius``,
* The sink inner accretion radius fraction in terms of the cut-off radius: ``f_acc``.
The ``use_fixed_cut_off_radius`` is mandatory and should be set to 0 or 1. If set to 1, the GEAR model will use a fixed cutoff radius equal to the value of ``cut_off_radius``. If not, the cutoff radius is allowed to vary according to the local gas density. In the code, the cutoff radius is always equal to the sink's smoothing length multiplied by a constant factor ``kernel_gamma`` - setting a fixed cutoff will fix the smoothing length at the appropriate value for all sinks.
The ``cut_off_radius`` parameter is optional, and is ignored if ``use_fixed_cut_off_radius`` is 0. If a fixed cut-off is used, every sink's smoothing length will be permanently set to this number divided by ``kernel_gamma`` on formation.
The ``f_acc`` parameter is also optional. Its default value is :math:`0.1`. Its value must respect :math:`0 \leq f_\text{acc} \leq 1` . It describes the inner radius :math:`f_{\text{acc}} \cdot r_{\text{cut-off}}` in which gas particles are swallowed without any further checks, as explained below.
The next three mandatory parameters are:
* the gas temperature threshold to form a sink when :math:`\rho_\text{threshold} < \rho_\text{gas} < \rho_\text{maximal}` :``temperature_threshold_K``,
* the minimal gas threshold density required to form a sink:``density_threshold_Hpcm3``,
* the maximal density at which the temperature check is not performed:``maximal_density_threshold_Hpcm3`` (Default: ``FLT_MAX``).
These three parameters govern the first two criteria of the sink formation scheme. If these criteria are not passed, sink particles are not created. If they are passed, the code performs further checks to form sink particles. Some of those criteria checks can be disabled, as explained below.
The next set of parameters deals with the sampling of the IMF and the representation of star particles:
* minimal mass of stars represented by discrete particles: ``minimal_discrete_mass_Msun``,
* mass of the stellar particle representing the continuous part of the IMF: ``stellar_particle_mass_Msun``,
* minimal mass of the first stars represented by discrete particles: ``minimal_discrete_mass_first_stars_Msun``,
* mass of the stellar particle (first stars) representing the continuous part of the IMF:: ``stellar_particle_mass_first_stars_Msun``.
With sink particles, star particles can represent either a single star or a population of stars in the low mass part of the IMF (continuous IMF sampling). The stars in the continuous part of the IMF are put together in a particle of mass ``stellar_particle_mass_Msun`` or ``stellar_particle_mass_first_stars_Msun``, while individual stars in the discrete part have their mass sampled from the IMF. The limit between the continuous and discrete sampling of the IMF is controlled by ``minimal_discrete_mass_Msun`` and ``minimal_discrete_mass_first_stars_Msun``.
The next set of parameters controls the sink formation scheme. More details are provided in the GEAR documentation. Here is a brief overview:
* whether or not the gas must be contracting: ``sink_formation_contracting_gas_criterion`` (default: 1),
* whether or not the gas smoothing length must be small enough: ``sink_formation_smoothing_length_criterion`` (default: 1),
* whether or not the gas must be in a Jeans unstable state: ``sink_formation_jeans_instability_criterion`` (default: 1),
* whether or not the gas must be in a bound state: ``sink_formation_bound_state_criterion`` (default: 1),
* whether or not a new sink can be formed in a region where its ``cut_off_radius`` and the one of an existing sink overlap: ``sink_formation_overlapping_sink_criterion`` (default: 1).
Those criteria are checked if the density and temperature criteria are successfully passed. They control the behaviour of the sink formation scheme. By default, they are all activated and set to ``1``.
The next parameter is ``disable_sink_formation`` (default: 0). It controls whether sinks are formed or not in the simulation. The main purpose is when sinks are put in initial conditions and sinks are not wanted to be added during the run. This parameter is set to ``0`` by default, i.e. sink formation is *enabled*.
The next set of parameters deals with the sink time-steps:
* Courant-Friedrich-Levy constant for the CFL-like time-step constraint:``CFL_condition``,
* age (in Myr) at which a sink is considered dead (no accretion) and without time-step limitations, except for 2-body encounters involving another young/old sink and gravity: ``timestep_age_threshold_unlimited_Myr`` (default: FLT_MAX),
* age (in Myr) at which sinks switch from young to old for time-stepping purposes:``timestep_age_threshold`` (default: FLT_MAX),
* maximal time-step length of young sinks (in Myr): ``max_timestep_young_Myr`` (default: FLT_MAX),
* maximal time-step length of old sinks (in Myr): ``max_timestep_old_Myr`` (default: FLT_MAX),
* number of times the IMF mass can be swallowed in a single time-step: ``n_IMF`` (default: FLT_MAX).
* tolerance parameter for SF timestep constraint: ``tolerance_SF_timestep`` (default: 0.5)
The last parameter is ``sink_minimal_mass_Msun``. This parameter is mainly intended for low-resolution simulations with :math:`m_\text{gas} > 100 \; M_\odot`. It prevents :math:`m_\text{sink} \ll m_\text{gas}` simulations when sinks spawn stars, which can lead to gravity run away.
The full section is:
.. code:: YAML
GEARSink:
use_fixed_cut_off_radius: 1 # Are we using a fixed cutoff radius? If we are, in GEAR the cutoff radius is fixed at the value specified below, and the sink smoothing length is fixed at this value divided by kernel_gamma. If not, the cutoff radius varies with the sink smoothing length as r_cut = h*kernel_gamma.
cut_off_radius: 1e-3 # Cut off radius of all the sinks in internal units. Ignored if use_fixed_cut_off_radius is 0.
f_acc: 0.1 # (Optional) Fraction of the cut_off_radius that determines if a gas particle should be swallowed wihtout additional check. It has to respect 0 <= f_acc <= 1. (Default: 0.1)
temperature_threshold_K: 100 # Max temperature (in K) for forming a sink when density_threshold_Hpcm3 <= density <= maximal_density_threshold_Hpcm3.
density_threshold_Hpcm3: 1e3 # Minimum gas density (in Hydrogen atoms/cm3) required to form a sink particle.
maximal_density_threshold_Hpcm3: 1e5 # If the gas density exceeds this value (in Hydrogen atoms/cm3), a sink forms regardless of temperature if all other criteria are passed. (Default: FLT_MAX)
sink_minimal_mass_Msun: 0. # (Optional) Sink minimal mass in Msun. This parameter prevents m_sink << m_gas in low resolution simulations. (Default: 0.0)
stellar_particle_mass_Msun: 20 # Mass of the stellar particle representing the low mass stars (continuous IMF sampling) (in solar mass)
minimal_discrete_mass_Msun: 8 # Minimal mass of stars represented by discrete particles (in solar mass)
stellar_particle_mass_first_stars_Msun: 20 # Mass of the stellar particle representing the low mass stars (continuous IMF sampling) (in solar mass). First stars
minimal_discrete_mass_first_stars_Msun: 8 # Minimal mass of stars represented by discrete particles (in solar mass). First stars
star_spawning_sigma_factor: 0.2 # Factor to rescale the velocity dispersion of the stars when they are spawned. (Default: 0.2)
sink_formation_contracting_gas_criterion: 1 # (Optional) Activate the contracting gas check for sink formation. (Default: 1)
sink_formation_smoothing_length_criterion: 1 # (Optional) Activate the smoothing length check for sink formation. (Default: 1)
sink_formation_jeans_instability_criterion: 1 # (Optional) Activate the two Jeans instability checks for sink formation. (Default: 1)
sink_formation_bound_state_criterion: 1 # (Optional) Activate the bound state check for sink formation. (Default: 1)
sink_formation_overlapping_sink_criterion: 1 # (Optional) Activate the overlapping sink check for sink formation. (Default: 1)
disable_sink_formation: 0 # (Optional) Disable sink formation. (Default: 0)
CFL_condition: 0.5 # Courant-Friedrich-Levy condition for time integration.
timestep_age_threshold_unlimited_Myr: 100. # (Optional) Age above which sinks no longer have time-step restrictions, except for 2-body encounters involving another young/old sink and gravity (in Mega-years). (Default: FLT_MAX)
timestep_age_threshold_Myr: 25. # (Optional) Age at which sink switch from young to old for time-stepping purposes (in Mega-years). (Default: FLT_MAX)
max_timestep_young_Myr: 1.0 # (Optional) Maximal time-step length of young sinks (in Mega-years). (Default: FLT_MAX)
max_timestep_old_Myr: 5.0 # (Optional) Maximal time-step length of old sinks (in Mega-years). (Default: FLT_MAX)
n_IMF: 2 # (Optional) Number of times the IMF mass can be swallowed in a single timestep. (Default: FLTM_MAX)
tolerance_SF_timestep: 0.5 # (Optional) Tolerance parameter for SF timestep constraint. (Default: 0.5)
.. warning::
Some parameter choices can greatly impact the outcome of your simulations. Think twice when choosing them.
Sink accretion radius
~~~~~~~~~~~~~~~~~~~~~
The most critical parameter is ``cut_off_radius``. As explained in the theory, to form a sink, the gas smoothing kernel edge :math:`\gamma_k h` (:math:`\gamma_k` is a kernel dependent constant) must be smaller than ``cut_off_radius`` (if this criterion is enabled). Therefore, the cut-off radius strongly depends on the resolution of your simulations. Moreover, if you use a minimal gas smoothing length `h`, and plan to use sink particles, consider whether the cut-off radius will meet the smoothing length criterion. If `h` never meets the aforementioned criterion, you will never form sinks and thus never have stars.
On the contrary, if you set a too high cut-off radius, then sinks will accrete a lot of gas particles and spawn a lot of stars in the same cell, which the code might not like and crash with the error:
``runner_others.c:runner_do_star_formation_sink():274: Too many stars in the cell tree leaf! The sorting task will not be able to perform its duties. Possible solutions: (1) The code need to be run with different star formation parameters to reduce the number of star particles created. OR (2) The size of the sorting stack must be increased in runner_sort.c.``
This problem can be mitigated by choosing a higher value of ``stellar_particle_mass_Msun`` and ``stellar_particle_mass_first_stars_Msun``, or higher values of ``minimal_discrete_mass_Msun`` and ``minimal_discrete_mass_first_stars_Msun``. Of course, this comes at the price of having fewer individual stars. Finally, all parameters will depend on your needs.
*If you do not want to change your parameters*, you can increase the ``sort_stack_size`` variable at the beginning ``runner_sort.c``. The default value is 10 in powers of 2 (so the stack size is 1024 particles). Increase it to the desired value. Be careful to not overestimate this.
Note that if a cutoff radius is not specified, and the radius is instead left to vary with the local gas density, the smoothing length criterion is always satisfied.
Guide to choose the the accretion radius or the density threshold
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
We provide some advice to help you set up the sink accretion radius or the threshold density appropriately.
First, you must choose either the sink accretion radius or the threshold density. Choosing the density might be easier based on your previous work or if you have an expected star formation density. Once you fix the density or the accretion radius, you can use the following formula to estimate the remaining parameter. In the code, the gas smoothing length is determined with:
.. math::
h = \eta \left( \frac{X_{\text{H}} m_B}{m_{\text{H}} n_{\text{H}}} \right)^{1/3} \, ,
where :math:`\eta` is a constant related to the number of neighbours in the kernel, :math:`X_{\text{H}}` is the hydrogen mass fraction, :math:`m_B` the gas particle's mass, :math:`m_{\text{H}}` the hydrogen particle mass and :math:`n_{\text{H}}` the hydrogen number density.
Let us provide an example. In GEAR, we do not model physical processes below the parsec scale. Hence, let us take :math:`h \sim 1` pc. In zoom-in simulations we have :math:`m_B \simeq 95 \; M_{\odot}`. The remaining parameters are :math:`\eta = 1.2348` and :math:`X_{\text{H}} = 0.76`. So, after inverting the formula, we find :math:`n_H \simeq 5500 \text{ hydrogen atoms/cm}^3`. In practice, we use :math:`n_H = 1000 \text{ hydrogen atoms/cm}^3`, close to the estimation, and an accretion radius :math:`r_{\text{acc}} = 10` pc. These values are slightly different for safety reasons, but they are consistent.
Remember that this was a way, among others, to determine good accretion radius and threshold density. It can help you with your first runs with sink particles.
Comment on star formation efficiency
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Notice that this model does not have parameters to control the star formation rate of the sink. The SFR is self-regulated by the gas/sink accretion and other feedback mechanisms. Supernovae tend to create bubbles of lower density at the site of star formation, removing the gas and preventing further gas accretion. However, the sink might run into this stack size problem by the time the first supernovae explode. Other pre-stellar feedback mechanisms could do the job earlier, though they are not implemented in GEAR.
.. note::
We provide a piece of general advice: do some calibration on low-resolution simulations. This will help to see what works and what does not work. Keep in mind that you might want to put a higher ``stellar_particle_mass_X_Msun`` at the beginning to avoid spawning too many stars. For the high-resolution simulations, you then can lower the particle's mass.
doc/RTD/source/SubgridModels/GEAR/sinks/sink_accretion_radius.png

506 KiB

doc/RTD/source/SubgridModels/GEAR/sinks/sink_imf.png

90.5 KiB

doc/RTD/source/SubgridModels/GEAR/sinks/sink_overlapping.png

665 KiB

doc/RTD/source/SubgridModels/GEAR/sinks/sink_scheme.png

372 KiB

.. Sink particles in GEAR model
Darwin Roduit, 24 November 2024
.. _sink_GEAR_timesteps:
Sink timesteps
~~~~~~~~~~~~~~
Sink particles interact with the surrounding gas through accretion. To accurately follow the local gas dynamics and assess the stability of gas hydrodynamics, we must impose timestep constraints on the sink particles.
First, we want to ensure sinks know the *local dynamics* and thus they obey a Courant-Friedrichs-Lewy (CFL)-like timestep constraint:
.. math::
\Delta t_s \leq \Delta t_\text{CFL} = C_\text{CFL} \frac{r_{\text{cut, min}}}{\sqrt{c_{s,s}^2 + \| \Delta \mathbf{v}_s \|^2}} \; ,
where :math:`r_\text{cut, min} = \min(r_{\text{cut}, s}, \gamma_k \min_j(h_j))` is the minimal cut-off radius between the sink :math:`s` and the gas neighbours:math:`j`, :math:`c_{s, s}` and :math:`\Delta \mathbf{v}_s` are the gas sound-speed and the relative gas-sink velocity at the sink location. The latter two are reconstructed at the sink location with SPH interpolation. The value of :math:`C_\text{CFL}` is given in the YAML parameter file with ``GEARSink:CFL_condition``.
Since sink particles accrete gas, they must anticipate the surrounding gas infall. We achieve this with a gas free-fall time criterion similar to `Grudic et al. (2021) <https://academic.oup.com/mnras/article/506/2/2199/6276745>`_:
.. math::
\Delta t_s \leq \Delta t_\text{ff} = \sqrt{ \frac{3 \pi}{32 G \rho_s} } \quad \text{with} \quad \rho_s = \frac{3 m_s}{4 \pi {r_{\text{cut, min}}}} \; ,
with :math:`m_s` the sink mass.
These constraints ensure smooth gas accretion by removing a few particles per timestep. This is important since, in SPH, gas particles serve as interpolation points. By preventing the removal of a large number of particles, we ensure the stability of the hydrodynamic computations.
Sink 2-body encounters
++++++++++++++++++++++
To accurately follow *sink mergers*, we implemented `Grudic et al. (2021) <https://academic.oup.com/mnras/article/506/2/2199/6276745>`_ two body timestep constraints between all sink particles:
.. math::
\Delta t_s \leq \Delta t_\text{2-body} = \frac{ t_\text{c, min} t_\text{dyn, min}}{t_\text{c, min} + t_\text{dyn, min}} \; ,
with
.. math::
\quad t_\text{c, min} = \min_{s \neq s'} \frac{ |\varphi^{-1}(r_{ss'}, \, H_\text{sink})| }{v_{ss'}} \quad \text{and} \quad t_\text{dyn, min} = \min_{s \neq s'} \sqrt{ \frac{ |\varphi'(r_{ss'}, \, H_\text{sink})|^{-1}} { G (m_s + m_{s'})} } \; ,
where :math:`r_{ss'}` is the sinks relative separation, :math:`v_{ss'}` the relative velocity, :math:`m_{s}` and :math:`m_{s'}` their masses and :math:`H_\text{sink}` is the sink (fixed) gravitational softening. The function :math:`\varphi(r, H)` is the potential corresponding to the Wendland C2 kernel density field (see `Schaller et al. (2024) <https://doi.org/10.1093/mnras/stae922>`_ section 4.1) and :math:`\varphi'(r, H) \equiv \frac{\mathrm{d} \varphi(r, H)}{\mathrm{d} r}` its derivative.
Timesteps per sink's age categories
+++++++++++++++++++++++++++++++++++
We also implemented maximal timesteps sizes depending on the sink age; :math:`\Delta t_\text{max,s}^\text{age}`. A sink can be young, old or dead. In the first two cases, the sink's timestep is :math:`\min(\Delta t_\text{max,s}^\text{age}, \Delta t_s)`. In the last case, we impose :math:`\Delta t_\text{2-body}` only if a dead sink is involved in a two-boy encounter with an alive sink. Otherwise, the sink has no timestep constraint (apart from gravity). The parameters controlling the transition between the young and old is ``GEARSink:timestep_age_threshold``, and the one between old and dead is ``GEARSink:timestep_age_threshold_unlimited_Myr``. The maximal timesteps are given by ``GEARSink:max_timestep_young_Myr`` and ``GEARSink:max_timestep_old_Myr``.
Notice that sink particles also satisfy a gravity timestep constraint, as do all gravitational particles in Swift.
Star formation constraint
+++++++++++++++++++++++++
Although the stars' masses are sampled from the IMF, the stars' metallicities are not. If a sink accretes mass, it can create many star particles simultaneously. However, these stars will all have the same metallicity, which does not represent the actual metals' evolution during the accretion. In such a situation, the galaxies' properties are affected and do not represent the underlying physics.
Another problem is that we can spawn many stars simultaneously, and the code may complain. Such a situation could be better. Although the constraints in the previous section will help, more is needed. Our solution is to introduce a new accretion criterion using the IMF properties. However, since our politics is that accretion should be feedback-regulated and not based on an arbitrary accretion rate, we reduce the sink time step to avoid limiting the star formation rate to an arbitrary value.
The new accretion criterion is the following. The swallowed gas and sink mass does not exceed ``n_IMF`` times the IMF mass (see the IMF sampling section), but make sure to swallow at least one particle: :math:`M_\text{swallowed} \leq n_\text{IMF} M_\text{IMF} \text{ or } M_\text{swallowed} = 0`.
Since we artificially restrict mass accretion, we keep track of the mass :math:`M_\text{eligible}` that would be swallowed without this criterion. Then, we compute the error :math:`\Delta M` between the restricted and unrestricted swallow. The absolute error is :math:`\Delta M = M_\text{swallowed} - M_\text{eligible}` and the relative error is :math:`| \Delta M | / M_\text{eligible}`.
When :math:`\Delta M < 0` (i.e. :math:`M_\text{swallowed} \neq M_\text{eligible}`), we know the accretion was restricted and we can apply another time-step contraint. To compute a timestep, we convert :math:`\Delta M` to accretion rate by dividing by :math:`\Delta t_\text{s, estimated} = \min(\Delta t_\text{CFL}, \, \Delta t_\text{ff}, \Delta t_\text{2-body})`. Hence, we have the constraint:
.. math::
\Delta t_s \leq \Delta t_\text{SF} = \eta \cfrac{M_\text{eligible} \Delta t_\text{s, estimated}}{\Delta M} \text{ if } \Delta M < 0 \; ,
where :math:`\eta` is a tolerance parameter. This parameter corresponds to ``GEARSink:tolerance_SF_timestep`` in the code.