Dark Matter Computational Programs

Interactive Python and Fortran programs exploring dark matter physics — rotation curves, halo profiles, thermal freeze-out, direct detection, gravitational lensing, N-body simulations, and the halo mass function. All programs run server-side with adjustable parameters.

Galaxy Rotation Curve Decomposition

Python

Decompose a galaxy rotation curve into disk, bulge, gas, and NFW dark matter halo contributions. Uses Freeman (1970) disk model and Hernquist bulge.

Parameters

Disk Mass (×10¹⁰ M☉)

Stellar disk mass

Disk Scale Length (kpc)

Exponential disk scale

NFW ρ₀ (M☉/pc³)

NFW central density

NFW Scale Radius (kpc)

NFW scale radius

rotation_curve.py118 lines

import numpy as np
import matplotlib.pyplot as plt
def _iv_kv_bessel(v, x):
    """Modified Bessel functions I_v, K_v for v=0,1 via series/integral."""
    x = np.asarray(x, dtype=float)
    if v == 0:
        I = np.ones_like(x); term = np.ones_like(x)
        for k in range(1, 40):
            term *= (x/2)**2 / k**2; I += term
        K = np.where(x > 0, -np.log(np.clip(x/2, 1e-30, None)) * I, 0.0)
        term2 = np.ones_like(x); hk = 0.0
        for k in range(1, 40):
            hk += 1.0/k; term2 *= (x/2)**2 / k**2; K += term2 * hk
        return I, K
    else:
        I = x/2 * np.ones_like(x); term = x/2
        for k in range(1, 40):
            term = term * (x/2)**2 / (k*(k+1)); I += term
        K = np.where(x > 0, 1.0/x + x/2 * np.log(np.clip(x/2, 1e-30, None)), 1e30)
        return I, K
def i0(x): return _iv_kv_bessel(0, x)[0]
def i1(x): return _iv_kv_bessel(1, x)[0]
def k0(x): return _iv_kv_bessel(0, x)[1]
def k1(x): return _iv_kv_bessel(1, x)[1]

# === Parameters (adjustable) ===
M_disk = 5.0e10    # Disk mass [solar masses]
R_d = 3.5          # Disk scale length [kpc]
M_bulge = 1.0e10   # Bulge mass [solar masses]
R_b = 0.5          # Bulge effective radius [kpc]
rho_0 = 0.01       # NFW central density [Msun/pc^3]
r_s = 20.0         # NFW scale radius [kpc]
M_gas = 1.0e10     # Gas mass [solar masses]
R_gas = 7.0        # Gas scale length [kpc]

# Constants
G = 4.302e-3  # pc (km/s)^2 / Msun
kpc = 1000.0  # pc

r = np.linspace(0.1, 50, 500)  # radius in kpc
r_pc = r * kpc

# --- Exponential disk (Freeman 1970) ---
y = r / (2 * R_d)
v2_disk = 4 * np.pi * G * (M_disk / (2 * np.pi * R_d**2 * kpc**2)) * R_d * kpc * y**2 * (
    i0(y) * k0(y) - i1(y) * k1(y)
)
v2_disk = np.abs(v2_disk) * kpc

# --- Hernquist bulge ---
M_enc_bulge = M_bulge * (r / (r + R_b))**2
v2_bulge = G * M_enc_bulge / r_pc

# --- Gas disk (exponential) ---
y_g = r / (2 * R_gas)
v2_gas = 4 * np.pi * G * (M_gas / (2 * np.pi * R_gas**2 * kpc**2)) * R_gas * kpc * y_g**2 * (
    i0(y_g) * k0(y_g) - i1(y_g) * k1(y_g)
)
v2_gas = np.abs(v2_gas) * kpc

# --- NFW dark matter halo ---
rho_0_pc = rho_0  # Msun/pc^3
r_s_pc = r_s * kpc
x = r_pc / r_s_pc
M_enc_nfw = 4 * np.pi * rho_0_pc * r_s_pc**3 * (np.log(1 + x) - x / (1 + x))
v2_halo = G * M_enc_nfw / r_pc

