Project 2: Science Background

COMP 536 | Short Projects

Author

Dr. Anna Rosen

Published

April 22, 2026

Course ethos: Glass-box modeling. Every function is physics-motivated and transparent. This document is your end-to-end reference for sampling positions from a Plummer sphere and assigning virial equilibrium velocities for star-cluster simulations.

Why the Plummer Model?

The Plummer model is a simple, analytic, finite-mass, non-singular density profile that approximates bound, relaxed stellar systems. It gives closed-form cumulative mass and an exact inverse-CDF for Monte Carlo sampling of radii.

Learning goals:

  1. Translate physics \(\to\) probability distributions \(\to\) samplers.
  2. Build a correct, testable sampler for a spherical density profile.
  3. Understand virial equilibrium and how to initialize a stable gravitational system.

1. Plummer Sphere for Spatial Positions

1.1 Density and Cumulative Mass

The Plummer density with total mass \(M\) and scale radius \(a\) is \[ \rho(r) = \frac{3M}{4\pi a^3}\left(1+\frac{r^2}{a^2}\right)^{-5/2}. \] The enclosed mass is \[ M(<r) = M\,\frac{r^3}{\big(r^2+a^2\big)^{3/2}}. \] Define the CDF for radius: \(F(r)\equiv M(<r)/M\in[0,1)\).

1.2 Inverse-Transform Sampling for the Radius

Set \(F(r)=u\) with \(u\sim\mathcal U(0,1)\). Solve for \(r\): \[ r = a\,\Big(u^{-2/3}-1\Big)^{-1/2}. \] > Draw \(u\in[10^{-12},1-10^{-12})\) to avoid infinities.

1.3 Sample Angles Uniformly on a Sphere

Draw \(\phi=2\pi v\) with \(v\sim\mathcal U(0,1)\) and \(\cos\theta=1-2w\) with \(w\sim\mathcal U(0,1)\). Then \[ x=r\sin\theta\cos\phi,\quad y=r\sin\theta\sin\phi,\quad z=r\cos\theta. \]

1.4 Scale by Half-Mass Radius

If you want a specific half-mass radius \(r_{1/2}\): \[ r_{1/2} = \frac{a}{\sqrt{2^{2/3}-1}}\;\;\Rightarrow\;\; a = r_{1/2}\,\sqrt{2^{2/3}-1} \approx 0.766\, r_{1/2}. \]

1.5 Typical Scales for Simulations

For computational simulations with \(N \leq 200\) particles, use much smaller scales to represent a cluster core:

Setting Scale Radius (\(a\)) Recommended Use
Testing (\(N \leq 50\)) 100 AU \(\approx\) 0.0005 pc Quick validation runs
Production (\(N = 100\)–\(200\)) 1000 AU \(\approx\) 0.005 pc Full analysis
Large \(N\) (\(N > 1000\)) 0.1–1 pc Requires vectorized code

1.6 Center-of-Mass Recentering

After sampling \(N\) stars with positions \(\{\boldsymbol{r}_i\}\) and masses \(\{m_i\}\), subtract the mass-weighted center of mass so that \[\sum_i m_i\boldsymbol{r}_i=\boldsymbol{0}.\]

2. Virial Equilibrium Velocities

For a Plummer sphere in virial equilibrium, assign velocities by drawing each component independently from a Gaussian distribution with radius-dependent (1D) velocity dispersion:

\[\sigma_{1D}^2(r) = \frac{GM}{6a} \cdot \left(1 + \frac{r^2}{a^2}\right)^{-1/2}\]

where \(M\) is the total mass and \(a\) is the Plummer scale radius.

Implementation:

For each particle at radius \(r\):

  1. Calculate \(\sigma_{1D}(r)\) using the formula above
  2. Draw \(v_x \sim \mathcal{N}(0, \sigma_{1D})\)
  3. Draw \(v_y \sim \mathcal{N}(0, \sigma_{1D})\)
  4. Draw \(v_z \sim \mathcal{N}(0, \sigma_{1D})\)

Validation check: The 3D velocity dispersion should satisfy: \[\sigma_{3D} = \sqrt{\langle v^2 \rangle} = \sqrt{3} \cdot \sigma_{1D}\]

After assigning velocities, subtract the mass-weighted center-of-mass velocity to ensure \(\sum_i m_i \vec{v}_i = 0\).

This initialization should produce a system near virial equilibrium:

  • Virial ratio \(Q\) close to 1
  • Virial residual \(\Delta\) below 0.01

3. Putting It All Together — Initial Conditions Sampler

Inputs: \(N\), \(a\) (Plummer scale radius), \(m\) (uniform mass per star).

Outputs: arrays of positions \((x_i, y_i, z_i)\) and initial velocities \((v_{x,i}, v_{y,i}, v_{z,i})\).

Algorithm:

  1. Set masses: \(m_i = m\) for all \(i\) (uniform masses). Total mass \(M = N \cdot m\).
  2. For each star: draw \(u \sim \mathcal{U}(0,1)\), compute \(r_i\) via inverse-CDF.
  3. Draw \(v, w \sim \mathcal{U}(0,1)\) \(\to\) compute \(\theta_i, \phi_i\) \(\to\) \((x_i, y_i, z_i)\).
  4. Subtract mass-weighted center of mass from all positions.
  5. For each star: compute \(\sigma_{1D}(r_i)\), draw \((v_{x,i}, v_{y,i}, v_{z,i}) \sim \mathcal{N}(0, \sigma_{1D})\).
  6. Subtract mass-weighted center-of-mass velocity.

