Astrophysics — Part VI

Part VI: Extragalactic Astronomy

From galaxy morphology and the luminosity function through active galactic nuclei, galaxy clusters, large-scale structure, and the physics of galaxy formation and evolution in a $\Lambda$CDM universe.

1. Galaxy Classification and Morphology

Galaxies are the fundamental building blocks of the visible universe. They span a remarkable range of sizes (dwarfs with $M_* \sim 10^6\,M_\odot$ to giant ellipticals with$M_* \sim 10^{12}\,M_\odot$), morphologies, colors, and star-formation histories. A systematic classification scheme is essential for understanding the physical processes that govern galaxy evolution.

1.1 The Hubble Sequence

Edwin Hubble’s 1926 classification scheme — the “tuning fork” diagram — remains the foundation of galaxy morphology. It divides galaxies into three broad families:

Ellipticals (E0–E7): Smooth, featureless light distributions classified by apparent ellipticity $\epsilon = 1 - b/a$, where$a$ and $b$ are the semi-major and semi-minor axes. An E0 galaxy appears circular; an E7 is highly elongated. These systems are pressure-supported (random stellar motions dominate over ordered rotation), gas-poor, and dominated by old stellar populations.
Lenticulars (S0): Transitional systems with a prominent bulge and a disk component but no spiral arms. They bridge the gap between ellipticals and spirals.
Spirals (Sa–Sd, SBa–SBd): Disk-dominated systems with spiral arms. The sequence from Sa to Sd reflects decreasing bulge-to-disk ratio, increasingly open spiral arms, and higher gas fractions. The “SB” prefix denotes barred spirals, which possess a linear stellar bar extending from the central bulge. The Milky Way is classified as SBbc.
Irregulars (Irr): Galaxies lacking a clear symmetric structure. They include Irr I (some structure, e.g., the Magellanic Clouds) and Irr II (truly chaotic morphologies, often the result of mergers or strong tidal interactions).

1.2 The Sérsic Surface Brightness Profile

A quantitative description of galaxy morphology uses the Sérsic (1963) surface brightness profile, which generalizes both the de Vaucouleurs law (for ellipticals) and the exponential disk profile (for spirals) into a single parametric form:

$$I(r) = I_e \exp\!\left[-b_n\left(\left(\frac{r}{r_e}\right)^{1/n} - 1\right)\right]$$

where $I_e$ is the surface brightness at the effective (half-light) radius $r_e$,$n$ is the Sérsic index, and $b_n \approx 2n - 1/3 + 4/(405n)$ ensures that $r_e$ encloses half the total luminosity.

de Vaucouleurs Profile (n = 4)

The $r^{1/4}$ law, empirically established by de Vaucouleurs (1948), describes elliptical galaxies and classical bulges. Setting $n = 4$ gives $b_4 \approx 7.67$:

$$I(r) = I_e \exp\!\left[-7.67\left(\left(\frac{r}{r_e}\right)^{1/4} - 1\right)\right]$$

Exponential Disk (n = 1)

The exponential profile describes galactic disks. With $n = 1$ and $b_1 \approx 1.68$, this reduces to $I(r) = I_0 \exp(-r/h)$ where $h = r_e/1.68$ is the disk scale length.

$$I(r) = I_0\,e^{-r/h}, \qquad h \equiv r_e / b_1$$

The total luminosity integrated over the full Sérsic profile is:

$$L_{\text{tot}} = 2\pi\,I_e\,r_e^2\,\frac{n\,e^{b_n}}{b_n^{2n}}\,\Gamma(2n)$$

where $\Gamma(2n)$ is the complete gamma function. Dwarf galaxies typically have$n \sim 0.5\text{--}1$, disk galaxies $n \sim 1$, and giant ellipticals$n \sim 4\text{--}10$. The Sérsic index correlates strongly with galaxy mass, luminosity, and stellar velocity dispersion.

1.3 Galaxy Color Bimodality

One of the most striking features of the galaxy population is its bimodal color distribution. When galaxies are plotted in a color–magnitude diagram (e.g., $u - r$ versus$M_r$), they separate cleanly into two populations:

Red Sequence

Massive, gas-poor, quiescent galaxies dominated by old stellar populations. Primarily ellipticals and S0s. The red sequence shows a tight color–magnitude relation (redder at higher luminosity) driven by the mass–metallicity relation. Typical colors: $u - r \gtrsim 2.2$.

Green Valley

A sparsely populated transition zone between the blue cloud and red sequence. Galaxies here are thought to be undergoing quenching — their star formation is being shut down by AGN feedback, gas exhaustion, or environmental processes. Residence time in the green valley is short ($\sim 1$–2 Gyr).

Blue Cloud

Actively star-forming galaxies with young, hot stellar populations. Predominantly spirals and irregulars. Broader distribution of colors reflecting a range of star formation rates and dust contents. Typical colors: $u - r \lesssim 2.0$.

1.4 The Fundamental Plane of Elliptical Galaxies

Elliptical galaxies obey a tight three-parameter scaling relation known as the Fundamental Plane(Djorgovski & Davis 1987; Dressler et al. 1987). In the space of effective radius $r_e$, central velocity dispersion $\sigma_0$, and mean surface brightness within $r_e$(denoted $\langle I \rangle_e$ or equivalently $\langle \mu \rangle_e$), ellipticals lie on a plane:

$$r_e \propto \sigma_0^{1.24}\,\langle I \rangle_e^{-0.82}$$

Observed Fundamental Plane relation in the r-band (Bernardi et al. 2003).

If ellipticals were perfectly homologous systems in virial equilibrium (with constant mass-to-light ratio), the virial theorem would predict:

$$r_e \propto \sigma_0^{2}\,\langle I \rangle_e^{-1} \qquad \text{(virial prediction)}$$

The departure from the virial prediction (the “tilt” of the Fundamental Plane) implies systematic variations in $M/L$ with mass: more massive ellipticals have higher $M/L$.

The Fundamental Plane has a remarkably small scatter ($\sim 15\%$ in $r_e$), making it a powerful distance indicator. Related two-parameter projections include the Faber–Jackson relation ($L \propto \sigma^4$) and the$D_n\text{--}\sigma$ relation.

2. The Galaxy Luminosity Function

The galaxy luminosity function $\Phi(L)$ describes the number density of galaxies as a function of their luminosity. It is one of the most fundamental statistical descriptions of the galaxy population, encoding information about galaxy formation efficiency, feedback processes, and the underlying dark matter halo mass function.

2.1 The Schechter Function

Schechter (1976) proposed a simple analytic form that provides an excellent fit to the observed luminosity function over a wide range of luminosities:

$$\Phi(L)\,dL = \phi^*\left(\frac{L}{L^*}\right)^{\alpha}\exp\!\left(-\frac{L}{L^*}\right)\frac{dL}{L^*}$$

Three parameters: $\phi^*$ (normalization, in Mpc$^{-3}$),$L^*$ (characteristic luminosity), and $\alpha$ (faint-end slope).

The physical interpretation of each parameter:

$\phi^*$: Overall normalization — the number density of galaxies at $L \sim L^*$. Typical value:$\phi^* \sim 1.2 \times 10^{-2}\,\text{Mpc}^{-3}$.
$L^*$: The characteristic (or “knee”) luminosity. For $L \gg L^*$, the exponential cutoff dominates, producing a sharp decline in the abundance of very luminous galaxies. In the r-band,$M^* \approx -20.83$ (corresponding to $L^* \approx 2 \times 10^{10}\,L_\odot$).
$\alpha$: The faint-end slope, governing the abundance of dwarf galaxies. For $\alpha < -1$, the number density diverges as$L \to 0$ (though the total luminosity density remains finite for $\alpha > -2$). Observed values range from $\alpha \approx -1.0$ to $-1.6$ depending on environment and band.

2.2 Magnitude Form