# --- Total ---
v2_total = v2_disk + v2_bulge + v2_gas + v2_halo

# --- Plot ---
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0f172a')
for ax in [ax1, ax2]:
    ax.set_facecolor('#1e293b')
    ax.tick_params(colors='#cbd5e1')
    for spine in ax.spines.values():
        spine.set_color('#334155')

# Left: rotation curve decomposition
ax1.plot(r, np.sqrt(np.maximum(v2_disk, 0)), '--', color='#60a5fa', label='Disk', linewidth=1.5)
ax1.plot(r, np.sqrt(np.maximum(v2_bulge, 0)), '--', color='#f59e0b', label='Bulge', linewidth=1.5)
ax1.plot(r, np.sqrt(np.maximum(v2_gas, 0)), '--', color='#34d399', label='Gas', linewidth=1.5)
ax1.plot(r, np.sqrt(np.maximum(v2_halo, 0)), '--', color='#c084fc', label='DM Halo (NFW)', linewidth=1.5)
ax1.plot(r, np.sqrt(np.maximum(v2_total, 0)), '-', color='#f87171', label='Total', linewidth=2.5)
ax1.axhline(y=np.sqrt(v2_total[400]), color='#f8717180', linestyle=':', linewidth=1)
ax1.set_xlabel('Radius [kpc]', color='#e2e8f0', fontsize=12)
ax1.set_ylabel('Circular Velocity [km/s]', color='#e2e8f0', fontsize=12)
ax1.set_title('Galaxy Rotation Curve Decomposition', color='#e2e8f0', fontsize=14, fontweight='bold')
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=10)
ax1.set_xlim(0, 50)
ax1.set_ylim(0, 300)
ax1.grid(True, alpha=0.15, color='#64748b')

# Right: enclosed mass
ax2.semilogy(r, M_disk * (1 - np.exp(-r/R_d) * (1 + r/R_d)), '--', color='#60a5fa', label='Disk', linewidth=1.5)
ax2.semilogy(r, M_enc_bulge, '--', color='#f59e0b', label='Bulge', linewidth=1.5)
ax2.semilogy(r, M_enc_nfw, '--', color='#c084fc', label='DM Halo', linewidth=1.5)
M_total_enc = M_disk * (1 - np.exp(-r/R_d)*(1+r/R_d)) + M_enc_bulge + M_enc_nfw
ax2.semilogy(r, M_total_enc, '-', color='#f87171', label='Total', linewidth=2.5)
ax2.set_xlabel('Radius [kpc]', color='#e2e8f0', fontsize=12)
ax2.set_ylabel('Enclosed Mass [M☉]', color='#e2e8f0', fontsize=12)
ax2.set_title('Enclosed Mass Profile', color='#e2e8f0', fontsize=14, fontweight='bold')
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#cbd5e1', fontsize=10)
ax2.set_xlim(0, 50)
ax2.grid(True, alpha=0.15, color='#64748b')

plt.tight_layout()
plt.savefig('rotation_curve.png', dpi=150, bbox_inches='tight', facecolor='#0f172a')
plt.close()

# Print results
v_flat = np.sqrt(v2_total[-1])
print(f"Asymptotic rotation velocity: {v_flat:.1f} km/s")
print(f"DM fraction at 50 kpc: {v2_halo[-1]/v2_total[-1]*100:.1f}%")
print(f"DM fraction at 8 kpc (solar): {v2_halo[80]/v2_total[80]*100:.1f}%")
print(f"Total mass within 50 kpc: {M_total_enc[-1]:.2e} Msun")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Equations

The total circular velocity at radius r:

$$v_{\text{tot}}^2(r) = v_{\text{disk}}^2 + v_{\text{bulge}}^2 + v_{\text{gas}}^2 + v_{\text{halo}}^2$$

NFW enclosed mass:

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