"""A module for holding blackhole emission models.
The class defined here should never be instantiated directly, there are only
ever instantiated by the parametric/particle child classes.
BlackholesComponent is a child class of Component.
"""
import numpy as np
from unyt import G, K, Lsun, Msun, c, cm, deg, km, s, yr
from unyt.physical_constants import stefan_boltzmann_constant
from synthesizer import exceptions
from synthesizer.components.component import Component
from synthesizer.synth_warnings import deprecated, deprecation
from synthesizer.units import (
Quantity,
accepts,
deepcopy_with_shared_units,
get_quantity_unit,
)
from synthesizer.utils import (
TableFormatter,
array_to_scalar,
scalar_to_array,
)
from synthesizer.utils.precision import convert_array_dtype
[docs]
class BlackholesComponent(Component):
"""The parent class for black hole components of a galaxy.
This class contains the attributes and spectra creation methods which are
common to both parametric and particle stellar components.
This should never be instantiated directly, instead it provides the common
functionality and attributes used by the child parametric and particle
BlackHole/s classes.
Attributes:
spectra (dict, Sed)
A dictionary containing black hole spectra.
mass (np.ndarray of float):
The mass of each blackhole.
accretion_rate (np.ndarray of float):
The accretion rate of each blackhole.
epsilon (np.ndarray of float):
The radiative efficiency of the blackhole.
accretion_rate_eddington (np.ndarray of float):
The accretion rate expressed as a fraction of the Eddington
accretion rate.
inclination (np.ndarray of float):
The inclination of the blackhole disc.
spin (np.ndarray of float):
The dimensionless spin of the blackhole.
bolometric_luminosity (np.ndarray of float):
The bolometric luminosity of the blackhole.
metallicity (np.ndarray of float):
The metallicity of the blackhole which is assumed for the line
emitting regions.
Attributes (For EmissionModels):
ionisation_parameter_blr (np.ndarray of float):
The ionisation parameter of the broad line region.
hydrogen_density_blr (np.ndarray of float):
The hydrogen density of the broad line region.
covering_fraction_blr (np.ndarray of float):
The fraction of disc emission intercepted and reprocessed by the
broad line region.
velocity_dispersion_blr (np.ndarray of float):
The velocity dispersion of the broad line region.
ionisation_parameter_nlr (np.ndarray of float):
The ionisation parameter of the narrow line region.
hydrogen_density_nlr (np.ndarray of float):
The hydrogen density of the narrow line region.
covering_fraction_nlr (np.ndarray of float):
The fraction of disc emission intercepted and reprocessed by the
narrow line region.
velocity_dispersion_nlr (np.ndarray of float):
The velocity dispersion of the narrow line region.
blue_bump_temperature (np.ndarray of float):
The temperature of the big blue bump, which sets where the bump
turns over at high energies. Used by the cloudy_agn incident
grids. Can be passed as a kwarg, otherwise it is calculated from
the mass and accretion rate (see calculate_blue_bump_temperature).
theta_torus (np.ndarray of float):
The angle of the torus.
torus_fraction (np.ndarray of float):
The fraction of the torus angle to 90 degrees.
transmission_fraction_escape (np.ndarray of float):
The deterministic fraction of the disc emission that escapes,
equal to ``1 - covering_fraction_blr - covering_fraction_nlr``.
transmission_fraction_nlr (np.ndarray of float):
The deterministic fraction of the disc emission transmitted
through the NLR, equal to ``covering_fraction_nlr``.
transmission_fraction_blr (np.ndarray of float):
The deterministic fraction of the disc emission transmitted
through the BLR, equal to ``covering_fraction_blr``.
random_transmission_fraction_escape (np.ndarray of float):
A one-hot random realization of the escaping disc fraction.
random_transmission_fraction_nlr (np.ndarray of float):
A one-hot random realization of the NLR-transmitted disc fraction.
random_transmission_fraction_blr (np.ndarray of float):
A one-hot random realization of the BLR-transmitted disc fraction.
"""
# Define class level Quantity attributes
accretion_rate = Quantity("mass_rate")
inclination = Quantity("angle")
bolometric_luminosity = Quantity("luminosity")
eddington_luminosity = Quantity("luminosity_solar")
shakura_sunyaev_temperature = Quantity("temperature")
mass = Quantity("mass")
# Share units rather than deep copying unyt's whole unit registry
__deepcopy__ = deepcopy_with_shared_units
@accepts(
mass=Msun.in_base("galactic"),
accretion_rate=Msun.in_base("galactic") / yr,
inclination=deg,
bolometric_luminosity=Lsun,
hydrogen_density_blr=cm**-3,
hydrogen_density_nlr=cm**-3,
velocity_dispersion_blr=km / s,
velocity_dispersion_nlr=km / s,
theta_torus=deg,
)
def __init__(
self,
fesc,
mass=None,
accretion_rate=None,
epsilon=0.1,
accretion_rate_eddington=None,
inclination=0.0 * deg,
spin=None,
bolometric_luminosity=None,
metallicity=None,
ionisation_parameter_blr=0.1,
hydrogen_density_blr=1e9 / cm**3,
covering_fraction_blr=0.1,
velocity_dispersion_blr=2000 * km / s,
ionisation_parameter_nlr=0.01,
hydrogen_density_nlr=1e4 / cm**3,
covering_fraction_nlr=0.1,
velocity_dispersion_nlr=500 * km / s,
theta_torus=10 * deg,
**kwargs,
):
"""Initialise the BlackholeComponent.
Where they're not provided missing quantities are automatically
calculated. Not all parameters need to be set for every emission model.
Args:
fesc (float):
The escape fraction of the blackhole.
mass (np.ndarray of float):
The mass of each blackhole.
accretion_rate (np.ndarray of float):
The accretion rate of each blackhole.
epsilon (np.ndarray of float):
The radiative efficiency of the blackhole.
accretion_rate_eddington (np.ndarray of float):
The accretion rate expressed as a fraction of the Eddington
accretion rate.
inclination (np.ndarray of float):
The inclination of the blackhole disc.
spin (np.ndarray of float):
The dimensionless spin of the blackhole.
bolometric_luminosity (np.ndarray of float):
The bolometric luminosity of the blackhole.
metallicity (np.ndarray of float):
The metallicity of the blackhole which is assumed for the line
emitting regions.
ionisation_parameter_blr (np.ndarray of float):
The ionisation parameter of the broad line region.
hydrogen_density_blr (np.ndarray of float):
The hydrogen density of the broad line region.
covering_fraction_blr (np.ndarray of float):
The fraction of disc emission intercepted and reprocessed by
the broad line region.
velocity_dispersion_blr (np.ndarray of float):
The velocity dispersion of the broad line region.
ionisation_parameter_nlr (np.ndarray of float):
The ionisation parameter of the narrow line region.
hydrogen_density_nlr (np.ndarray of float):
The hydrogen density of the narrow line region.
covering_fraction_nlr (np.ndarray of float):
The fraction of disc emission intercepted and reprocessed by
the narrow line region.
velocity_dispersion_nlr (np.ndarray of float):
The velocity dispersion of the narrow line region.
theta_torus (np.ndarray of float):
The angle of the torus.
**kwargs (dict):
Any other parameter for the emission models can be provided as
kwargs.
"""
# Initialise the parent class
Component.__init__(self, "BlackHoles", fesc, **kwargs)
# Save the black hole properties
self.mass = mass
self.accretion_rate = accretion_rate
self.epsilon = epsilon
self.accretion_rate_eddington = accretion_rate_eddington
self.spin = spin
self.bolometric_luminosity = bolometric_luminosity
self.metallicity = metallicity
# Below we attach all the possible attributes that could be needed by
# the emission models.
# Set BLR attributes
self.ionisation_parameter_blr = ionisation_parameter_blr
self.hydrogen_density_blr = hydrogen_density_blr
self.covering_fraction_blr = covering_fraction_blr
self.velocity_dispersion_blr = velocity_dispersion_blr
# Set NLR attributes
self.ionisation_parameter_nlr = ionisation_parameter_nlr
self.hydrogen_density_nlr = hydrogen_density_nlr
self.covering_fraction_nlr = covering_fraction_nlr
self.velocity_dispersion_nlr = velocity_dispersion_nlr
self.covering_fraction = None
self.escape_fraction = None
self.transmission_fraction_escape = None
self.transmission_fraction_nlr = None
self.transmission_fraction_blr = None
self.random_transmission_fraction_escape = None
self.random_transmission_fraction_nlr = None
self.random_transmission_fraction_blr = None
# If covering fractions are provided, compute both the deterministic
# transmission weights and a random one-hot realization used by the
# UnifiedAGN emission model.
if self.covering_fraction_blr is not None:
# Convert the covering_fraction_blr and covering_fraction_nlr to
# arrays to allow us to use the same logic for both parametric and
# particle blackholes.
covering_fraction_blr = scalar_to_array(covering_fraction_blr)
# Validate both covering fractions are set and have matching
# lengths.
if covering_fraction_nlr is None:
raise exceptions.InconsistentArguments(
"covering_fraction_nlr must be provided when "
"covering_fraction_blr is set"
)
covering_fraction_nlr = scalar_to_array(covering_fraction_nlr)
if len(covering_fraction_blr) != len(covering_fraction_nlr):
raise exceptions.InconsistentArguments(
"covering_fraction_blr "
f"(length {len(covering_fraction_blr)}) and "
"covering_fraction_nlr "
f"(length {len(covering_fraction_nlr)}) must have the"
"same length"
)
# Validate each covering fraction before combining them.
if np.any(
(covering_fraction_blr < 0.0) | (covering_fraction_blr > 1.0)
):
raise exceptions.InconsistentArguments(
"covering_fraction_blr must be between 0.0 and 1.0"
)
if np.any(
(covering_fraction_nlr < 0.0) | (covering_fraction_nlr > 1.0)
):
raise exceptions.InconsistentArguments(
"covering_fraction_nlr must be between 0.0 and 1.0"
)
# Calculate the total covering fraction and corresponding escape
# fraction.
self.covering_fraction = (
covering_fraction_blr + covering_fraction_nlr
)
self.escape_fraction = 1.0 - self.covering_fraction
# Validate that covering fractions don't exceed unity.
if np.any(self.covering_fraction > 1.0):
raise exceptions.InconsistentArguments(
"Sum of BLR and NLR covering fractions cannot exceed 1.0"
)
# Define number of blackholes.
N = len(covering_fraction_blr)
# Deterministic transmission fractions mirror the physical
# covering fractions on the component.
transmission_fraction_escape = self.escape_fraction
transmission_fraction_nlr = covering_fraction_nlr
transmission_fraction_blr = covering_fraction_blr
# Loop over blackholes and decide whether the disc emission
# escapes (none), or is transmitted through the BLR (blr) or NLR
# (nlr).
random_transmission_index = np.empty(N, dtype=int)
for i, (
blr_covering_fraction_,
nlr_covering_fraction_,
) in enumerate(zip(covering_fraction_blr, covering_fraction_nlr)):
# Define the probabilities for each option based on the
# covering fractions.
probabilities = [
blr_covering_fraction_,
nlr_covering_fraction_,
1.0 - blr_covering_fraction_ - nlr_covering_fraction_,
]
# Randomly choose the scenario using these probabilities.
random_transmission_index[i] = np.random.choice(
3,
p=probabilities,
)
# Random transmission fractions are a one-hot realization of the
# three possible transmission paths.
random_transmission_fraction_escape = np.zeros(N)
random_transmission_fraction_nlr = np.zeros(N)
random_transmission_fraction_blr = np.zeros(N)
# For each corresponding scenario set the transmission fraction to
# unity.
random_transmission_fraction_escape[
random_transmission_index == 2
] = 1.0
random_transmission_fraction_nlr[
random_transmission_index == 1
] = 1.0
random_transmission_fraction_blr[
random_transmission_index == 0
] = 1.0
# convert to scalars if only one value
if N == 1:
self.transmission_fraction_escape = array_to_scalar(
transmission_fraction_escape
)
self.transmission_fraction_nlr = array_to_scalar(
transmission_fraction_nlr
)
self.transmission_fraction_blr = array_to_scalar(
transmission_fraction_blr
)
self.random_transmission_fraction_escape = array_to_scalar(
random_transmission_fraction_escape
)
self.random_transmission_fraction_nlr = array_to_scalar(
random_transmission_fraction_nlr
)
self.random_transmission_fraction_blr = array_to_scalar(
random_transmission_fraction_blr
)
else:
self.transmission_fraction_escape = (
transmission_fraction_escape
)
self.transmission_fraction_nlr = transmission_fraction_nlr
self.transmission_fraction_blr = transmission_fraction_blr
self.random_transmission_fraction_escape = (
random_transmission_fraction_escape
)
self.random_transmission_fraction_nlr = (
random_transmission_fraction_nlr
)
self.random_transmission_fraction_blr = (
random_transmission_fraction_blr
)
# The inclination of the black hole disc
self.inclination = (
inclination if inclination is not None else 0.0 * deg
)
# The angle of the torus
self.theta_torus = theta_torus
self.torus_fraction = (self.theta_torus / (90 * deg)).value
# Check to make sure that both accretion rate and bolometric luminosity
# haven't been provided because that could be confusing.
if (self.accretion_rate is not None) and (
self.bolometric_luminosity is not None
):
raise exceptions.InconsistentArguments(
"""Both accretion rate and bolometric luminosity provided but
that is confusing. Provide one or the other!"""
)
if (self.accretion_rate_eddington is not None) and (
self.bolometric_luminosity is not None
):
raise exceptions.InconsistentArguments(
"""Both accretion rate (in terms of Eddington) and bolometric
luminosity provided but that is confusing. Provide one or
the other!"""
)
# Ensure that both accretion_rate and accretion_rate_eddington are not
# provided, this is also confusing.
if (self.accretion_rate is not None) and (
self.accretion_rate_eddington is not None
):
raise exceptions.InconsistentArguments(
"""Both accretion rate and accretion rate in terms of
Eddington provided but that is confusing. Provide one or the
other!"""
)
# If mass calculate the Eddington luminosity.
if self.mass is not None:
self.calculate_eddington_luminosity()
# If accretion_rate_eddington provided calculate the accretion rate.
if (
self.accretion_rate_eddington is not None
and self.eddington_luminosity is not None
):
self.calculate_accretion_rate()
# If mass, accretion_rate, and epsilon provided calculate the
# bolometric luminosity.
if (
self.mass is not None
and self.accretion_rate is not None
and self.epsilon is not None
):
self.calculate_bolometric_luminosity()
# If mass, accretion_rate, and epsilon provided calculate the
# Shakura-Sunyaev temperature.
if (
self.mass is not None
and self.accretion_rate is not None
and self.epsilon is not None
):
self.calculate_shakura_sunyaev_temperature()
# If mass and accretion_rate provided, and no blue_bump_temperature
# was passed as a kwarg (singular or plural), calculate the blue bump
# temperature.
if (
self.mass is not None
and self.accretion_rate is not None
and getattr(self, "blue_bump_temperature", None) is None
and getattr(self, "blue_bump_temperatures", None) is None
):
self.calculate_blue_bump_temperature()
# If mass, accretion_rate, and epsilon provided calculate the
# Eddington ratio.
if (
self.mass is not None
and self.accretion_rate is not None
and self.epsilon is not None
):
self.calculate_eddington_ratio()
# If mass, accretion_rate, and epsilon provided calculate the
# accretion rate in units of the Eddington accretion rate. This is the
# bolometric_luminosity / eddington_luminosity.
if (
self.mass is not None
and self.accretion_rate is not None
and self.epsilon is not None
):
self.calculate_accretion_rate_eddington()
# If inclination, calculate the cosine of the inclination, required by
# some models (e.g. AGNSED).
if self.inclination is not None:
self.cosine_inclination = np.cos(
self.inclination.to_value("radian")
)
[docs]
def calculate_accretion_rate(self):
"""Calculate the black hole accretion rate from the eddington ratio.
This is used when accretion rate is provided in terms of the Eddington
ratio and not the explicit accretion rate.
Returns:
unyt_array:
The black hole accretion rate
"""
self.accretion_rate = convert_array_dtype(
self.accretion_rate_eddington
* self.eddington_luminosity
/ (self.epsilon * c**2),
self.mass.dtype,
)
return self.accretion_rate
[docs]
def calculate_bolometric_luminosity(self):
"""Calculate the black hole bolometric luminosity.
Returns:
unyt_array:
The black hole bolometric luminosity
"""
# Derived properties keep the precision of the masses
self.bolometric_luminosity = convert_array_dtype(
(self.epsilon * self.accretion_rate * c**2).to(
get_quantity_unit(self, "bolometric_luminosity"),
),
self.mass.dtype,
)
return self.bolometric_luminosity
[docs]
def calculate_eddington_luminosity(self):
"""Calculate the eddington luminosity of the black hole.
Returns:
unyt_array
The black hole eddington luminosity in solar luminosities
"""
# The Eddington luminosity is given by:
# L_Edd = 4*pi*G*mp*c*M/sigma_thompson = 1.257e38 * M/Msun erg/s
# Converting to solar luminosities:
# L_Edd = 1.257e38 / 3.828e33 = 3.284e4 Lsun/Msun
self.eddington_luminosity = 3.284e4 * self._mass * Lsun
return self.eddington_luminosity
[docs]
def calculate_eddington_ratio(self):
"""Calculate the eddington ratio of the black hole.
Returns:
unyt_array:
The black hole eddington ratio
"""
# Dividing the luminosities cancels their units. NOTE: to_value makes
# a copy, which is fine for these small per-blackhole arrays.
self.eddington_ratio = convert_array_dtype(
(self.bolometric_luminosity / self.eddington_luminosity).to_value(
"dimensionless",
),
self.mass.dtype,
)
return self.eddington_ratio
[docs]
def calculate_shakura_sunyaev_temperature(self):
"""Calculate the black hole Shakura-Sunyaev temperature.
This is the characteristic temperature of a Shakura and Sunyaev
(1973) thin accretion disc with its inner edge at 6 GM/c^2,
T* = [3 G M Mdot / (8 pi sigma r_in^3)]^(1/4).
Returns:
unyt_array:
The black hole Shakura-Sunyaev temperature
"""
# Calculate the Shakura-Sunyaev temperature
self.shakura_sunyaev_temperature = (
2.24e9 * self._accretion_rate ** (1 / 4) * self._mass**-0.5
)
return self.shakura_sunyaev_temperature
@property
def bb_temperature(self):
"""Get the Shakura-Sunyaev temperature.
Deprecated alias for shakura_sunyaev_temperature.
Returns:
unyt_array:
The black hole Shakura-Sunyaev temperature.
"""
deprecation(
"`bb_temperature` is deprecated and will be removed in version "
"1.4.0. Use `shakura_sunyaev_temperature` instead."
)
return self.shakura_sunyaev_temperature
@bb_temperature.setter
def bb_temperature(self, value):
"""Set the Shakura-Sunyaev temperature.
Deprecated alias for shakura_sunyaev_temperature.
Args:
value (unyt_array):
The black hole Shakura-Sunyaev temperature.
"""
deprecation(
"`bb_temperature` is deprecated and will be removed in version "
"1.4.0. Use `shakura_sunyaev_temperature` instead."
)
self.shakura_sunyaev_temperature = value
[docs]
@deprecated(
"is deprecated and will be removed in version 1.4.0. "
"Use calculate_shakura_sunyaev_temperature instead."
)
def calculate_bb_temperature(self):
"""Calculate the Shakura-Sunyaev temperature.
Deprecated alias for calculate_shakura_sunyaev_temperature.
Returns:
unyt_array:
The black hole Shakura-Sunyaev temperature.
"""
return self.calculate_shakura_sunyaev_temperature()
[docs]
def calculate_blue_bump_temperature(self):
"""Calculate the black hole blue bump temperature.
This is the temperature of the big blue bump in the AGN continuum,
f_nu ~ nu^alpha_uv exp(-h nu / k T), so it sets where the bump turns
over at high energies. It is the first parameter of cloudy's agn
command, and an axis of the cloudy_agn incident grids in syncretize.
The temperature is calculated following Inayoshi and Maiolino (2025)
(https://ui.adsabs.harvard.edu/abs/arXiv:2409.07805/abstract).
Returns:
unyt_array:
The black hole blue bump temperature.
"""
# The radius and the inner edge of the disc
radius = 20 * G * self.mass / c**2
inner_radius = 6 * G * self.mass / c**2
# Calculate the temperature of the disc at this radius
temperature = (
3
* G
* self.mass
* self.accretion_rate
/ (8 * np.pi * stefan_boltzmann_constant * radius**3)
* (1 - np.sqrt(inner_radius / radius))
) ** (1 / 4)
# Derived properties keep the precision of the masses
self.blue_bump_temperature = convert_array_dtype(
temperature.to(K),
self.mass.dtype,
)
return self.blue_bump_temperature
[docs]
def calculate_accretion_rate_eddington(self):
"""Calculate the black hole accretion in units of the Eddington rate.
Returns:
unyt_array
The black hole accretion rate in units of the Eddington rate.
"""
# Dividing the luminosities cancels their units. NOTE: to_value makes
# a copy, which is fine for these small per-blackhole arrays.
self.accretion_rate_eddington = convert_array_dtype(
(self.bolometric_luminosity / self.eddington_luminosity).to_value(
"dimensionless",
),
self.mass.dtype,
)
return self.accretion_rate_eddington
def __str__(self):
"""Return a string representation of the particle object.
Returns:
table (str):
A string representation of the particle object.
"""
# Intialise the table formatter
formatter = TableFormatter(self)
return formatter.get_table("Black Holes")
[docs]
def calculate_circular_velocity(self, radial_distance):
r"""Calculate the circular velocity.
v_c(r) = \\sqrt{\frac{G M(<r)}{r}}.
Args:
radial_distance (np.ndarray of float):
The distance from the blackhole.
Returns:
unyt_array:
The circular velocity at the radial distance.
"""
return np.sqrt(G * self.mass / radial_distance)
[docs]
def calculate_schwarzschild_radius(self):
r"""Calculate the Schwarzschild radius of the blackhole(s).
R_{\rm s} = \frac{2GM}{c^2}
Returns:
unyt_array:
The Schwarzschild radius.
"""
return 2 * G * self.mass / c**2
[docs]
def calculate_integrated_ionising_luminosity(self):
"""Calculates the integrated ionising luminosity of the blackhole(s).
This requires that the disc_incident spectra be available.
Returns:
unyt_array:
The ionising photon production rate (s^-1).
"""
if "disc_incident" in self.spectra.keys():
return self.spectra[
"disc_incident"
].calculate_ionising_photon_production_rate()
else:
raise exceptions.MissingSpectraType(
"It is necessary to first calculate the disc_incident spectra "
"before calculating the ionising luminosity"
)
[docs]
def calculate_ionisation_parameter(
self, radial_distance, hydrogen_density
):
r"""Calculate the ionisation parameter.
Calculates the ionising parameter (U) at radial_distance for the given
hydrogen_density.
U = \frac{Q_{\mathrm{H}}}{4\pi r^{2} n_{\mathrm{H}} c}
Args:
radial_distance (np.ndarray of float):
The distance from the blackhole.
hydrogen_density (np.ndarray of float):
The hydrogen density at radial_distance.
Returns:
unyt_array:
The ionisation parameter.
"""
ionising_luminosity = self.calculate_ionising_luminosity()
return ionising_luminosity / (
4 * np.pi * radial_distance**2 * hydrogen_density * c
)