Fit a galaxy rotation curve with curve_fit: the dark matter in NGC 3198 (Python)

The problem

The rotation curve of the spiral galaxy NGC 3198 stays near 150 km/s out to 44 kpc, where its stars and gas alone would give 65 km/s. The gap is the dark matter. The pulls \(v^2/r\) of superposed masses add, so squared velocities add:

\[ v^2(r) = \Upsilon\, v_{\mathrm{disk}}|v_{\mathrm{disk}}| + v_{\mathrm{gas}}|v_{\mathrm{gas}}| + v_{\mathrm{halo}}^2(r), \qquad v_{\mathrm{halo}}^2(r) = 4\pi G \rho_0 r_c^2 \left[1 - \frac{r_c}{r}\arctan\frac{r}{r_c}\right]. \]

The halo is pseudo-isothermal, density \(\rho_0/(1 + r^2/r_c^2)\). \(\Upsilon\) is the stellar mass-to-light ratio at 3.6 µm: \(v_{\mathrm{disk}}\) is tabulated for \(\Upsilon = 1\), and \(v^2\) scales with mass. The data are SPARC (Lelli et al. 2016), a 111 kB archive of 175 galaxies, pinned by its checksum. You want the parameters from curve_fit, with bounds, and the dark matter fraction at the last point with its error.

Swap in another SPARC file, or your own curve with its disk and gas components.

The code

import warnings
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

warnings.filterwarnings("ignore", message="IProgress not found")   # tqdm, imported by pooch, wants ipywidgets
import pooch

# ---- data: replace the file name with your own galaxy (any *_rotmod.dat of SPARC)
pooch.get_logger().setLevel("WARNING")
(path,) = pooch.retrieve(
    "https://astroweb.case.edu/SPARC/Rotmod_LTG.zip",
    known_hash="sha256:0a80cc90714828cc28b7dd57923576714d209f2490328c087c4a4ad607faf588",
    processor=pooch.Unzip(members=["NGC3198_rotmod.dat"]),
)
r, v, ev, vgas, vdisk, vbul = np.loadtxt(path, usecols=range(6), unpack=True)   # kpc, km/s

# ---- model
G = 4.30091e-6                                     # kpc (km/s)^2 / Msun

def v_halo2(r, rc, rho0):                          # rho0 in Msun/pc^3, 1e9 pc^3 per kpc^3
    return 4 * np.pi * G * rho0 * 1e9 * rc**2 * (1 - rc / r * np.arctan(r / rc))

def model(r, ups, rc, rho0):                       # v|v| keeps the sign of each pull
    return ups * vdisk * np.abs(vdisk) + vgas * np.abs(vgas) + v_halo2(r, rc, rho0)

# ---- fit
y, sigma = v**2, 2 * v * ev                        # error of a square, to first order
popt, pcov = curve_fit(model, r, y, p0=[0.5, 2.0, 0.05], sigma=sigma, absolute_sigma=True,
                       bounds=([0, 0.01, 0], np.inf))
perr = np.sqrt(np.diag(pcov))
corr = pcov / np.outer(perr, perr)
chi2_dof = np.sum(((y - model(r, *popt)) / sigma) ** 2) / (r.size - popt.size)

# ---- dark matter fraction inside the last point: the halo's share of v^2
def f(p):
    return v_halo2(r[-1], *p[1:]) / model(r, *p)[-1]

J = np.empty(popt.size)                            # d f / d parameter
for k in range(popt.size):
    dp = np.zeros_like(popt)
    dp[k] = 1e-6 * popt[k]
    J[k] = (f(popt + dp) - f(popt - dp)) / (2 * dp[k])
f_err = np.sqrt(J @ pcov @ J)

# ---- report and plot
for name, val, err, unit in zip(["Upsilon", "r_c", "rho_0"], popt, perr, ["Msun/Lsun", "kpc", "Msun/pc^3"]):
    print(f"{name:7s} = {val:6.3f} ± {err:.3f} {unit}")
print(f"chi2 / dof = {chi2_dof:.2f}, correlations: Upsilon-r_c {corr[0, 1]:.2f}, "
      f"Upsilon-rho_0 {corr[0, 2]:.2f}, r_c-rho_0 {corr[1, 2]:.2f}")
v_bary = np.sqrt(popt[0] * vdisk**2 + vgas**2)
print(f"at {r[-1]:.2f} kpc: fit {np.sqrt(model(r, *popt))[-1]:.0f} km/s, stars and gas {v_bary[-1]:.0f} km/s, "
      f"dark matter fraction {100 * f(popt):.1f} ± {100 * f_err:.1f} %")