Since astronomical observations are typically reported in magnitudes, we often write the luminosity function in terms of absolute magnitude $M$. Using $L/L^* = 10^{0.4(M^*-M)}$ and$dL/L^* = 0.4\ln(10)\cdot 10^{0.4(M^*-M)}\,dM$, the Schechter function becomes:

$$\Phi(M)\,dM = 0.4\ln(10)\,\phi^*\,\left[10^{0.4(M^*-M)}\right]^{\alpha+1}\exp\!\left(-10^{0.4(M^*-M)}\right)\,dM$$

2.3 Integrated Quantities

The total number density of galaxies (integrating from $L = 0$ to $\infty$):

$$n_{\text{tot}} = \int_0^\infty \Phi(L)\,dL = \phi^*\,\Gamma(\alpha + 1)$$

This integral converges only for $\alpha > -1$; for steeper faint-end slopes, one must impose a minimum luminosity cutoff.

The total luminosity density (luminosity per unit comoving volume):

$$j = \int_0^\infty L\,\Phi(L)\,dL = \phi^*\,L^*\,\Gamma(\alpha + 2)$$

This converges for $\alpha > -2$. The observed r-band luminosity density is$j \approx 1.8 \times 10^8\,h\,L_\odot\,\text{Mpc}^{-3}$.

The luminosity density can be converted to a stellar mass density using a mass-to-light ratio ($M/L$) appropriate for the stellar population. Current estimates give a cosmic stellar mass density of $\rho_* \approx (4\text{--}5) \times 10^8\,M_\odot\,\text{Mpc}^{-3}$, corresponding to a baryon-to-star conversion efficiency of roughly 5–7%.

Simulation: Schechter Luminosity Function

The following Python code plots the Schechter galaxy luminosity function for several parameter sets, illustrating the effects of varying the faint-end slope $\alpha$, along with the cumulative number counts and luminosity density integrand.

Python

schechter_luminosity_function.py150 lines

