In [2]:
import numpy as np
import matplotlib.pyplot as plt
mpl.rcParams['font.size'] = 10
In [23]:
G = 6.67259e-8 # cm^3 s^-2 g^-1
M6 = 1.e6 * 1.989e33 # g [10^8 solar masses]
t100 = 100.0 * 3.156e7 # s [100 year timescale for ghost halo]
pc_to_cm = 3.0857e18 # cm / pc
km_to_cm = 1.e5 # cm / km
1. Assuming mass is virialized
$$ V^2 = \frac{GM}{R} $$2. Associated disk radius
$$ V = \frac{2 \pi R}{T_{\rm circ}} $$Therefore
$$ R^3 = G M \left( \frac{T_{\rm circ}}{2 \pi} \right)^2 $$
In [24]:
## Calculate it
gm_term = G * M6 # cm^3 s^-2
t_term = t100 / (2*np.pi) # s
r_cubed = gm_term * t_term**2 # cm^3
Rtest = np.power(r_cubed, 1./3.) / pc_to_cm # pc
In [25]:
print Rtest
So the final result is
$$ R = 0.01\ {\rm pc}\ \left( \frac{M}{10^6 M_\odot} \right)^{1/3}\ \left( \frac{\delta t}{100\ {\rm years}} \right)^{2/3} $$What is the velocity of this region?
In [26]:
velocity = np.power(G*M6/(Rtest*pc_to_cm), 0.5) # cm/s
velocity /= km_to_cm # km/s
print velocity
So that's $$ V = 3000\ {\rm km/s}\ \left( \frac{M}{10^8 M_{\odot}} \right)^{1/2}\ \left( \frac{R}{0.05\ {\rm pc}} \right)^{-1/2} $$
In [27]:
# Radius associated with t_max timescale [pc]
def calc_R( m, t ):
gm = G * M6 * m
tt = t100*t / (2.0*np.pi)
rc = gm * tt**2
return np.power(rc, 1./3.) / pc_to_cm
In [28]:
mgrid = np.logspace(0, 3, 100) # 10^6 to 10^9 Msun SMBH
plt.plot(mgrid*1.e8, calc_R(mgrid, 1.0), 'k-', lw=2 )
plt.semilogx()
plt.xlabel(r"$M_{SMBH}$ [$M_{\odot}$]")
plt.ylabel(r"$R$ [pc]")
Out[28]:
In [29]:
tgrid = np.logspace(-2, 1, 100) # 1 to 1000 years
plt.plot(tgrid*100.0, calc_R(1.0, tgrid), 'k-', lw=2)
plt.semilogx()
plt.xlabel(r"$\delta t$ [years]")
plt.ylabel(r"$R$ [pc]")
Out[29]:
Taking the famous "alpha prescription", an estimate for the disk viscosity is
$$ \nu = \alpha c_s h\ \ (1) $$where $\alpha \leq 1$ (Eq 4.34, Accretion Power in Astrophysics)
The other quantities are $$ c_s^2 = \frac{kT}{m_p}\ \ (2) $$ and $$ h^2 = \frac{kT R^3}{G M m_p}\ \ (3) $$ (with help from Section 5.7.2 of Choudhuri, Physics of Fluids and Plasmas)
The "viscous timescale" -- or amount of time it takes for a "disc annulus to move a radial distance $R$" -- is $$ t_{\rm visc} \sim R^2 / \nu $$ (Eq 5.13, Accretion Power in Astrophysics)
In [30]:
## Constants for calculation
k = 1.380658e-16 # erg/K
T0 = 1.e7 # K
mp = 1.6726231e-24 # g
R0 = 3.0857e18 # cm (1 pc)
yr = 3.e7 # s/yr
In [31]:
def cs( T=1.0 ):
return np.sqrt( k * T0 * T / mp )
def h( M=1.0, R=1.0, T=1.0 ):
ftop = k * T0*T * (R0*R)**3
fbot = G * M6*M * mp
return np.sqrt( ftop / fbot )
ALPHA = 0.1
def visc( M=1.0, R=1.0, T=1.0, alpha=ALPHA ):
Cs = cs(T)
H = h(M, R, T)
return alpha * Cs * H
def tvisc( M=1.0, R=1.0, T=1.0, alpha=ALPHA ):
nu = visc(M, R, T, alpha)
return (R0*R)**2 / nu
In [32]:
## Get fiducial values
print "Fiducial c_s", cs()*1.e-5, "km/s"
print "Fiducial h", h()/R0, "pc"
print "Fiducial viscosity", visc()
print "Fiducial t_visc", tvisc()/yr, "years"
Taking the fiducial values above:
$$ c_s = 300\ \left(\frac{T}{10^7\ {\rm K}}\right)^{1/2}\ {\rm km/s}$$$$ h = 4\ \left(\frac{T}{10^7\ {\rm K}}\right)^{1/2} \left(\frac{R}{{\rm pc}}\right)^{3/2} \left(\frac{M}{10^6\ M_\odot}\right)^{-1/2}\ {\rm pc} $$$$ t_{\rm visc} = \frac{8000}{\alpha}\ \left(\frac{M}{10^6\ M_\odot}\right)^{1/2} \left(\frac{R}{{\rm pc}}\right)^{1/2} \left(\frac{T}{10^6\ {\rm K}}\right)^{-1}\ {\rm years}$$
In [38]:
## Calculate tvisc for a number of black hole masses
rad = np.logspace(-7, 0, 100)
MBH = np.array([0.1, 1.0, 10.0, 100.0]) # 10^6 - 10^8 Msun
col = dict(zip(MBH, ['0.0','0.3','0.4','0.6']))
labels = dict(zip(MBH, ["$10^5$ M$_\odot$","$10^6$ M$_\odot$","$10^7$ M$_\odot$","$10^8$ M$_\odot$"]))
for MM in MBH:
plt.plot(np.log10(rad), np.log10(tvisc(MM, rad)/yr), lw=2, color=col[MM], label=labels[MM])
plt.legend(loc="upper left", frameon=False)
plt.ylim(0,6)
rtmax = calc_R(MBH/100.0, 1.0)
tvmax = tvisc(MBH,rtmax)/yr
plt.plot( np.log10(rtmax), np.log10(tvmax), 'ko' )
plt.axhline(np.log10(t100/yr), color='k', ls=':')
plt.xlabel(r"$\log$ Radius [pc]", size=12)
plt.ylabel(r"$\log\ t_{\rm visc}\ $ [years]", size=12)
#plt.savefig('variability_timescales_fig1.pdf', format='pdf')
The timescale for the break frequency in the PSD is likely associated with the inner edge of the disk (McHardy 2010, L.N.Ph.) $$ T_B = 1/\nu_B \propto M/\dot{m}_E $$
Using NGC 4051 as template, which has a break frequency $\nu_B \sim 8 \times 10^{-4}$ and $\dot{m}_E \sim 0.3$ $$ \nu_B = 800\ \dot{m}_E \left(\frac{M}{M_\odot}\right)^{-1}\ {\rm Hz} $$ or
In [34]:
1.0 / (800.0/1.e6)
Out[34]:
What radius does that correspond to?
In [35]:
RB = calc_R( 0.01, 1250./t100 ) #pc
print(RB)
That's several schwarzschild radii, for a $10^6\ M_\odot$ black hole ($r_s \approx 10^{-7}\ M_6$ pc).
$$ R_B \approx 6\ r_s\ \dot{m}_E^{-2/3} $$
In [36]:
def rs( m6 ):
return 1.e-7 * m6 # pc
In [39]:
## Redo above plot in units of schwarzschild radii
rs_MBH = dict(zip(MBH,rs(MBH)))
for MM in MBH:
plt.plot(np.log10(rad/rs_MBH[MM]), np.log10(tvisc(MM, rad)/yr), lw=2, color=col[MM], label=labels[MM])
plt.legend(loc="lower right", frameon=False)
plt.xlim(0,8)
rtmax = calc_R(MBH, 1.0)
tvmax = tvisc(MBH,rtmax)/yr
plt.plot( np.log10(rtmax/rs(MBH)), np.log10(tvmax), 'ko' )
plt.axhline(np.log10(t100/yr), color='k', ls=':')
rb_sub_edd = 6.0 * np.power(0.3, -2./3.)
plt.axvline(np.log10(6.0), color='r', ls='-')
plt.axvline(np.log10(rb_sub_edd), color='r', ls='--')
plt.xlabel(r"$\log\ r_s$", size=12)
plt.ylabel(r"$\log\ t_{\rm visc}\ $ [years]", size=12)
#plt.savefig('variability_timescales_fig2.pdf', format='pdf')
In [19]:
eV = 1.602e-12 # erg/eV
print "Temp (1 keV) ~", 1.e3 * eV / k, "K"
So this should be a good rule of thumb to remember
$$ T \approx 10^7\ K\ \left(\frac{E}{{\rm keV}}\right) $$
In [19]: