notebook.community

Edit and run

Note:

Need to check the calcuation of X-ray luminosity Still don't know what the density in Freeland & Wilcots 2011 means. total number? electron? ion? Assuming the number is total number density.

Galaxy Group Profiles for Jet Simulations





In [ ]:



    


%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
from ipywidgets import interact, FloatSlider
from scipy.special import hyp2f1
from scipy.special import gamma

kpc = 3.086e+21 #cm
mH = 1.6737352238051868e-24 #g
Msun = 2E33 #g

sim_mu			= 0.60			# mean molecular weight

# Ratio of number of electrons to all particles
# 75% Hydrogen, 25% Helium
e_ratio = 0.875/1.688
# Data from Table 1 of Freeland & Wilcots 2011
R_FW2011 = np.array([300, 2000, 270, 700, 15])
n_IGM_FW2011 = np.array([3E-3, 5E-4, 2E-4, 5E-4, 2E-3])
n_IGM_err_FW2011 = np.array([2E-3, 4E-4, 1E-4, 2E-4, 1E-3])

ne_FW2011 = n_IGM_FW2011*e_ratio
ne_err_FW2011 = n_IGM_err_FW2011*e_ratio

# Ratio of number of electrons to hydrogen
eH_ratio = 0.875/0.75

# Data point from Pisano et al., 2004
R_Pisano2004 = np.array([650])
nH_Pisano2004 = np.array([10**(-5.2)])

ne_Pisano2004 = nH_Pisano2004*eH_ratio

$\beta$-model of density

$$ \rho(r) = \frac{\rho_0}{[1+(\frac{r}{R_c})^2]^{\frac{3}{2}\beta}} $$

Integration of the square of the above 3D profile along line-of-sight will give x-ray surface brightness profile

$$ S(R) = S_0 [ 1 + (\frac{R}{R_c})^2 ]^{-3\beta+0.5} $$

in which

$$ S_0 = 2 ~c_0 ~n_0^2 ~\Lambda(T) $$

$ \Lambda(T)$ is the x-ray cooling function of temperature T, $ c_0 $ is a constant depending on $\beta$.

$ c_0 = \frac{\frac{\sqrt{\pi}}{2}\Gamma(3\beta-\frac{1}{2})}{\Gamma(3\beta)} $

Total x-ray luminosity

The x-ray luminosity is calculated from the integration of X-ray surface brightness $$ \begin{align*} L_x &= \int_0^{r_{max}} S(r)~ 2\pi r ~\mathrm{d} r \\ &= 2\pi c_0~n_0^2 R_c^3 \Lambda \frac{1}{-3\beta+3/2} \left\{ \left[ 1+(\frac{r_{max}}{R_c})^2 \right]^{-3\beta+3/2} -1 \right\} &\textrm{if} \quad \beta \neq 0.5 \\ &= 2\pi c_0~n_0^2 R_c^3 \Lambda \frac{1}{-3\beta+3/2} \left\{ \ln \left[ 1+(\frac{r_{max}}{R_c})^2 \right]\right\} &\textrm{if} \quad \beta = 0.5 \end{align*} $$





In [ ]:



    


def Lx(rmax, rc, beta, n0, metallicity=0.5):
    Lambda = 2.95E-23*metallicity + 7.86E-24
    # X-ray 0.1~2.4keV metal cooling rate = 2.95E-23 erg cm^3 s^{-1}
    #             primordial cooling rate = 7.86E-24 erg cm^3 s^{-1}
    L = 2*np.pi*n0**2*(rc*kpc)**3*Lambda*np.sqrt(np.pi)/2*gamma(3*beta-0.5)/gamma(3*beta)
    if beta == 0.5:
        return L*np.log(1+(rmax/rc)**2)
    else:
        return L*1/(1.5-3*beta)*((1+(rmax/rc)**2)**(-3*beta+1.5)-1)

Two $\beta$-models





In [ ]:



    


rhoCore1	= 1E-3
beta1 		= 0.5
rCore1		= 30

rhoCore2	= 1
beta2 		= 0.8
rCore2		= 1

rr = np.logspace(-1, 3, 100)

density1 = rhoCore1*(1.0 + (rr/rCore1)**2)**(-1.5*beta1)
density2 = rhoCore2*(1.0 + (rr/rCore2)**2)**(-1.5*beta2)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

line1, = ax.plot(rr, density1, ls='--', 
         label=r'group $r_c=%i$kpc, $\rho_c$=%.3fcm$^{-3}$ $\beta$=%.1f' % (rCore1, rhoCore1, beta1))
line2, = ax.plot(rr, density2, ls=':', 
         label=r'galaxy $r_c=%i$kpc, $\rho_c=%.3f$cm$^{-3}$ $\beta$=%.1f' % (rCore2, rhoCore2, beta2))
line3, = ax.plot(rr, (density1+density2), label='total density')
Mgas1 = rhoCore1*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta1, 2.5, -(rr[-1]/rCore1)**2)
text1 = ax.text(0.6, 0.1, r'$M_{gas,group}=%.2e M_{\odot}$' % Mgas1, transform=ax.transAxes)
Mgas2 = rhoCore2*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta2, 2.5, -(rr[-1]/rCore2)**2)
text2 = ax.text(0.6, 0.03, r'$M_{gas,galaxy}=%.2e M_{\odot}$' % Mgas2, transform=ax.transAxes)

ax.errorbar(R_FW2011, n_IGM_FW2011, yerr=n_IGM_err_FW2011, fmt='.', capsize=5, label='Freeland & Wilcots 2011', c='k')

#ax.legend(loc=1)
ax.set_xlim(rr[0], rr[-1])
plt.semilogy()
plt.semilogx()
plt.grid(ls=':')
plt.xlabel('radius (kpc)')

def update(log_rhoCore1=-3, beta1=0.5, rCore1=30,
           log_rhoCore2=0, beta2=0.8, rCore2=1, log_rmax=3.5):
    rhoCore1 = 10**log_rhoCore1
    rhoCore2 = 10**log_rhoCore2
    rr = np.logspace(-1, log_rmax, 100)
    density1 = rhoCore1*(1.0 + (rr/rCore1)**2)**(-1.5*beta1)
    density2 = rhoCore2*(1.0 + (rr/rCore2)**2)**(-1.5*beta2)
    line1.set_xdata(rr)
    line1.set_ydata(density1)
    label1 = r'group $r_c=%i$kpc, $\rho_c$=%.3fcm$^{-3}$ $\beta$=%.2f' % (rCore1, rhoCore1, beta1)
    line1.set_label(label1)
    
    line2.set_xdata(rr)
    line2.set_ydata(density2)
    label2 = r'galaxy $r_c=%i$kpc, $\rho_c$=%.3fcm$^{-3}$ $\beta$=%.2f' % (rCore2, rhoCore2, beta2)
    line2.set_label(label2)
    
    line3.set_xdata(rr)
    line3.set_ydata((density1+density2))
    Mgas1 = rhoCore1*mH/Msun*kpc**3*4*np.pi*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta1, 2.5, -(rr[-1]/rCore1)**2)
    text1.set_text(r'$M_{gas,group}=%.2e M_{\odot}$' % Mgas1)
    Mgas2 = rhoCore2*mH/Msun*kpc**3*4*np.pi*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta2, 2.5, -(rr[-1]/rCore2)**2)
    text2.set_text(r'$M_{gas,galaxy}=%.2e M_{\odot}$' % Mgas2)
    ax.set_xlim(rr[0], rr[-1])
    ax.legend(loc=3, fontsize=8)
    fig.canvas.draw()

interact(update, log_rhoCore1=(-4, -1, 0.1), beta1=(0.0, 1, 0.01), rCore1=(1,100,1),
                 log_rhoCore2=(-2, 2, 0.1), beta2=(0.3, 1, 0.01), rCore2=(1,100,1), log_rmax=(2, 4, 0.1));

Use one $\beta$-model to fit two $\beta$-models





In [ ]:



    


rhoCore1	= 1E-3
beta1 		= 0.5
rCore1		= 30

rhoCore2	= 1
beta2 		= 0.8
rCore2		= 1

rr = np.logspace(-1, 3, 100)

density1 = rhoCore1*(1.0 + (rr/rCore1)**2)**(-1.5*beta1)
density2 = rhoCore2*(1.0 + (rr/rCore2)**2)**(-1.5*beta2)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