4. Validation & Diagnostics

4.1 Plummer Checks

  • CDF check: Verify that \(u_i = r_i^3 / (r_i^2 + a^2)^{3/2}\) is uniformly distributed. Plot the empirical CDF against the theoretical uniform CDF.
  • Shell counts: Number density in spherical shells: \(n \propto (1 + r^2/a^2)^{-5/2}\).
  • Half-mass radius: Compute \(r_{1/2,\rm sample}\) and compare to target.

4.2 Center-of-Mass Checks

  • \(|r_{\rm cm}| < 0.001 \times a\) (position)
  • \(|v_{\rm cm}| < 0.001 \times \sigma_{3D}\) (velocity)

4.3 Virial Check

  • Initial virial residual \(\Delta < 0.01\) (equivalently \(|Q - 1| < 0.01\))

5. Common Pitfalls (and Fixes)

  1. Angular sampling bias. Fix: Sample \(\cos\theta \sim \mathcal{U}(-1,1)\), \(\phi \sim \mathcal{U}(0, 2\pi)\).
  2. Radius divergence at \(u \to 1\). Fix: Clamp \(u\) away from 0 and 1 using np.clip(u, 1e-12, 1-1e-12).
  3. Forgetting center-of-mass recentering. Fix: Subtract mass-weighted COM from both positions and velocities.
  4. Using wrong \(M\) in velocity dispersion. Fix: \(M\) is the total mass of all particles: \(M = N \cdot m\).
  5. Inconsistent units. Fix: All quantities in AU, \(M_\odot\), yr with \(G = 4\pi^2\).

6. Energy and Virial Diagnostics for N-body Systems

6.1 Energy Components

For an N-body gravitational system, monitor these quantities throughout your simulation:

Kinetic Energy (sum over all particles): \[ E_K = \frac{1}{2} \sum_{i=1}^{N} m_i v_i^2 \] where \(v_i^2 = v_{x,i}^2 + v_{y,i}^2 + v_{z,i}^2\) for star \(i\).

Gravitational Potential Energy (sum over unique pairs): \[ W = -\sum_{i<j} \frac{Gm_im_j}{r_{ij}} \] Critical: Count each pair only once using \(i<j\), not \(i \neq j\) which would double-count.

Total Energy: \[ E_{\rm tot} = E_K + W \]

6.2 The Virial Theorem

For a self-gravitating system in equilibrium, the time-averaged kinetic and potential energies satisfy: \[ \langle 2E_K + W \rangle = 0 \]

Define the Virial ratio: \[ Q = \frac{2E_K}{|W|} \]

Define the Virial residual: \[ \Delta = \frac{|2E_K + W|}{|W|} = |Q - 1| \]

Physical Interpretation:

  • \(Q \approx 1\): System in virial equilibrium (bound and stable)
  • \(Q > 1\): System is super-virial (will expand)
  • \(Q < 1\): System is sub-virial (will contract)
  • \(\Delta < 0.01\): Near-equilibrium initialization quality target

6.3 Using Energy Conservation for Debugging

Monitor all three components (\(E_K\), \(W\), \(E_{\rm tot}\)) separately:

  • If \(E_{\rm tot}\) drifts but \(E_K\) is stable \(\to\) potential energy calculation error
  • If \(E_{\rm tot}\) drifts but \(W\) is stable \(\to\) kinetic energy or velocity update error
  • If both drift proportionally \(\to\) force calculation or integration scheme error
  • Sudden jumps \(\to\) likely self-force bug or sign error

Expected Conservation by Method:

  • Euler: \(E_{\rm tot}\) drifts systematically (energy error grows linearly with time)
  • RK4: \(E_{\rm tot}\) drifts slowly but monotonically (very accurate per step, but error accumulates)
  • Leapfrog: \(E_{\rm tot}\) oscillates with bounded amplitude \(\propto \Delta t^2\) — no secular drift. This is the hallmark of a symplectic integrator.

For a properly initialized cluster, target \(\Delta < 0.01\) (equivalently \(|Q - 1| < 0.01\)).

7. Optional Extension Ideas

These extensions are not required but can deepen your understanding:

  • Kroupa IMF mass sampling: Replace uniform masses with a realistic mass distribution using a piecewise power-law Initial Mass Function. See Kroupa (2001) for the standard two-segment form: \(\alpha_1 = 1.3\) for \(m \in [0.08, 0.5)\,M_\odot\) and \(\alpha_2 = 2.3\) for \(m \in [0.5, 150]\,M_\odot\). Implement using inverse-CDF sampling within each segment.

  • Mass segregation: Sample positions conditionally on mass, e.g., \(r(m) \sim a(m)\,(u^{-2/3}-1)^{-1/2}\) with \(a(m) = a_0\,(m/\langle m\rangle)^{-\beta}\).

  • Adaptive timestepping: Reduce \(\Delta t\) when particles are close, increase when they’re far apart.

8. Quick-Reference Formulas

Plummer radius inverse: \[ r(u)= a\,\big(u^{-2/3}-1\big)^{-1/2},\quad u\sim\mathcal U(0,1). \]

Half-mass scaling: \(a = r_{1/2}\,\sqrt{2^{2/3}-1}\approx 0.766\,r_{1/2}\).

Velocity dispersion: \[\sigma_{1D}^2(r) = \frac{GM}{6a}\left(1 + \frac{r^2}{a^2}\right)^{-1/2}\]

Virial ratio: \(Q = 2E_K / |W|\), target \(Q \approx 1\).

Code Verification Checklist

Before you perform production runs.

End of guide.