fig, ax = plt.subplots(figsize=(7, 4), dpi=110)
ax.errorbar(r, v, yerr=ev, fmt="o", color="#1f2a44", ms=4, capsize=2, lw=1)
ax.plot(r, np.sqrt(model(r, *popt)), color="#c8553d", lw=1.8)
ax.plot(r, np.sqrt(v_halo2(r, *popt[1:])), color="#2a7f9e", lw=1.8)
ax.plot(r, np.sqrt(popt[0]) * vdisk, color="#8a8f98", lw=1.2)
ax.plot(r, vgas, color="#8a8f98", lw=1.2)
for x, yy, name, color in [(45, 150, "fit", "#c8553d"), (45, 135, "halo", "#2a7f9e"),
                           (17, 88, "stars", "#8a8f98"), (24, 34, "gas", "#8a8f98")]:
    ax.text(x, yy, name, color=color, va="center")
ax.text(49, 8, f"dark matter inside {r[-1]:.0f} kpc: {100 * f(popt):.0f} ± {100 * f_err:.0f} %", ha="right")
ax.set(xlabel="r / kpc", ylabel="v / (km/s)", xlim=(0, 50), ylim=(0, 175))
ax.spines[["top", "right"]].set_visible(False)
plt.show()

Upsilon =  0.514 ± 0.044 Msun/Lsun
r_c     =  2.759 ± 0.332 kpc
rho_0   =  0.049 ± 0.012 Msun/pc^3
chi2 / dof = 1.26, correlations: Upsilon-r_c 0.97, Upsilon-rho_0 -0.98, r_c-rho_0 -1.00
at 44.08 kpc: fit 150 km/s, stars and gas 65 km/s, dark matter fraction 81.2 ± 0.8 %

Rotation speed of NGC 3198 in km/s against radius in kpc. The data stay flat near 150 km/s, the fit follows them, the stellar disk peaks near 100 km/s at 7 kpc and falls, the gas stays below 50 km/s, and the dark halo rises to 135 km/s: 81 % of v squared at 44 kpc.

The knobs

Fitting \(v^2\) makes the model a plain sum of components, linear in \(\Upsilon\). Each SPARC component is the velocity its mass would give alone, signed by its pull. The gas of NGC 3198 has a central hole, and a flat ring, unlike a spherical shell, pulls the matter inside it outward, toward its near side; where that pull wins, the gas velocity is negative. A square would turn it inward; \(v|v|\) keeps it. A bulge adds ups_bul * vbul * np.abs(vbul) with its own mass-to-light ratio; NGC 3198 has none.

The answer turns on \(\Upsilon\). A heavier disk leaves room for a smaller, denser halo, hence the printed correlations of 0.97 or more in absolute value. The cell fixes \(\Upsilon\) at 0.3, 0.5, 0.8, and the maximum disk, scaled until disk and gas alone carry the observed speed where the disk peaks. The other knob is the halo. The NFW halo (Navarro, Frenk, and White 1996), the one simulations of cold dark matter produce, has a density rising as \(1/r\) toward the center, a cusp where the pseudo-isothermal halo has a flat core; the cell fits it as v_nfw2.

print(f"gas velocity negative from {r[vgas < 0][0]:.2f} to {r[vgas < 0][-1]:.2f} kpc")
i = np.argmax(vdisk)
ups_max = (y[i] - vgas[i] * abs(vgas[i])) / vdisk[i] ** 2      # disk and gas fill v^2 where the disk peaks
over = r[ups_max * vdisk * np.abs(vdisk) + vgas * np.abs(vgas) > y]
print(f"maximum disk: Upsilon = {ups_max:.2f} at {r[i]:.2f} kpc, stars and gas above the data at {over.size} radii, "
      f"{over[0]:.2f} to {over[-1]:.2f} kpc")
for ups in [0.3, 0.5, 0.8, ups_max]:
    p, _ = curve_fit(lambda r, rc, rho0: model(r, ups, rc, rho0), r, y, p0=[2.0, 0.05],
                     sigma=sigma, absolute_sigma=True, bounds=([0.01, 0], np.inf))
    q = np.r_[ups, p]
    chi2 = np.sum(((y - model(r, *q)) / sigma) ** 2) / (r.size - p.size)
    print(f"Upsilon = {ups:4.2f} fixed: chi2 / dof = {chi2:5.2f}, dark matter fraction = {100 * f(q):.1f} %")

def v_nfw2(r, rs, rhos):                           # NFW halo, rhos in Msun/pc^3, x = r / r_s
    x = r / rs
    return 4 * np.pi * G * rhos * 1e9 * rs**3 * (np.log1p(x) - x / (1 + x)) / r

def model_nfw(r, ups, rs, rhos):
    return ups * vdisk * np.abs(vdisk) + vgas * np.abs(vgas) + v_nfw2(r, rs, rhos)

p, _ = curve_fit(model_nfw, r, y, p0=[0.5, 10, 0.01], sigma=sigma, absolute_sigma=True, bounds=([0, 0.01, 0], np.inf))
chi2 = np.sum(((y - model_nfw(r, *p)) / sigma) ** 2) / (r.size - p.size)
print(f"NFW halo: Upsilon = {p[0]:.2f}, r_s = {p[1]:.1f} kpc, chi2 / dof = {chi2:.2f}")