line1, = ax.plot(rr, density1, ls='--', 
         label=r'group $r_c=%i$kpc, $\rho_c$=%.3fcm$^{-3}$ $\beta$=%.1f' % (rCore1, rhoCore1, beta1))
line2, = ax.plot(rr, density2, ls=':', 
         label=r'galaxy $r_c=%i$kpc, $\rho_c=%.3f$cm$^{-3}$ $\beta$=%.1f' % (rCore2, rhoCore2, beta2))
line3, = ax.plot(rr, (density1+density2), label='total density')
Mgas1 = rhoCore1*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta1, 2.5, -(rr[-1]/rCore1)**2)
text1 = ax.text(0.6, 0.93, r'$M_{gas,group}=%.2e M_{\odot}$' % Mgas1, transform=ax.transAxes)
Mgas2 = rhoCore2*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta2, 2.5, -(rr[-1]/rCore2)**2)
text2 = ax.text(0.6, 0.85, r'$M_{gas,galaxy}=%.2e M_{\odot}$' % Mgas2, transform=ax.transAxes)

rhoCore, rCore, beta = (0.5, 2, 0.6)

density = rhoCore*(1.0 + (rr/rCore)**2)**(-1.5*beta)
label = r'$\beta$-model $r_c=%i$kpc, $\rho_c=%.3f$cm$^{-3}$ $\beta$=%.1f' % (rCore, rhoCore, beta)
line, = ax.plot(rr, density, label=label)
Mgas = rhoCore*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta, 2.5, -(rr[-1]/rCore)**2)
text = ax.text(0.6, 0.77, r'$M_{gas,\beta}=%.2e M_{\odot}$' % Mgas, transform=ax.transAxes)

ax.errorbar(R_FW2011, n_IGM_FW2011, yerr=n_IGM_err_FW2011, fmt='.', capsize=5, label='Freeland & Wilcots 2011', c='k')

ax.legend(fontsize=7)
ax.set_xlim(rr[0], rr[-1])
plt.semilogy()
plt.semilogx()
plt.grid(ls=':')
plt.xlabel('radius (kpc)')

def update(log_rhoCore=0, beta=0.6, rCore=1):
    rhoCore = 10**log_rhoCore
    density = rhoCore*(1.0 + (rr/rCore)**2)**(-1.5*beta)
    line.set_ydata(density)
    label = r'$\beta$-model $r_c=%i$kpc, $\rho_c=%.3f$cm$^{-3}$ $\beta$=%.1f' % (rCore, rhoCore, beta)
    line.set_label(label)
    
    Mgas = rhoCore*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta, 2.5, -(rr[-1]/rCore)**2)
    text.set_text(r'$M_{gas,\beta}=%.2e M_{\odot}$' % Mgas)
    ax.legend(loc=3, fontsize=7)
    fig.canvas.draw()

interact(update, log_rhoCore=(-4, 0, 0.1), beta=(0.0, 1, 0.01), rCore=(1,100,1));

Compare the x-ray flux to the ROSAT detection limit

The x-ray luminosity $L_x$ is calculated until the end of the radius array rr.

Assuming distance to the object is D, the flux is calculated as $F_x = \frac{L_x}{4\pi D^2}$





In [ ]:



    


%matplotlib notebook
# Distance to the object
D = 500*1E3 # 500 Mpc

# ROSAT PSPC flux detection limit in the 0.1-2.4 keV energy band
# From Voges et al., 1999, caption of Fig. 13
# F_limit ~= 0.05 counts s^-1 ~= 5E-13 erg cm^-2 s^-1
F_ROSAT = 5E-13 # erg cm^-2 s^-1

rr = np.logspace(-1, 3.5, 100)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

rhoCore, rCore, beta = (0.5, 2, 0.6)

density = rhoCore*(1.0 + (rr/rCore)**2)**(-1.5*beta)
label = r'$\beta$-model $r_c=%i$kpc, $\rho_c=%.3f$cm$^{-3}$ $\beta$=%.1f' % (rCore, rhoCore, beta)
line, = ax.plot(rr, density, label=label)
Mgas = rhoCore*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta, 2.5, -(rr[-1]/rCore)**2)
text = ax.text(0.6, 0.93, r'$M_{gas,\beta}=%.2e M_{\odot}$' % Mgas, transform=ax.transAxes)

