#!/usr/bin/env python3
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi
# Date: 2023-04-12 22:41:33
"""Tools for assigning MBH."""
import numpy as np
from astropy.table import Table
from scipy.interpolate import CubicSpline
[docs]
def get_log_mbh_shankar2019(log_mstar):
"""Return equation 5 from Shankar+2019."""
x = log_mstar - 11
log_mbh = 7.574 + 1.946 * x - 0.306 * x**2 - 0.011 * x**3
return log_mbh
[docs]
def get_log_mbh_sg2016(log_mstar):
"""Return equation 3 from Shankar+2019."""
return 8.54 + 1.18 * (log_mstar - 11)
[docs]
def get_log_mbh_tanaka2024(log_mstar):
"""Return equation 10 from Tanaka+2024, see also Figure 12."""
return 0.67 * log_mstar + 1.31
[docs]
def get_log_mbh_pacucci2023(log_mstar):
"""Return the first equation from Pacucci+2023 abstract."""
return -2.43 + 1.06 * log_mstar
[docs]
def get_log_mbh_suh2020(log_mstar, alpha=1.64, beta=10.29, dalpha=0.07, dbeta=0.04, add_error=False):
"""Return the first equation of Suh+2020 Sec. 5."""
if add_error:
alpha += np.random.normal(scale=dalpha)
beta += np.random.normal(scale=dbeta)
return alpha * log_mstar - beta
[docs]
def get_log_mbh_reines_volonteri2015(log_mstar, alpha=7.45, beta=1.05):
"""Return the first equation of Reines & Volonteri 2015 Sec. 4.1."""
return alpha + beta * (log_mstar - 11)
[docs]
def get_log_mbh_kormendy_ho2013(log_mstar, alpha=8.56, beta=1.58, from_rv15=False):
"""Return Kormendy & Ho 2013.
If from rv15 = True, then return the "scaled" K&H13 relation from their
Fig. 8
"""
if not from_rv15:
return alpha + beta * (log_mstar - 11)
# NOTE: estimated visually from RV15 data on 20251214. Values are
# approximate ones but reproduce the line in Fig. 8
c = -4.20
m = 1.17
return m * log_mstar + c
[docs]
def get_log_mbh_haring_rix2004(log_mstar, a=8.20, b=1.12, da=0.10, db=0.06, add_error=False):
"""Return Häring & Rix from Angela."""
if add_error:
a += np.random.normal(scale=da, size=log_mstar.size)
b += np.random.normal(scale=db, size=log_mstar.size)
return a + b * (log_mstar - 11)
[docs]
def get_log_mbh_const(log_mstar, A=500):
"""Return logMBH according to Haring+Rix."""
return log_mstar - np.log10(A)
[docs]
def get_log_mbh_zou2024(log_mstar, z=0.0, is_tng300=False):
"""
Return Zou+2024 MBH-Mstar(z).
The MBH-Mstar is based on a hybrid approach using the observed BHAR from
earlier Zou work combined with mergers from IllustrisTNG. The Mbh-Mstar(z)
is then built by following seeded black hole masses starting from z=4.0.
Refer to Fig. 4 and Table 1 of Zou+2024:
https://ui.adsabs.harvard.edu/abs/2024ApJ...976....6Z
"""
if not np.all(z == np.clip(z, 0, 4)):
raise ValueError("Redshift must be in [0, 4]")
assert np.all(z[0] == z)
z = z[0]
# Zou defines zbins with deltaz = 0.5
zs = np.arange(0.0, 4.0 + 0.1, 0.5)
print(zs)
zclosest = zs[np.argmin(np.abs(zs - z))]
print(f"Found {zclosest=} corresponding to {z}")
def get_mstar_mbh(table, z):
x1 = table["logMstar_min"][is_tng300::2]
x2 = table["logMstar_max"][is_tng300::2]
x = (x1 + x2) / 2
y = [table[i + is_tng300][f"logMBH(z={z:.1f})"] for i in range(0, len(table), 2)]
return np.array([x, y])
table = Table.read("data/mbh_mstar_zou2024/table1.txt", format="ascii")
x, y = get_mstar_mbh(table, z)
# remove masked values...
select = np.isfinite(y)
x = x[select]
y = y[select]
return np.interp(log_mstar, x, y, left=np.nan, right=np.nan)
[docs]
GRAHAM2023_ROWS = [
r"E/Es,e",
r"E/ES,e/S0(dust=Y)",
r"BCG∗",
r"major_mergers",
r"S",
r"S_(w/o_Circinus)",
r"S0/Es,b_(dust=Y)",
]
[docs]
def get_log_mbh_graham2023(log_mstar, which="E/Es,e"):
"""
Return a MBH relation from Graham+2023 Table 1: https://arxiv.org/pdf/2305.03242
"""
row = GRAHAM2023_ROWS.index(which)
data = np.loadtxt("data/graham2023/table1.txt", usecols=(2, 4, 5))
A, C, B = data[row, :]
return A * (log_mstar - C) + B
[docs]
def get_delta_log_mbh_shankar2019(log_mstar, low=0.0, high=np.inf):
"""
Return equation 6 from Shankar+2019.
The reference for the relation is: Shankar F. et al., 2016, MNRAS, 460, 3119
"""
delta_mbh = 0.32 - 0.1 * (log_mstar - 12)
delta_mbh = np.clip(delta_mbh, low, high)
return delta_mbh
[docs]
def get_spline(Mst, Mbh):
"""
Take the Mstar and Mbh values at z=0 and creates a cubic spline.
The cubic spline extrapolates with dy/dx of the final knot and d2y/dx2 = 0.
"""
BinWidth = 0.25
Mbins = np.arange(int(min(Mst)), int(max(Mst)) + 1 + BinWidth, BinWidth)
Marr = (Mbins[1:] + Mbins[:-1]) / 2
nbins = Marr.size
meanMbh = np.array(
[
np.log10(np.mean(10 ** Mbh[np.where((Mst >= Mbins[i]) & (Mst < Mbins[i + 1]))]))
for i in range(nbins)
]
)
spline = CubicSpline(Marr[np.isfinite(meanMbh)], meanMbh[np.isfinite(meanMbh)], bc_type="natural")
"""
Now we extend the spline so that the extrapolation is a C1 continuous
linear relation
"""
# Get First point
x0 = spline.x[0]
y0 = spline(x0)
dydx0 = spline(x0, nu=1)
newx = np.nextafter(x0, x0 - 1)
newy = y0 + dydx0 * (newx - x0)
newCoeffs = np.array([0, 0, dydx0, newy]) # [d^3y/dx^3, d^2y/dx^2, dy/dx, y]
spline.extend(newCoeffs[..., None], np.r_[x0])
# Get last point
x1 = spline.x[-1]
y1 = spline(x1)
dydx1 = spline(x1, nu=1)
newx = np.nextafter(x1, x1 + 1)
newy = y1 + dydx1 * (newx - x1)
newCoeffs = np.array([0, 0, dydx1, newy]) # [d^3y/dx^3, d^2y/dx^2, dy/dx, y]
spline.extend(newCoeffs[..., None], np.r_[x1])
return spline
[docs]
SPLINE: dict[float, CubicSpline] = {}
[docs]
FUN_LOG_MBH_CONTINUITY_NEW2 = {}
[docs]
def get_log_mbh_continuity_new2(
log_mstar, z, filename="data/AGILE_Mstar_Mbh/DECODEmeans_Mstar_Mbh_Zou2024cc.dat"
):
"""
Return logMBH according to the continuity equation with Zou+24 p(lambda).
Steps to generate the required input file:
1) run opt/AGILE_Mstar_Mbh/Modules/ERDFs/genERDFdata.py
2) run opt/AGILE_Mstar_Mbh/DECODEmeans.py
Zou2024c refers to the case of assuming Zou+2024 from the parameter maps
and extrapolating above z>4 and logM<9.5.
The extrapolation scheme assumed is simple boundary extrapolation where z
and logM are clipped to their minimum and maximum values.
"""
global FUN_LOG_MBH_CONTINUITY_NEW2
if filename not in FUN_LOG_MBH_CONTINUITY_NEW2:
with open(filename) as f:
f.readline()
line = f.readline()
z1, z2, zn = map(float, [line.split()[-1][1:-1].split(",")[i] for i in (0, 1, -1)])
zvec = np.arange(z1, zn + (z2 - z1) / 2, z2 - z1)
mvec = np.loadtxt(filename)[:, 0]
MBH = np.loadtxt(filename)[:, 1:]
from scipy.interpolate import RegularGridInterpolator
FUN_LOG_MBH_CONTINUITY_NEW2[filename] = RegularGridInterpolator((mvec, zvec), MBH)
# clip
log_mstar = np.atleast_1d(log_mstar)
z = np.atleast_1d(z)
select = (log_mstar >= 8.5) & (log_mstar < 12.5) & (z >= 0.0) & (z < 5.5)
ret = np.full_like(log_mstar, np.nan)
ret[select] = FUN_LOG_MBH_CONTINUITY_NEW2[filename]((log_mstar[select], z[select]))
return np.squeeze(ret)
[docs]
REFERENCES = (
"H&R04",
"K&H13",
"R&V15",
"S&G16",
"Shankar+16",
"Suh+19",
"Pacucci+23",
"Tanaka+24",
"continuity",
"continuity_new",
"continuity_new2",
"Zou+24",
)
[docs]
def get_log_mbh(log_mstar, reference="Shankar+16", z=0.0):
"""Return logMBH according to the scaling relation from the reference."""
x = log_mstar
return {
"H&R04": get_log_mbh_haring_rix2004(log_mstar),
"K&H13": 8.56 + 1.58 * (x - 11),
"R&V15": 7.45 + 1.05 * (x - 11),
"S&G16": 8.54 + 1.18 * (x - 11),
"Shankar+16": 7.574 + 1.946 * (x - 11) - 0.306 * ((x - 11) ** 2) - 0.011 * ((x - 11) ** 3),
"Suh+19": 1.47 * x - 8.56,
"Tanaka+24": get_log_mbh_tanaka2024(log_mstar),
"Pacucci+23": get_log_mbh_pacucci2023(log_mstar),
"continuity_new2": get_log_mbh_continuity_new2(log_mstar, z=np.zeros_like(log_mstar)),
}[reference]
[docs]
def get_occupation_fraction(logMstar):
"""
Return the black hole occupation fraction.
The occupation fraction is the percentage of host galaxies that are
expected to host SMBH, which could be below unity for dwarf galaxies
(Miller+15, Burke+25, Zou+). The updated AGN duty cycle is then the product
occupation fraction and the previous AGN duty cycle.
Parameters
----------
logMstar: float or array_like
Log10 of host galaxy stellar mass in Msun.
"""
table = Table.read("data/from_zou/table2.csv")
x = table["logm"]
y = table["median"]
select = (x.min() <= logMstar) & (logMstar <= x.max())
if not np.all(select):
raise ValueError(f"logMstar must be within {x.min()} and {x.max()}")
return np.interp(logMstar, x, y)
if __name__ == "__main__":
import matplotlib.pyplot as plt
[docs]
log_mstar = np.linspace(9.0, 12.0)
for idx, reference in enumerate(REFERENCES):
plt.plot(
log_mstar,
get_log_mbh(log_mstar, reference),
label=reference,
lw=2.0,
ls=["-", "--", "-.", ":"][idx % 4],
)
plt.legend(fontsize="x-small")
plt.xlabel(r"$\log M_{\rm star}$ [Msun]")
plt.ylabel(r"$\log M_{\rm BH}$ [Msun]")
plt.ylim(5.0, None)
plt.savefig("mbh_mstar_all.pdf")