gas velocity negative from 1.93 to 3.54 kpc
maximum disk: Upsilon = 1.02 at 6.74 kpc, stars and gas above the data at 13 radii, 0.32 to 6.43 kpc
Upsilon = 0.30 fixed: chi2 / dof =  2.50, dark matter fraction = 85.2 %
Upsilon = 0.50 fixed: chi2 / dof =  1.23, dark matter fraction = 81.5 %
Upsilon = 0.80 fixed: chi2 / dof =  2.36, dark matter fraction = 76.7 %
Upsilon = 1.02 fixed: chi2 / dof =  9.34, dark matter fraction = 74.1 %
NFW halo: Upsilon = 0.46, r_s = 14.8 kpc, chi2 / dof = 1.54

The free fit lands at \(\Upsilon = 0.51\), near the 0.5 of stellar-population models at 3.6 µm, and the fixed fits get worse on both sides. The maximum disk, \(\Upsilon = 1.02\), overshoots the data with stars and gas alone at 13 radii inside 6.5 kpc, at a \(\chi^2\) per degree of freedom of 9.34. The NFW halo reaches 1.54 against 1.26 for the cored one: keep the pseudo-isothermal halo for this galaxy. Between \(\Upsilon = 0.3\) and 0.8 the dark matter fraction moves from 85.2 % to 76.7 %, ten times the ± 0.8 % the code reports. That ± is the statistical error of this model, nothing more.

The \(\chi^2\) per degree of freedom of 1.26 is a rough check only: the error bars of a rotation curve include the asymmetry between the two sides of the disk, which is not noise. The fraction is a ratio of \(v^2\), not of masses. For the spherical halo \(v^2 r/G\) is exactly the mass inside \(r\), for the flat disk and gas only close to it, so 81 % is the halo's share of the mass inside 44 kpc to a good approximation. It is not the fraction of the whole galaxy, whose halo extends far past the last point.

Pitfalls

The cell makes each of the three mistakes once, in the order of the paragraphs below it.

def v_linear(r, ups, rc, rho0):                    # pitfall 1: speeds added
    return np.sqrt(ups) * vdisk + vgas + np.sqrt(v_halo2(r, rc, rho0))
p, _ = curve_fit(v_linear, r, v, p0=[0.5, 2.0, 0.05], sigma=ev, absolute_sigma=True, bounds=([0, 0.01, 0], np.inf))
chi2 = np.sum(((v - v_linear(r, *p)) / ev) ** 2) / (r.size - p.size)
share = v_halo2(r[-1], *p[1:]) / v_linear(r, *p)[-1] ** 2
print(f"speeds added: Upsilon = {p[0]:.2f}, chi2 / dof = {chi2:.2f}, halo share of v^2 = {100 * share:.0f} %")
p, _ = curve_fit(model, r, y, p0=[0.1, 20, 0.001], sigma=sigma, absolute_sigma=True, maxfev=20000)   # pitfall 2: no bounds
chi2 = np.sum(((y - model(r, *p)) / sigma) ** 2) / (r.size - p.size)
print(f"no bounds:    r_c = {p[1]:.2f} kpc, chi2 / dof = {chi2:.2f}")
# pitfall 3: only the diagonal of pcov
print(f"error of the fraction: {100 * np.sqrt(np.sum((J * perr) ** 2)):.1f} % from the diagonal, {100 * f_err:.1f} % from pcov")

speeds added: Upsilon = 0.27, chi2 / dof = 0.82, halo share of v^2 = 21 %
no bounds:    r_c = -2.76 kpc, chi2 / dof = 1.26
error of the fraction: 5.1 % from the diagonal, 0.8 % from pcov

Adding the velocities instead of their squares. Write the model as sqrt(ups) * vdisk + vgas + sqrt(v_halo2) and fit \(v\): the fit converges with a \(\chi^2\) per degree of freedom of 0.82, better than the right model. But \(\Upsilon\) drops to 0.27 and the halo's share of \(v^2\) at the last point to 21 %, against 81 % from the code. A fit that looks good with parameters that make no sense is the symptom. The cause is the rule stated above the model: speeds do not add. Use the sum of signed squares of the code.

A negative core radius without bounds. The halo velocity is even in \(r_c\), so the data cannot tell \(r_c\) from \(-r_c\). Without bounds curve_fit uses Levenberg-Marquardt, and from the start [0.1, 20, 0.001] it lands on \(r_c = -2.76\) kpc with the same \(\chi^2\) and no warning. \(\Upsilon\) has no mirror, since the model is linear in it, but an unbounded fit can still find a separate, unphysical minimum below zero: the dwarf galaxy DDO 168 of the same archive lands at \(\Upsilon = -2.41\). Keep the bounds of the code, with the lower bound of \(r_c\) above zero so that \(r/r_c\) never divides by zero.

The fraction's error from the diagonal alone. Propagate only the three standard errors and the fraction comes out at ± 5.1 % instead of ± 0.8 %, more than six times too large. The parameters move together, as the correlations printed by the code say, and along that direction the fraction barely changes; the diagonal does not know that. Put the whole pcov into \(J \Sigma J^\mathsf{T}\), as the code does.