#!/usr/bin/env python3
"""
schechter_luminosity_function.py
Compute and plot the Schechter galaxy luminosity function for
different parameter sets, along with cumulative number counts.
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gammaincc, gamma as gamma_func
from scipy.integrate import quad

# ------------------------------------------------------------------
# Schechter function in luminosity space
# Phi(L) dL = phi_star (L/L_star)^alpha exp(-L/L_star) d(L/L_star)
# ------------------------------------------------------------------
def schechter_L(L, phi_star, L_star, alpha):
    """Schechter function in luminosity units (per Mpc^3 per L_star)."""
    x = L / L_star
    return phi_star * x**alpha * np.exp(-x) / L_star

# ------------------------------------------------------------------
# Schechter function in absolute magnitude space
# Phi(M) = 0.4 ln(10) phi_star [10^{0.4(M*-M)}]^(alpha+1)
#           * exp(-10^{0.4(M*-M)})
# ------------------------------------------------------------------
def schechter_M(M, phi_star, M_star, alpha):
    """Schechter function per magnitude per Mpc^3."""
    x = 10.0**(0.4 * (M_star - M))
    return 0.4 * np.log(10) * phi_star * x**(alpha + 1) * np.exp(-x)

# ------------------------------------------------------------------
# Cumulative number density  n(>L) = phi_star * Gamma(alpha+1, L/L_star)
# using the upper incomplete gamma function
# ------------------------------------------------------------------
def cumulative_n(L_over_Lstar, phi_star, alpha):
    """Cumulative number density of galaxies brighter than L."""
    # scipy gammaincc = regularized upper incomplete gamma
    # Gamma(a, x) = gamma_func(a) * gammaincc(a, x)
    a = alpha + 1.0
    return phi_star * gamma_func(a) * gammaincc(a, L_over_Lstar)

# ------------------------------------------------------------------
# Parameter sets (typical values from surveys)
# ------------------------------------------------------------------
params = [
    {"phi_star": 1.2e-2, "L_star": 1.0, "M_star": -20.83, "alpha": -1.07,
     "label": r"SDSS r-band ($\alpha=-1.07$)", "color": "#818cf8", "ls": "-"},
    {"phi_star": 1.6e-2, "L_star": 1.0, "M_star": -20.04, "alpha": -1.25,
     "label": r"Blue galaxies ($\alpha=-1.25$)", "color": "#38bdf8", "ls": "--"},
    {"phi_star": 0.5e-2, "L_star": 1.0, "M_star": -21.2,  "alpha": -0.5,
     "label": r"Red sequence ($\alpha=-0.5$)", "color": "#f87171", "ls": "-."},
    {"phi_star": 1.0e-2, "L_star": 1.0, "M_star": -20.5,  "alpha": -1.6,
     "label": r"Steep faint-end ($\alpha=-1.6$)", "color": "#4ade80", "ls": ":"},
]

# ------------------------------------------------------------------
# Figure 1: Schechter function in magnitude space
# ------------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(13, 10))
fig.patch.set_facecolor('#0f172a')
for ax in axes.flat:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='white')
    for spine in ax.spines.values():
        spine.set_color('#334155')

M_arr = np.linspace(-24, -14, 500)

ax1 = axes[0, 0]
for p in params:
    phi = schechter_M(M_arr, p["phi_star"], p["M_star"], p["alpha"])
    ax1.semilogy(M_arr, phi, label=p["label"], color=p["color"],
                 ls=p["ls"], lw=2)
ax1.set_xlabel("Absolute Magnitude M", fontsize=12, color='white')
ax1.set_ylabel(r"$\Phi(M)$ [Mpc$^{-3}$ mag$^{-1}$]", fontsize=12, color='white')
ax1.set_title("Schechter Luminosity Function", fontsize=14, color='white', fontweight='bold')
ax1.set_ylim(1e-7, 1e-1)
ax1.set_xlim(-24, -14)
ax1.invert_xaxis()
ax1.legend(fontsize=9, facecolor='#1e293b', edgecolor='#475569', labelcolor='white')
ax1.grid(alpha=0.2)

# ------------------------------------------------------------------
# Figure 2: In luminosity space (log-log)
# ------------------------------------------------------------------
ax2 = axes[0, 1]
L_arr = np.logspace(-3, 2, 500)  # in units of L_star
for p in params:
    phi = schechter_L(L_arr, p["phi_star"], 1.0, p["alpha"])
    ax2.loglog(L_arr, phi * 1.0, color=p["color"], ls=p["ls"], lw=2)
ax2.set_xlabel(r"$L / L^*$", fontsize=12, color='white')
ax2.set_ylabel(r"$\Phi(L)\;L^*$ [Mpc$^{-3}$]", fontsize=12, color='white')
ax2.set_title("Luminosity Space (log-log)", fontsize=14, color='white', fontweight='bold')
ax2.set_xlim(1e-3, 30)
ax2.set_ylim(1e-8, 1e0)
ax2.grid(alpha=0.2)

# ------------------------------------------------------------------
# Figure 3: Cumulative number density n(>L)
# ------------------------------------------------------------------
ax3 = axes[1, 0]
L_frac = np.logspace(-2, 1.5, 300)
for p in params:
    n_cum = np.array([cumulative_n(ll, p["phi_star"], p["alpha"]) for ll in L_frac])
    ax3.loglog(L_frac, n_cum, color=p["color"], ls=p["ls"], lw=2)
ax3.axvline(1.0, color='white', ls=':', alpha=0.4, lw=1)
ax3.text(1.1, 2e-2, r"$L = L^*$", color='white', fontsize=10)
ax3.set_xlabel(r"$L_{\min}/L^*$", fontsize=12, color='white')
ax3.set_ylabel(r"$n(>L)$ [Mpc$^{-3}$]", fontsize=12, color='white')
ax3.set_title("Cumulative Number Density", fontsize=14, color='white', fontweight='bold')
ax3.grid(alpha=0.2)

# ------------------------------------------------------------------
# Figure 4: Luminosity density integrand  L * Phi(L)
# ------------------------------------------------------------------
ax4 = axes[1, 1]
for p in params:
    phi = schechter_L(L_arr, p["phi_star"], 1.0, p["alpha"])
    lum_density = L_arr * phi * 1.0  # L * Phi(L) * L_star
    ax4.semilogx(L_arr, lum_density, color=p["color"], ls=p["ls"], lw=2)
ax4.set_xlabel(r"$L / L^*$", fontsize=12, color='white')
ax4.set_ylabel(r"$L\,\Phi(L)$ [Mpc$^{-3}$]", fontsize=12, color='white')
ax4.set_title("Luminosity Density Integrand", fontsize=14, color='white', fontweight='bold')
ax4.set_xlim(1e-3, 30)
ax4.grid(alpha=0.2)

plt.tight_layout(pad=2.0)
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.show()
print("Schechter luminosity function plots saved to output.png")

# ------------------------------------------------------------------
# Print summary table
# ------------------------------------------------------------------
print("\n" + "="*65)
print("Schechter Function Parameter Summary")
print("="*65)
for p in params:
    a = p["alpha"]
    # Total number density = phi_star * Gamma(alpha+1)  [for alpha > -1]
    if a + 1 > 0:
        n_tot = p["phi_star"] * gamma_func(a + 1)
        print(f"  {p['label']:40s}  n_tot = {n_tot:.4e} Mpc^-3")
    else:
        print(f"  {p['label']:40s}  n_tot diverges (alpha <= -1)")
    # Luminosity density = phi_star * L_star * Gamma(alpha+2)
    lum_dens = p["phi_star"] * 1.0 * gamma_func(a + 2)
    print(f"  {'':40s}  j = {lum_dens:.4e} L* Mpc^-3")
print("="*65)

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation Notes

The upper-left panel shows $\Phi(M)$ for four parameter sets: the SDSS r-band field luminosity function ($\alpha = -1.07$), blue star-forming galaxies ($\alpha = -1.25$), the red sequence ($\alpha = -0.5$, shallow), and an extreme steep faint end ($\alpha = -1.6$). The cumulative panel (lower-left) shows how the total number density depends strongly on $\alpha$at low luminosities. The luminosity density integrand (lower-right) peaks near $L \sim L^*$ for all cases, showing that most of the light comes from galaxies near the characteristic luminosity.

3. Active Galactic Nuclei

Active Galactic Nuclei (AGN) are compact, extremely luminous regions at the centers of galaxies powered by accretion onto supermassive black holes (SMBHs). AGN emit across the entire electromagnetic spectrum and can outshine their entire host galaxy by factors of 100 or more. They represent some of the most energetic persistent phenomena in the universe, with luminosities reaching$L \sim 10^{46}\text{--}10^{48}\,\text{erg\,s}^{-1}$.

3.1 The Unified Model of AGN

The AGN unified model (Antonucci 1993; Urry & Padovani 1995) posits that the apparent diversity of AGN types arises primarily from orientation effects — different viewing angles to a common underlying structure consisting of:

Central SMBH: Mass $M_{\bullet} \sim 10^6\text{--}10^{10}\,M_\odot$. The gravitational engine powering the AGN. The Schwarzschild radius is$r_S = 2GM_\bullet/c^2 \approx 3 \times 10^{13}\,(M_\bullet/10^8\,M_\odot)\,\text{cm}$.
Accretion Disk: A geometrically thin, optically thick disk (Shakura & Sunyaev 1973) extending from the innermost stable circular orbit (ISCO) to$\sim 10^3\,r_S$. The disk temperature profile follows$T(r) \propto r^{-3/4}$, peaking at UV/soft X-ray wavelengths. The accretion luminosity is:
$$L_{\text{acc}} = \eta\,\dot{M}\,c^2$$
where the radiative efficiency $\eta \approx 0.06\text{--}0.42$ depends on BH spin.
Broad-Line Region (BLR): Dense ($n_e \sim 10^{9}\text{--}10^{11}\,\text{cm}^{-3}$), fast-moving gas clouds at $r \sim 0.01\text{--}0.1\,\text{pc}$ from the SMBH, producing broad emission lines ($\text{FWHM} \sim 2{,}000\text{--}10{,}000\,\text{km\,s}^{-1}$) through photoionization by the accretion disk continuum.
Dusty Torus: A geometrically thick, axially symmetric structure of dust and gas at $r \sim 0.1\text{--}10\,\text{pc}$, obscuring the central engine and BLR for observers at high inclination angles.
Narrow-Line Region (NLR): Lower-density gas ($n_e \sim 10^{3}\text{--}10^{6}\,\text{cm}^{-3}$) at $r \sim 10\text{--}1{,}000\,\text{pc}$, producing narrow emission lines ($\text{FWHM} \sim 200\text{--}900\,\text{km\,s}^{-1}$). The NLR is visible from all viewing angles.
Relativistic Jet: Collimated outflows of plasma at relativistic speeds ($\Gamma \sim 5\text{--}50$), launched along the BH spin axis. Only $\sim 10\%$ of AGN are radio-loud with powerful jets.

3.2 AGN Classification by Orientation

Face-on (Type 1)

When viewed along the polar axis (face-on to the torus opening), both the BLR and NLR are visible. The spectrum shows broad and narrow emission lines superimposed on a strong featureless continuum.

• Seyfert 1: Low-luminosity AGN with broad + narrow lines
• Quasar (QSO): High-luminosity counterpart, $M_B < -23$
• Blazar (BL Lac): Jet pointed directly at observer; featureless continuum

Edge-on (Type 2)

When viewed through the dusty torus (edge-on), the BLR is obscured. Only narrow emission lines are visible in the optical spectrum.

• Seyfert 2: Narrow lines only; BLR hidden by torus
• Narrow-line radio galaxy (NLRG): Radio-loud Type 2
• Type 2 QSO: Luminous but obscured quasar

3.3 Quasars

Quasars (quasi-stellar objects) are the most luminous AGN, with bolometric luminosities reaching$L_{\text{bol}} \sim 10^{46}\text{--}10^{48}\,\text{erg\,s}^{-1}$. They were first identified in 1963 by Maarten Schmidt, who recognized that the puzzling emission lines of 3C 273 were familiar hydrogen lines redshifted to $z = 0.158$. The implied luminosity was$\sim 4 \times 10^{46}\,\text{erg\,s}^{-1}$ — hundreds of times the luminosity of the entire Milky Way, emanating from a region smaller than the Solar System.

The Eddington luminosity sets an upper limit on the luminosity for spherically symmetric accretion:

$$L_{\text{Edd}} = \frac{4\pi\,G\,M_\bullet\,m_p\,c}{\sigma_T} \approx 1.26 \times 10^{38}\left(\frac{M_\bullet}{M_\odot}\right)\,\text{erg\,s}^{-1}$$

A quasar at $L \sim 10^{46}\,\text{erg\,s}^{-1}$ requires$M_\bullet \gtrsim 10^8\,M_\odot$ if accreting at the Eddington rate.

3.4 The M–sigma Relation

One of the most important discoveries in extragalactic astronomy is the tight correlation between the mass of the central SMBH and the stellar velocity dispersion $\sigma$ of the host galaxy’s bulge (Ferrarese & Merritt 2000; Gebhardt et al. 2000):

$$M_\bullet = (1.9 \pm 0.6) \times 10^8\,M_\odot \left(\frac{\sigma}{200\,\text{km\,s}^{-1}}\right)^{4.24 \pm 0.41}$$

The M–$\sigma$ relation (Gültekin et al. 2009), with intrinsic scatter$\sim 0.3$ dex. The exponent is approximately 4–5 depending on the sample.

This remarkably tight scaling relation — connecting a black hole with$r_S \sim 10^{-5}\,\text{pc}$ to the kinematics of stars at$\sim 1\text{--}10\,\text{kpc}$ — implies a fundamental co-evolution of SMBHs and their host galaxies, mediated by AGN feedback.

3.5 AGN Feedback

AGN feedback refers to the energy and momentum injection from the active nucleus into the surrounding gas, regulating (and ultimately quenching) star formation in the host galaxy. There are two principal modes:

Quasar (Radiative) Mode

High-accretion-rate AGN ($\dot{M}/\dot{M}_{\text{Edd}} \gtrsim 0.01$) drive powerful winds and outflows through radiation pressure on dust and line-driven forces. These winds can reach velocities of $v \sim 0.1c$ (ultra-fast outflows, UFOs) and sweep gas out of the galaxy, suppressing star formation on galaxy-wide scales. Energy coupling:$\dot{E}_{\text{wind}} \sim 0.05\,L_{\text{bol}}$.

Radio (Kinetic) Mode

Low-accretion-rate AGN ($\dot{M}/\dot{M}_{\text{Edd}} \lesssim 0.01$) in massive ellipticals inject energy via relativistic jets that inflate cavities (bubbles) in the hot ICM. This “maintenance mode” feedback prevents the hot gas from cooling and fueling new star formation. Jet power: $P_{\text{jet}} \sim 10^{43}\text{--}10^{46}\,\text{erg\,s}^{-1}$. X-ray observations of cavities in clusters (e.g., Perseus, MS 0735) directly measure the mechanical work done by the jet.

3.6 Reverberation Mapping

Reverberation mapping is the premier technique for measuring BLR sizes and SMBH masses in AGN. The BLR gas is photoionized by the variable continuum emission from the accretion disk. Changes in the continuum luminosity are echoed in the broad emission lines after a time delay $\tau$ equal to the light-travel time across the BLR:

$$R_{\text{BLR}} = c\,\tau$$

Combined with the BLR velocity dispersion from the line width:

$$M_\bullet = f\,\frac{c\,\tau\,(\Delta V)^2}{G}$$

where $\Delta V$ is the velocity width of the broad line (FWHM or$\sigma_{\text{line}}$) and $f$ is a dimensionless factor of order unity encoding the BLR geometry and kinematics (typically $f \approx 4\text{--}5$ when using$\sigma_{\text{line}}$).

An empirical scaling discovered from reverberation mapping relates the BLR radius to the AGN luminosity:

$$R_{\text{BLR}} \approx 33\,\text{lt-days}\left(\frac{\lambda L_\lambda(5100\,\text{\AA})}{10^{44}\,\text{erg\,s}^{-1}}\right)^{0.53}$$

This $R\text{--}L$ relation enables single-epoch (or “virial”) BH mass estimates for large quasar samples using just a spectrum, enabling demographic studies of SMBH growth across cosmic time.

4. Galaxy Clusters

Galaxy clusters are the most massive gravitationally bound structures in the universe, with total masses of $M \sim 10^{14}\text{--}10^{15}\,M_\odot$. They contain hundreds to thousands of galaxies, but the galaxies themselves represent only $\sim 2\text{--}5\%$ of the total mass. The dominant baryonic component ($\sim 12\text{--}15\%$) is the hot intracluster medium (ICM), and the remainder ($\sim 80\text{--}85\%$) is dark matter.

4.1 Virial Mass

The simplest mass estimate for a cluster uses the virial theorem applied to the galaxy velocities. For a self-gravitating system in virial equilibrium ($2K + W = 0$):

$$M_{\text{vir}} = \frac{3\,\sigma_v^2\,R_v}{G}$$

where $\sigma_v$ is the line-of-sight velocity dispersion of the member galaxies and$R_v$ is the virial radius (related to the projected harmonic mean radius). For rich clusters, $\sigma_v \sim 800\text{--}1{,}500\,\text{km\,s}^{-1}$ and$R_v \sim 1\text{--}3\,\text{Mpc}$.

Derivation: The virial theorem for a self-gravitating system states$2\langle K \rangle + \langle W \rangle = 0$. The total kinetic energy is$K = \frac{1}{2}M\langle v^2 \rangle = \frac{3}{2}M\sigma_v^2$ (assuming isotropic velocities so $\langle v^2 \rangle = 3\sigma_v^2$). The gravitational potential energy for a uniform sphere is $W = -3GM^2/(5R_v)$, but the factor $3/5$ is replaced by a profile-dependent constant. The standard convention absorbs this into the definition of $R_v$.

4.2 The Intracluster Medium (ICM)

The ICM is a hot, tenuous plasma filling the gravitational potential well of the cluster. It is heated to the virial temperature by gravitational infall and shocks:

$$k_B T_{\text{vir}} \sim \frac{G\,M\,\mu\,m_p}{2\,R_v} \sim 5\text{--}15\,\text{keV}$$

corresponding to temperatures $T \sim 10^7\text{--}10^8\,\text{K}$.

At these temperatures, the ICM is fully ionized and emits thermal bremsstrahlung (free-free emission) in the X-ray band. The X-ray luminosity scales as:

$$L_X \propto \int n_e^2\,\Lambda(T)\,dV \propto n_e^2\,T^{1/2}\,V$$

where $\Lambda(T) \propto T^{1/2}$ is the bremsstrahlung cooling function and$n_e \sim 10^{-3}\text{--}10^{-1}\,\text{cm}^{-3}$ is the electron density.

4.3 Hydrostatic Mass Estimation

If the ICM is in hydrostatic equilibrium within the cluster gravitational potential, the pressure gradient balances gravity:

$$\frac{dP}{dr} = -\rho_g\,\frac{G\,M(r)}{r^2}$$

For an ideal gas with $P = n\,k_B\,T = \rho_g\,k_B\,T / (\mu\,m_p)$, solving for the enclosed mass gives the hydrostatic mass equation:

$$M_{\text{HSE}}(r) = -\frac{k_B\,T\,r}{G\,\mu\,m_p}\left(\frac{d\ln\rho_g}{d\ln r} + \frac{d\ln T}{d\ln r}\right)$$

where $\mu \approx 0.6$ is the mean molecular weight of the fully ionized ICM. X-ray observations provide both $\rho_g(r)$ (from the surface brightness profile) and $T(r)$ (from spatially resolved spectroscopy).

Hydrostatic mass bias: Simulations consistently show that hydrostatic masses underestimate the true mass by $\sim 10\text{--}30\%$ because non-thermal pressure support from turbulence, bulk motions, magnetic fields, and cosmic rays violates the assumption of pure hydrostatic equilibrium. This bias, parametrized as $b \equiv 1 - M_{\text{HSE}}/M_{\text{true}}$, is a critical systematic uncertainty in cluster cosmology.

4.4 The Sunyaev–Zel’dovich Effect

The Sunyaev–Zel’dovich (SZ) effect is the inverse Compton scattering of cosmic microwave background (CMB) photons by the hot electrons of the ICM. As CMB photons traverse the cluster, they gain energy from the hot electrons, distorting the CMB blackbody spectrum. The thermal SZ effect produces a characteristic decrement at frequencies below $\sim 218\,\text{GHz}$ and an increment above.

$$\frac{\Delta T}{T_{\text{CMB}}} = f(x)\,y$$

where $x = h\nu/(k_B T_{\text{CMB}})$ and the Comptonization parameter is:

$$y = \int \frac{k_B\,T_e}{m_e\,c^2}\,n_e\,\sigma_T\,dl$$

In the Rayleigh–Jeans limit ($x \ll 1$), $f(x) \approx -2$, giving $\Delta T/T_{\text{CMB}} \approx -2y$. The integrated SZ signal$Y = \int y\,d\Omega$ is proportional to the total thermal energy of the ICM and serves as a nearly redshift-independent mass proxy.

The redshift independence of the SZ effect (because CMB photons at higher redshift are hotter, compensating the $(1+z)^4$ surface brightness dimming) makes SZ surveys uniquely powerful for detecting clusters at high redshift. Surveys by the South Pole Telescope (SPT), Atacama Cosmology Telescope (ACT), and the Planck satellite have discovered thousands of clusters via the SZ effect.

4.5 Gravitational Lensing by Clusters

The enormous mass of galaxy clusters deflects light from background sources, providing a direct measure of the total (dark + baryonic) projected mass distribution. Two regimes are distinguished:

Strong Lensing

Near the cluster core ($r \lesssim 100\text{--}300\,\text{kpc}$), the surface mass density exceeds the critical value $\Sigma_{\text{cr}}$, producing multiple images, giant arcs, and Einstein rings. The Einstein radius is:

$$\theta_E = \sqrt{\frac{4GM(\theta_E)}{c^2}\frac{D_{LS}}{D_L\,D_S}}$$

Typical Einstein radii for massive clusters: $\theta_E \sim 20\text{--}50''$.

Weak Lensing

At larger radii ($r \gtrsim 300\,\text{kpc}$), the lensing signal is weak: background galaxies are slightly sheared (distorted tangentially around the cluster). The shear$\gamma(\theta)$ is related to the projected mass distribution:

$$\gamma_t(r) = \frac{\bar{\Sigma}(<r) - \Sigma(r)}{\Sigma_{\text{cr}}} = \frac{\Delta\Sigma(r)}{\Sigma_{\text{cr}}}$$

Stacking many background galaxies provides a statistical mass measurement out to several Mpc.

4.6 Mass–Observable Relations and Cluster Cosmology

Galaxy clusters are powerful cosmological probes because their abundance as a function of mass and redshift is sensitive to both the growth rate of structure and the expansion history. The cluster mass function constrains $\Omega_m$ and $\sigma_8$ (the amplitude of matter fluctuations on 8 $h^{-1}$ Mpc scales). Key scaling relations link the observable to mass:

X-ray luminosity–mass:

$$L_X \propto M^{4/3}\,E(z)^{7/3}$$

Temperature–mass:

$$k_B T \propto M^{2/3}\,E(z)^{2/3}$$

SZ signal–mass:

$$Y_{\text{SZ}} \propto M^{5/3}\,E(z)^{2/3}$$

where $E(z) = H(z)/H_0$. These self-similar scaling relations follow from simple virial arguments assuming gravitational heating only.

5. Large-Scale Structure and the Cosmic Web

On scales of tens to hundreds of megaparsecs, matter in the universe is organized into a complex network known as the cosmic web: a pattern of filaments, walls (sheets), nodes (clusters), and voids that emerged from the gravitational amplification of tiny primordial density fluctuations seeded during cosmic inflation.

5.1 The Matter Power Spectrum

The statistical properties of the density field $\delta(\mathbf{x}) = (\rho - \bar{\rho})/\bar{\rho}$are characterized by the power spectrum $P(k)$, defined through the Fourier transform of the two-point correlation function:

$$\langle \tilde{\delta}(\mathbf{k})\,\tilde{\delta}^*(\mathbf{k}')\rangle = (2\pi)^3\,P(k)\,\delta_D(\mathbf{k} - \mathbf{k}')$$

In configuration space, the equivalent statistic is the two-point correlation function:

$$\xi(r) = \langle \delta(\mathbf{x})\,\delta(\mathbf{x} + \mathbf{r}) \rangle = \int_0^\infty \frac{k^2\,P(k)}{2\pi^2}\,\frac{\sin(kr)}{kr}\,dk$$

The primordial power spectrum is nearly scale-invariant ($P_{\text{prim}}(k) \propto k^{n_s}$with $n_s \approx 0.965$). It is modified by the transfer function $T(k)$, which encodes the effects of radiation pressure, Silk damping, and the matter–radiation equality epoch: $P(k) = A\,k^{n_s}\,T^2(k)\,D^2(a)$, where $D(a)$ is the linear growth factor.

5.2 Baryon Acoustic Oscillations (BAO)

Before recombination ($z \approx 1{,}100$), photons and baryons were tightly coupled into a photon–baryon fluid. Gravity drove acoustic (sound) waves through this fluid, with the sound speed$c_s \approx c/\sqrt{3}$. At recombination, the photons decoupled and the sound waves froze, leaving a characteristic scale imprinted in the matter distribution:

$$r_s = \int_0^{t_{\text{rec}}} \frac{c_s(t)}{a(t)}\,dt \approx 147\,\text{Mpc}$$

The comoving sound horizon at recombination — the BAO standard ruler.

This $\sim 150\,\text{Mpc}$ scale appears as a subtle bump in the galaxy correlation function $\xi(r)$ and as oscillatory wiggles in $P(k)$. By measuring the BAO scale at different redshifts (both in the transverse direction and along the line of sight), we obtain geometric constraints on the angular diameter distance $D_A(z)$ and the Hubble parameter $H(z)$:

$$\Delta\theta = \frac{r_s}{(1+z)\,D_A(z)} \qquad \text{(transverse)}, \qquad \Delta z = \frac{r_s\,H(z)}{c} \qquad \text{(radial)}$$

BAO measurements from galaxy surveys (SDSS, BOSS, DESI) provide some of the most precise constraints on the dark energy equation of state and the Hubble constant, independently of the CMB.

5.3 Cosmic Web Morphology

The cosmic web consists of four distinct environments, classified by the eigenvalues of the tidal tensor (or Hessian of the gravitational potential):

Nodes (Clusters)

Dense knots where filaments intersect. Three positive eigenvalues of the tidal tensor (collapse along all three axes). Contain $\sim 5\%$ of mass,$\lesssim 1\%$ of volume.

Filaments

Elongated structures connecting clusters. Two positive eigenvalues (collapse along two axes, expansion along one). Contain $\sim 40\%$ of mass, $\sim 5\text{--}10\%$of volume. Galaxies flow along filaments toward clusters.

Walls (Sheets)

Planar structures forming the boundaries of voids. One positive eigenvalue (collapse along one axis). Contain $\sim 10\text{--}15\%$ of mass. The Great Wall (Geller & Huchra 1989) and the Sloan Great Wall are famous examples.

Voids

Large, underdense regions with $\delta \lesssim -0.8$. Three negative eigenvalues (expansion along all axes). Contain $\lesssim 10\%$ of mass but$\sim 80\%$ of volume. Typical diameters: 20–50 Mpc.

5.4 The Press–Schechter Mass Function

The Press–Schechter formalism (1974) provides an analytic prediction for the comoving number density of dark matter halos as a function of mass. It is based on the spherical collapse model and the assumption that regions where the smoothed density field exceeds a critical threshold $\delta_c \approx 1.686$will collapse to form bound halos:

$$n(M)\,dM = \sqrt{\frac{2}{\pi}}\,\frac{\bar{\rho}}{M^2}\,\frac{\delta_c}{\sigma(M)}\left|\frac{d\ln\sigma}{d\ln M}\right|\exp\!\left(-\frac{\delta_c^2}{2\sigma^2(M)}\right)\,dM$$

where $\sigma(M)$ is the rms variance of the density field smoothed on a scale$R = (3M/4\pi\bar{\rho})^{1/3}$, and $\delta_c = 1.686$ for an Einstein–de Sitter universe.

Derivation sketch: The fraction of mass in halos above mass $M$ is identified with the probability that the smoothed density field $\delta_R$ (assumed Gaussian) exceeds $\delta_c$:

$$F(>M) = \frac{1}{2}\,\text{erfc}\!\left(\frac{\delta_c}{\sqrt{2}\,\sigma(M)}\right)$$

The famous “fudge factor of 2” (multiplying by 2 to account for underdense regions that are eventually accreted) was later justified rigorously by the excursion set theory of Bond et al. (1991).

Modern calibrations from N-body simulations (e.g., Tinker et al. 2008; Despali et al. 2016) refine the functional form, but the Press–Schechter prediction captures the essential exponential suppression of massive halos and the power-law rise at low masses.

5.5 N-body Simulations

Numerical N-body simulations are indispensable for studying the nonlinear evolution of cosmic structure. They follow the gravitational dynamics of $N$ dark matter particles (typically$10^9\text{--}10^{12}$) in an expanding universe. Landmark simulations include:

Millennium Simulation (2005):$10^{10}$ particles in a $500\,h^{-1}\,\text{Mpc}$ box. Pioneered semi-analytic galaxy formation models applied to merger trees.
Bolshoi (2011):$8.6 \times 10^9$ particles in a $250\,h^{-1}\,\text{Mpc}$ box with$1\,h^{-1}\,\text{kpc}$ force resolution.
IllustrisTNG (2018): A suite of cosmological magneto-hydrodynamic simulations with baryonic physics (star formation, supernova and AGN feedback, magnetic fields). TNG300 spans $300\,\text{Mpc}$ with$2 \times 2500^3$ resolution elements.
AbacusSummit (2021): Over 150 simulations with$6912^3$ particles each, designed for DESI survey analyses.

6. Galaxy Formation and Evolution

In the standard $\Lambda$CDM cosmological model, galaxies form hierarchically: small dark matter halos collapse first, then merge to build up larger structures. Baryons fall into these halos, cool radiatively, and form stars. The interplay between gravity, gas dynamics, star formation, and feedback processes shapes the observable galaxy population.

6.1 Hierarchical Structure Formation

The bottom-up assembly of structure in a CDM universe proceeds through a sequence of halo mergers. The extended Press–Schechter formalism (Lacey & Cole 1993) provides the merger rate:

$$\frac{d^2 N_{\text{merger}}}{dM_2\,dt} = \sqrt{\frac{2}{\pi}}\,\frac{M_1}{M_2^2}\left|\frac{d\delta_c}{dt}\right|\frac{d\sigma_2}{dM_2}\frac{1}{(\sigma_2^2 - \sigma_1^2)^{3/2}}\exp\!\left[-\frac{\delta_c^2(\sigma_2^2 - \sigma_1^2)}{2\sigma_1^2\sigma_2^2}\right]$$

where $M_1$ is the descendant halo mass, $M_2 < M_1$ is the progenitor mass, and $\sigma_i = \sigma(M_i)$. Mergers are classified as major (mass ratio$\gtrsim 1:3$) or minor ($\lesssim 1:10$).

6.2 The NFW Halo Profile

N-body simulations show that dark matter halos in the CDM paradigm are well described by the Navarro–Frenk–White (1996, 1997) density profile:

$$\rho_{\text{NFW}}(r) = \frac{\rho_0}{\left(\dfrac{r}{r_s}\right)\left(1 + \dfrac{r}{r_s}\right)^2}$$

The enclosed mass within radius $r$:

$$M_{\text{NFW}}(<r) = 4\pi\,\rho_0\,r_s^3\left[\ln\!\left(1 + \frac{r}{r_s}\right) - \frac{r/r_s}{1 + r/r_s}\right]$$

The profile has two parameters: the scale radius $r_s$ (where the logarithmic slope equals $-2$) and the characteristic density $\rho_0$. Equivalently, one specifies the virial mass $M_{200}$ and the concentration$c = R_{200}/r_s$, with typical values $c \sim 3\text{--}15$depending on mass and redshift.

The inner slope $\rho \propto r^{-1}$ and outer slope $\rho \propto r^{-3}$are universal features of CDM halos. The concentration–mass relation$c(M) \propto M^{-0.1}$ reflects the fact that less massive halos collapse earlier when the universe was denser.

6.3 Gas Cooling and Star Formation

Baryons initially follow the dark matter, but once inside a halo, they must cool radiatively to reach the high densities needed for star formation. The cooling time is:

$$t_{\text{cool}} = \frac{3\,n\,k_B\,T}{2\,n_e\,n_H\,\Lambda(T)} \sim \frac{3\,k_B\,T}{2\,n\,\Lambda(T)}$$

where $\Lambda(T)$ is the cooling function (dominated by bremsstrahlung at$T > 10^7\,\text{K}$ and metal-line cooling at $10^4 < T < 10^7\,\text{K}$).

Galaxies form efficiently only in halos where $t_{\text{cool}} < t_{\text{ff}}$(the free-fall time). This condition is satisfied for halo masses$M \sim 10^{10}\text{--}10^{12}\,M_\odot$, naturally explaining why the galaxy formation efficiency peaks at $M_{\text{halo}} \sim 10^{12}\,M_\odot$ (the Milky Way mass). In more massive halos, cooling is inefficient (long cooling times in the hot virialized atmosphere), while in less massive halos, supernova feedback expels gas efficiently.

6.4 The Cosmic Star Formation Rate Density: Madau–Dickinson Plot

The cosmic star formation rate density (SFRD) — the total mass of stars formed per unit comoving volume per unit time — has been measured across most of cosmic history. The Madau & Dickinson (2014) compilation shows a characteristic rise-and-fall pattern:

$$\dot{\rho}_*(z) = 0.015\,\frac{(1+z)^{2.7}}{1 + \left[(1+z)/2.9\right]^{5.6}} \quad M_\odot\,\text{yr}^{-1}\,\text{Mpc}^{-3}$$

The SFRD rises steeply from early times, peaks at $z \approx 1.5\text{--}2$(“cosmic noon”, roughly 10 Gyr ago), and declines by a factor of $\sim 10$to the present day. The integrated stellar mass density is$\rho_* = \int_0^\infty \dot{\rho}_*(z)\,(1-R)\,|dt/dz|\,dz \approx 5 \times 10^8\,M_\odot\,\text{Mpc}^{-3}$, where $R \approx 0.4$ is the return fraction (mass recycled by stellar winds and supernovae).

6.5 Downsizing

A key observational puzzle is downsizing (Cowie et al. 1996): the most massive galaxies formed the bulk of their stars earliest and over the shortest timescales, despite residing in halos that assembled last in the hierarchical CDM framework. Specifically:

• Massive ellipticals ($M_* \gtrsim 10^{11}\,M_\odot$) have old stellar populations (ages $\gtrsim 10\,\text{Gyr}$) with enhanced $[\alpha/\text{Fe}]$ ratios, indicating short star-formation timescales ($\lesssim 1\,\text{Gyr}$).
• The characteristic mass at which galaxies transition from blue to red (the quenching mass) decreases with cosmic time.
• At $z \sim 2\text{--}3$, the most actively star-forming galaxies are the most massive.

Downsizing is naturally explained by mass-dependent quenching: in the most massive halos, AGN feedback efficiently suppresses star formation at early times, while lower-mass galaxies continue forming stars until their gas is consumed or removed by later processes.

6.6 Quenching Mechanisms

Several physical processes can suppress (quench) star formation in galaxies. These broadly divide into internal (mass-dependent) and external (environment-dependent) mechanisms:

Internal (Mass Quenching):

• AGN feedback: Energy and momentum from the central SMBH heats or expels gas (see Section 3.5). The dominant mechanism for massive galaxies ($M_* \gtrsim 10^{10.5}\,M_\odot$).
• Virial shock heating: In halos above a critical mass ($M_{\text{halo}} \gtrsim 10^{12}\,M_\odot$), infalling gas is shock-heated to the virial temperature, forming a stable hot atmosphere that cannot cool efficiently (Dekel & Birnboim 2006).
• Morphological quenching: The growth of a central stellar bulge (via mergers) stabilizes the gas disk against fragmentation, suppressing star formation even in the presence of gas (Martig et al. 2009).

External (Environmental Quenching):

• Ram pressure stripping: As a galaxy moves through the hot ICM of a cluster, the ram pressure $P_{\text{ram}} = \rho_{\text{ICM}}\,v^2$ strips cold gas from the disk. Effective in cluster cores where$P_{\text{ram}} > 2\pi\,G\,\Sigma_*\,\Sigma_{\text{gas}}$ (Gunn & Gott 1972).
• Strangulation (starvation): The hot gas halo of a satellite galaxy is stripped upon infall, cutting off the supply of fresh gas for star formation. The galaxy then slowly exhausts its remaining cold gas over 2–5 Gyr.
• Tidal stripping: Gravitational tides from the host halo strip mass from the outer regions of satellite galaxies.
• Harassment: Repeated high-speed tidal encounters with other cluster members heat the stellar disk and strip gas.

6.7 The Stellar-to-Halo Mass Relation

A key outcome of galaxy formation theory is the relationship between stellar mass and halo mass. Abundance matching — equating the cumulative galaxy stellar mass function with the halo mass function — reveals that galaxy formation efficiency $\epsilon = M_*/(f_b\,M_{\text{halo}})$ is a strongly peaked function of halo mass:

$$\frac{M_*}{M_{\text{halo}}} = 2\,\epsilon_0\left[\left(\frac{M_{\text{halo}}}{M_1}\right)^{-\beta} + \left(\frac{M_{\text{halo}}}{M_1}\right)^{\gamma}\right]^{-1}$$

The efficiency peaks at $M_{\text{halo}} \sim 10^{12}\,M_\odot$ with$M_*/M_{\text{halo}} \approx 0.03$ (about 20% of the available baryons converted to stars). Below this mass, supernova feedback suppresses star formation; above it, AGN feedback and hot-halo cooling suppression dominate.

Simulation: NFW Profile and Cluster Mass

The following Python code computes and plots the NFW density profile, enclosed mass, circular velocity, and hydrostatic mass bias for a massive galaxy cluster halo with$M_{200} = 10^{15}\,M_\odot$ and concentration $c = 5$.

Python

nfw_cluster_profile.py197 lines

#!/usr/bin/env python3
"""
nfw_cluster_profile.py
Compute and plot NFW density profile, enclosed mass, circular velocity,
and hydrostatic mass estimate for a galaxy cluster dark matter halo.
"""
import numpy as np
import matplotlib.pyplot as plt

# ------------------------------------------------------------------
# Physical constants (CGS)
# ------------------------------------------------------------------
G_cgs   = 6.674e-8      # cm^3 g^-1 s^-2
Msun    = 1.989e33       # g
Mpc_cm  = 3.0857e24      # cm
kpc_cm  = 3.0857e21      # cm
kB      = 1.381e-16      # erg K^-1
mp      = 1.673e-24      # g
mu_mol  = 0.6            # mean molecular weight for ICM

# ------------------------------------------------------------------
# NFW profile parameters
# ------------------------------------------------------------------
M200     = 1.0e15        # M200 in solar masses (rich cluster)
c200     = 5.0           # concentration parameter
z_cl     = 0.1           # cluster redshift

# Critical density at redshift z (flat LCDM, H0=70, Omega_m=0.3)
H0       = 70.0          # km/s/Mpc
Omega_m  = 0.3
Omega_L  = 0.7
H_z      = H0 * np.sqrt(Omega_m * (1 + z_cl)**3 + Omega_L)  # km/s/Mpc
H_z_cgs  = H_z * 1e5 / Mpc_cm  # s^-1
rho_crit = 3.0 * H_z_cgs**2 / (8.0 * np.pi * G_cgs)  # g cm^-3
rho_crit_Msun_kpc3 = rho_crit / Msun * kpc_cm**3  # Msun kpc^-3

# R200: radius enclosing 200 * rho_crit
R200_cm  = (3.0 * M200 * Msun / (4.0 * np.pi * 200.0 * rho_crit))**(1.0/3.0)
R200_kpc = R200_cm / kpc_cm
R200_Mpc = R200_kpc / 1000.0

# Scale radius and characteristic density
r_s_kpc  = R200_kpc / c200
r_s_cm   = r_s_kpc * kpc_cm

rho_0    = 200.0 * rho_crit * c200**3 / (3.0 * (np.log(1 + c200) - c200/(1 + c200)))
rho_0_Msun_kpc3 = rho_0 / Msun * kpc_cm**3

print("="*60)
print("NFW Halo Parameters")
print("="*60)
print(f"  M200          = {M200:.2e} Msun")
print(f"  c200          = {c200}")
print(f"  Redshift      = {z_cl}")
print(f"  H(z)          = {H_z:.1f} km/s/Mpc")
print(f"  rho_crit(z)   = {rho_crit:.3e} g/cm^3")
print(f"  R200          = {R200_kpc:.0f} kpc  ({R200_Mpc:.2f} Mpc)")
print(f"  r_s           = {r_s_kpc:.0f} kpc")
print(f"  rho_0         = {rho_0_Msun_kpc3:.3e} Msun/kpc^3")
print("="*60)

# ------------------------------------------------------------------
# NFW profile functions
# ------------------------------------------------------------------
def nfw_density(r_kpc):
    """NFW density in Msun/kpc^3."""
    x = r_kpc / r_s_kpc
    return rho_0_Msun_kpc3 / (x * (1 + x)**2)

def nfw_mass(r_kpc):
    """Enclosed mass within r in Msun."""
    x = r_kpc / r_s_kpc
    return 4.0 * np.pi * rho_0_Msun_kpc3 * r_s_kpc**3 * (np.log(1 + x) - x / (1 + x))

def nfw_v_circ(r_kpc):
    """Circular velocity in km/s."""
    M_enc = nfw_mass(r_kpc) * Msun  # grams
    r_cm  = r_kpc * kpc_cm
    v = np.sqrt(G_cgs * M_enc / r_cm)  # cm/s
    return v / 1e5  # km/s

# ------------------------------------------------------------------
# ICM temperature profile (isothermal beta-model for illustration)
# ------------------------------------------------------------------
T_icm     = 8.0e7   # 8 keV cluster, in Kelvin
beta_icm  = 0.67
r_c_kpc   = 250.0   # core radius

def icm_density_profile(r_kpc):
    """Beta-model electron density (cm^-3)."""
    n_e0 = 3.0e-3  # central electron density
    return n_e0 * (1.0 + (r_kpc / r_c_kpc)**2)**(-1.5 * beta_icm)

def hydrostatic_mass(r_kpc, dr=1.0):
    """Hydrostatic mass estimate at radius r (Msun)."""
    # M_hse = -kT r / (G mu m_p) * (d ln rho / d ln r + d ln T / d ln r)
    # For isothermal: d ln T / d ln r = 0
    r_cm = r_kpc * kpc_cm
    # Numerical derivative of ln(rho) w.r.t. ln(r)
    rho1 = icm_density_profile(r_kpc - dr/2)
    rho2 = icm_density_profile(r_kpc + dr/2)
    dlnrho_dlnr = (np.log(rho2) - np.log(rho1)) / (np.log(r_kpc + dr/2) - np.log(r_kpc - dr/2))
    M_hse = -kB * T_icm * r_cm / (G_cgs * mu_mol * mp) * dlnrho_dlnr
    return M_hse / Msun  # solar masses

# ------------------------------------------------------------------
# Radial arrays
# ------------------------------------------------------------------
r_arr = np.logspace(np.log10(10), np.log10(3000), 500)  # kpc

# Compute profiles
rho_arr   = nfw_density(r_arr)
M_arr     = nfw_mass(r_arr)
vc_arr    = nfw_v_circ(r_arr)
M_hse_arr = np.array([hydrostatic_mass(r) for r in r_arr])

# ------------------------------------------------------------------
# Four-panel plot
# ------------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(13, 10))
fig.patch.set_facecolor('#0f172a')
for ax in axes.flat:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='white')
    for spine in ax.spines.values():
        spine.set_color('#334155')

# Panel 1: Density profile
ax1 = axes[0, 0]
ax1.loglog(r_arr, rho_arr, color='#818cf8', lw=2.5, label='NFW density')
ax1.axvline(r_s_kpc, color='#fbbf24', ls='--', lw=1.5, alpha=0.7, label=f'$r_s$ = {r_s_kpc:.0f} kpc')
ax1.axvline(R200_kpc, color='#f87171', ls=':', lw=1.5, alpha=0.7, label=f'$R_{{200}}$ = {R200_kpc:.0f} kpc')
ax1.set_xlabel('r [kpc]', fontsize=12, color='white')
ax1.set_ylabel(r'$\rho$ [M$_\odot$ kpc$^{-3}$]', fontsize=12, color='white')
ax1.set_title('NFW Density Profile', fontsize=14, color='white', fontweight='bold')
ax1.legend(fontsize=9, facecolor='#1e293b', edgecolor='#475569', labelcolor='white')
ax1.grid(alpha=0.2)

# Panel 2: Enclosed mass
ax2 = axes[0, 1]
ax2.loglog(r_arr, M_arr, color='#34d399', lw=2.5, label='NFW enclosed mass')
ax2.loglog(r_arr, M_hse_arr, color='#f97316', lw=2, ls='--', label='Hydrostatic mass')
ax2.axhline(M200, color='#f87171', ls=':', lw=1.5, alpha=0.7, label=f'$M_{{200}}$ = {M200:.1e} M$_\\odot$')
ax2.axvline(R200_kpc, color='#f87171', ls=':', lw=1, alpha=0.4)
ax2.set_xlabel('r [kpc]', fontsize=12, color='white')
ax2.set_ylabel(r'$M(<r)$ [M$_\odot$]', fontsize=12, color='white')
ax2.set_title('Enclosed Mass Profile', fontsize=14, color='white', fontweight='bold')
ax2.legend(fontsize=9, facecolor='#1e293b', edgecolor='#475569', labelcolor='white')
ax2.grid(alpha=0.2)

# Panel 3: Circular velocity
ax3 = axes[1, 0]
ax3.semilogx(r_arr, vc_arr, color='#38bdf8', lw=2.5)
i_max = np.argmax(vc_arr)
ax3.plot(r_arr[i_max], vc_arr[i_max], 'o', color='#fbbf24', ms=8, zorder=5)
ax3.annotate(f'$v_{{max}}$ = {vc_arr[i_max]:.0f} km/s\nat r = {r_arr[i_max]:.0f} kpc',
             xy=(r_arr[i_max], vc_arr[i_max]), xytext=(r_arr[i_max]*3, vc_arr[i_max]*0.85),
             color='#fbbf24', fontsize=10,
             arrowprops=dict(arrowstyle='->', color='#fbbf24'))
ax3.set_xlabel('r [kpc]', fontsize=12, color='white')
ax3.set_ylabel('$v_c$ [km/s]', fontsize=12, color='white')
ax3.set_title('Circular Velocity', fontsize=14, color='white', fontweight='bold')
ax3.grid(alpha=0.2)

# Panel 4: Hydrostatic mass bias
ax4 = axes[1, 1]
bias = (M_hse_arr - M_arr) / M_arr
# Only plot where both are positive
mask = (M_hse_arr > 0) & (r_arr > 50)
ax4.semilogx(r_arr[mask], bias[mask], color='#e879f9', lw=2.5)
ax4.axhline(0, color='white', ls=':', alpha=0.3)
ax4.axhline(-0.2, color='#f87171', ls='--', alpha=0.5, label='Typical bias ~ -20%')
ax4.axvline(R200_kpc, color='#f87171', ls=':', lw=1, alpha=0.4)
ax4.set_xlabel('r [kpc]', fontsize=12, color='white')
ax4.set_ylabel(r'$(M_{hse} - M_{NFW})/M_{NFW}$', fontsize=12, color='white')
ax4.set_title('Hydrostatic Mass Bias', fontsize=14, color='white', fontweight='bold')
ax4.set_ylim(-0.8, 0.5)
ax4.legend(fontsize=9, facecolor='#1e293b', edgecolor='#475569', labelcolor='white')
ax4.grid(alpha=0.2)

plt.tight_layout(pad=2.0)
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.show()
print("NFW cluster profile plots saved to output.png")

# ------------------------------------------------------------------
# Print physical quantities
# ------------------------------------------------------------------
print(f"\nPeak circular velocity: {vc_arr[i_max]:.1f} km/s at r = {r_arr[i_max]:.0f} kpc")
print(f"M_NFW(R200) = {nfw_mass(R200_kpc):.3e} Msun (should ~ M200)")
r500_approx = R200_kpc * 0.65
print(f"Approximate R500 ~ {r500_approx:.0f} kpc")
print(f"M_NFW(R500) = {nfw_mass(r500_approx):.3e} Msun")
print(f"M_hse(R500) = {hydrostatic_mass(r500_approx):.3e} Msun")
bias_500 = (hydrostatic_mass(r500_approx) - nfw_mass(r500_approx)) / nfw_mass(r500_approx)
print(f"Hydrostatic mass bias at R500: {bias_500:.2%}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation Notes

The upper-left panel shows the NFW density profile with its characteristic$\rho \propto r^{-1}$ cusp at small radii and $\rho \propto r^{-3}$fall-off at large radii. The scale radius $r_s$ and virial radius $R_{200}$are marked. The upper-right panel compares the true NFW enclosed mass with the hydrostatic mass estimate from the isothermal beta-model ICM profile. The lower-left panel shows the circular velocity curve, peaking at $r \approx 2.2\,r_s$. The lower-right panel quantifies the hydrostatic mass bias $(M_{\text{HSE}} - M_{\text{NFW}})/M_{\text{NFW}}$, illustrating the systematic underestimate typical of X-ray cluster mass measurements.

7. Summary: Key Equations of Extragalactic Astronomy

Sérsic Profile

$$I(r) = I_e \exp\!\left[-b_n\left(\left(\frac{r}{r_e}\right)^{1/n} - 1\right)\right]$$

Schechter Luminosity Function

$$\Phi(L)\,dL = \phi^*\left(\frac{L}{L^*}\right)^{\alpha}e^{-L/L^*}\frac{dL}{L^*}$$

Eddington Luminosity

$$L_{\text{Edd}} = \frac{4\pi\,G\,M_\bullet\,m_p\,c}{\sigma_T} \approx 1.26 \times 10^{38}\left(\frac{M_\bullet}{M_\odot}\right)\,\text{erg\,s}^{-1}$$

M–sigma Relation

$$M_\bullet \propto \sigma^{4\text{--}5}$$

Hydrostatic Mass

$$M_{\text{HSE}}(r) = -\frac{k_B\,T\,r}{G\,\mu\,m_p}\left(\frac{d\ln\rho}{d\ln r} + \frac{d\ln T}{d\ln r}\right)$$

SZ Comptonization Parameter

$$y = \int \frac{k_B\,T_e}{m_e\,c^2}\,n_e\,\sigma_T\,dl$$

NFW Density Profile

$$\rho(r) = \frac{\rho_0}{(r/r_s)(1 + r/r_s)^2}$$

Press–Schechter Mass Function

$$n(M)\,dM = \sqrt{\frac{2}{\pi}}\frac{\bar{\rho}}{M^2}\frac{\delta_c}{\sigma}\left|\frac{d\ln\sigma}{d\ln M}\right|e^{-\delta_c^2/2\sigma^2}\,dM$$

Cosmic Star Formation Rate Density

$$\dot{\rho}_*(z) = 0.015\,\frac{(1+z)^{2.7}}{1 + [(1+z)/2.9]^{5.6}} \quad M_\odot\,\text{yr}^{-1}\,\text{Mpc}^{-3}$$

← Part V: Gravitational Waves Course Home Part VII: Observational →