L = Lx(rr[-1], rCore, beta, rhoCore)
Fx = L/4/np.pi/D**2/kpc**2
text2 = ax.text(0.6, 0.85, r'$L_{X}=%.2e erg\cdot s^{-1}$' % L, transform=ax.transAxes)
text3 = ax.text(0.6, 0.77, r'$F_{X,%iMpc}=%.3f F_{X,ROSAT}$' % (D/1E3, Fx/F_ROSAT), transform=ax.transAxes)

ax.errorbar(R_FW2011, ne_FW2011, yerr=ne_err_FW2011, fmt='.', capsize=5, label=r'Freeland & Wilcots 2011 ($n_e$)', c='k')

ax.legend(fontsize=7)
ax.set_xlim(rr[0], rr[-1])
plt.semilogy()
plt.semilogx()
plt.grid(ls=':')
plt.xlabel('radius (kpc)')
plt.ylabel(r'$n_e (cm^{-3})$')

def update(log_rhoCore=-2.5, beta=0.3, rCore=50):
    rhoCore = 10**log_rhoCore
    density = rhoCore*(1.0 + (rr/rCore)**2)**(-1.5*beta)
    line.set_ydata(density)
    label = r'$\beta$-model $r_c=%i$kpc, $n_c=%.3f$cm$^{-3}$ $\beta$=%.1f' % (rCore, rhoCore, beta)
    line.set_label(label)
    
    Mgas = rhoCore*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta, 2.5, -(rr[-1]/rCore)**2)
    text.set_text(r'$M_{gas,\beta}=%.2e M_{\odot}$' % Mgas)
    L = Lx(rr[-1], rCore, beta, rhoCore)
    Fx = L/4/np.pi/D**2/kpc**2
    text2.set_text(r'$L_{X}=%.2e erg\cdot s^{-1}$' % L)
    text3.set_text(r'$F_{X,%iMpc}=%.3f F_{X,ROSAT}$' % (D/1E3, Fx/F_ROSAT))
    ax.legend(loc=3, fontsize=7)
    fig.canvas.draw()

interact(update, log_rhoCore=(-4, 0, 0.1), beta=(0.0, 1, 0.01), rCore=(1,3000,1));



    











In [ ]:



    


%matplotlib inline
plt.rcParams['figure.dpi'] = 150
# cosmic background density: 1 protron per 4 cubic meters
n_cosmicBG = 2.5E-7

#nCores = [ 1E-2, 1E-3, 1E-4, 1E-5]
nCores = [1E-2]

betas = [0.3, 0.5, 0.7]
#betas  		= [0.7]
#rCores		= [100, 30, 10]
rCores = [10]

# Ratio of number of all particles to electrons
e_ratio = 1.93
# Data from Table 1 of Freeland & Wilcots 2011
R = np.array([300, 2000, 270, 700, 15])
n_IGM = np.array([3E-3, 5E-4, 2E-4, 5E-4, 2E-3])/e_ratio
n_IGM_err = np.array([2E-3, 4E-4, 1E-4, 2E-4, 1E-3])/e_ratio

rr = np.logspace(-1, 3, 100)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.semilogy()
#ax2 = ax.twinx()
plt.semilogy()
for nCore in nCores:
    for beta in betas:
        for rCore in rCores:
            density = nCore*(1.0 + (rr/rCore)**2)**(-1.5*beta)
            Mgas = nCore*mH/Msun*kpc**3*1/3*rr[-1]**3*hyp2f1(1.5, 1.5*beta, 2.5, -(rr[-1]/rCore)**2)
            L = Lx(rr[-1], rCore, beta, nCore)
            line, = ax.plot(rr, density, ls='--', 
                     label=r'$r_c=%.1f$kpc, $n_c$=%.3fcm$^{-3}$, $\beta$=%.1f, $M_{gas}=%.1e M_{\odot}, L_{X,Z=0.5}^{0.1-2.4keV}=%.1e erg/s$' %
                            (rCore, nCore, beta, Mgas, L))
            #e_x = 3E-23*rhoCore**2*(1.0 + (rr/rCore)**2)**(-3*beta+0.5)*rCore**kpc/4/np.pi/D/D
            #ax2.plot(rr, e_x)

ax.errorbar(R_FW2011, ne_FW2011, yerr=ne_err_FW2011, fmt='.', capsize=5, 
            label=r'Freeland & Wilcots 2011 (n -> $n_e$)', c='k')
ax.errorbar(R_Pisano2004, ne_Pisano2004, yerr=ne_Pisano2004, fmt='.', lolims=True, 
            label=r'Pisano et al. 2004 ($n_H$ -> $n_e$)', c='grey')
ax.fill_between(rr, n_cosmicBG, 1E-8, color='grey', alpha=0.3)

ax.legend(loc=3, fontsize=7)
ax.set_xlim(rr[0], rr[-1])
ax.set_ylim(1E-8, 2E-1)

plt.semilogx()


plt.grid(ls=':')
plt.xlabel('radius (kpc)')
plt.ylabel(r'$n_e$ (cm$^{-3}$)')



    











In [ ]:



    


# Testing the integration of beta-model 
from scipy.special import hyp2f1
b = 0.8
rmax = 300
rc = 10
print(1/3*rmax**3*hyp2f1(1.5, 1.5*b, 2.5, -(rmax/rc)**2))

from scipy.integrate import quad
quad(lambda r: r**2*(1.0 + (r/rc)**2)**(-1.5*b), 0, rmax)



    











In [ ]:



    


# Data from Table 1 of Freeland & Wilcots 2011
R = np.array([300, 2000, 270, 700, 15])
n_IGM = np.array([3E-3, 5E-4, 2E-4, 5E-4, 2E-3])
n_IGM_err = np.array([2E-3, 4E-4, 1E-4, 2E-4, 1E-3])

plt.errorbar(R, n_IGM, yerr=n_IGM_err, fmt='x', capsize=5)
plt.show()



    











In [ ]:



    


sim_Tcore		= 1.16E7
sim_Tout		= 1.16E7
sim_rhoCore		= 1*1.67E-27		# g/cm3 = N particles/cm3 * 1.67E-24 = 1.91 * N electrons/cm3 * 1.67E-24
sim_rhoFloor	= 1.67E-28
sim_mu			= 0.61			# mean molecular weight
sim_windVel     = 0.0
sim_gammaICM	= 1.66666
sim_bzAmbient	= 0.0 #1.8623E-6
sim_densityProfile = "betacore"
sim_densityBeta = 0.83			# for the beta-model density profile
sim_rCore		= 3.1E22		# 1.26 arcmin *20.81 kpc/arcmin
sim_rCoreT		= 3.1E22		# 60kpc

kpc = 3.086e+21
kB = 1.38E-16
R = 8.3144725E7
mH = 1.6737352238051868e-24 #g

distance = np.logspace(-1, 3, 100)*kpc
e = 1.602176562e-09

densityBG = sim_rhoCore*(1.0 + (distance/sim_rCore)**2)**(-1.5*sim_densityBeta)
densityBG = np.clip(densityBG, 1.67E-28, 1)
tempBG = sim_Tout*(1.0+(distance/sim_rCoreT)**3)/(sim_Tout/sim_Tcore+(distance/sim_rCoreT)**3)
pressureBG = densityBG*tempBG*R/sim_mu
entropyBG = pressureBG*((densityBG/mH*sim_mu)**(-sim_gammaICM))/e



    











In [ ]:



    


%matplotlib inline
#plt.style.use('dark_background')
keV = 1.1604525E7
rho = 1.67E-24
T = 1
P = 1E-10
entropy = 1E5

plt.plot(distance/kpc, densityBG/rho, label='density / $(1.67x10^{-24} g/cm^3)$')
plt.plot(distance/kpc, tempBG/keV/T, ls='--', label='temperature / $keV$')
plt.plot(distance/kpc, pressureBG/P, ls=':', label='pressure / $(10^{-10} dyn / cm^2)$')
#plt.plot(distance/kpc, entropyBG/entropy, ls=':', label='entropy / 1E5 $keV*cm^2$')
plt.semilogy()
plt.semilogx()
plt.grid(ls=':')
plt.xlabel('radius (kpc)')
plt.xlim(0,1500)
plt.legend(loc=1)
plt.title("Galaxy Group Profile")
#plt.savefig("ICM_profile.pdf")



    











In [ ]:

Content source: yihaochen/FLASHtools

Similar notebooks:

notebook.community | gallery | about