Numerical Simulation of the BeeTheory Dark Mass Model

A complete, reproducible account of the numerical integration, baryonic velocity decomposition, χ² minimisation, and implementation choices behind the BeeTheory dark mass simulation.

This technical page explains how to reproduce the BeeTheory numerical simulation of the Milky Way’s hidden mass. It describes the observational data, baryonic model, wave-based density equation, numerical integration, fitting method, convergence tests, and reference code.

The goal is simple: start from a visible Milky Way disk, apply the BeeTheory wave-density model, compute the effective hidden mass, and compare the resulting circular velocity curve with Gaia-era observations.

0. What We Are Computing and Why

The Milky Way rotation curve is the circular velocity Vc(R) of stars as a function of their distance R from the Galactic Center. It is measured far more precisely today than the total mass distribution we can directly see.

The deficit between the observed velocity and what visible baryonic matter predicts is the hidden mass problem. Standard models invoke an invisible particle halo, usually cold dark matter. BeeTheory proposes a different interpretation: every visible mass element radiates a wave field that decays exponentially in three dimensions, and the accumulated field behaves like hidden mass.

The simulation does three things:

1. Computes V_c^bar(R) from the visible disk and bulge using Freeman’s analytical disk formula.
2. Numerically integrates the BeeTheory density ρ_dark(r; ℓ, λ) at each radius, then converts it into V_c^DM(R).
3. Minimises χ²(ℓ, λ) against the Gaia-era Milky Way rotation curve to find representative best-fit parameters.

The simulation is designed to be reproducible. It can run in JavaScript or Python without a specialized astrophysics library.

Notation used throughout

• R is the cylindrical galactocentric radius in the disk plane, in kpc.
• r is the spherical galactocentric radius.
• z is height above the disk.
• ℓ is the BeeTheory coherence length, in kpc.
• λ is the dimensionless wave–mass coupling.
• G is used in units of kpc km² s⁻² M⊙⁻¹.

1. Observational Data — Gaia 2024

The dataset is based on the Gaia-era Milky Way rotation curve. It uses radii R, circular velocities V_c, and uncertainties σ. The original technical version used 16 data points covering R = 4–27.3 kpc.

The important observational facts are:

• V_c(R⊙ = 8 kpc) ≈ 230 km/s, the anchor point near the Sun’s orbit.
• V_c is approximately flat from about 5 to 20 kpc.
• V_c declines in the outer measured disk, reaching about 173 ± 17 km/s near 27.3 kpc in the referenced data.

This requires an effective hidden mass distribution that rises strongly at intermediate radii, then weakens at larger radii.

We exclude the innermost Galaxy because the central bar and non-circular motions dominate there. A simplified axisymmetric model is not reliable inside that region.

3. The BeeTheory Dark Density Equation

3.1 Physical Postulate

Every visible mass element dV at position r′ in the galactic disk generates a gravitational wave field that creates an effective mass density at a field point r:

The total dark density at any point (R,z) is the superposition over all disk elements. Since the source is a disk, the volume density term becomes a surface density term:

The full double integral is:

3.3 Analytical Verification of the Asymptotic Limit

For R_d ≪ r ≪ ℓ, the expression simplifies. The disk is much smaller than the radius, and the coherence length is still much larger than the radius.

Because:

This asymptotic density behaves as ρ ∝ r⁻², which gives M(<r) ∝ r and therefore a flat rotation curve.

5. χ² Minimisation and Parameter Fitting

5.1 Objective Function

The fit minimises the reduced chi-squared:

with N = 16 data points and p = 2 free parameters, ℓ and λ. This gives 14 degrees of freedom.

5.2 Two-Pass Grid Search

A grid search is used rather than gradient descent because the landscape has a long, curved degeneracy valley between λ and ℓ.

• Pass 1: coarse grid over ℓ and λ.
• Pass 2: local refinement around the coarse minimum.

The representative best-fit region is:

5.3 The Degeneracy Valley

The χ² landscape is not a symmetric bowl. It forms an elongated valley because the leading density normalization depends strongly on the coupling strength and only weakly on the coherence length inside the flat-rotation regime.

The outer decline of the rotation curve breaks this degeneracy because smaller ℓ values suppress the effective density earlier.

Data beyond 30 kpc, including globular clusters, stellar streams, halo stars, and satellite galaxies, are essential for constraining ℓ more tightly.

7. Full Reference Code

The following JavaScript reference implementation reproduces the main simulation logic. It is intended for technical validation and should be placed in a proper script environment, not directly inside a standard WordPress content block unless custom scripts are allowed.

// Physical constants
const G = 4.302e-3;      // kpc km² s⁻² M☉⁻¹
const Sig0 = 800.0;      // M☉ pc⁻²
const Rd = 2.6;          // kpc
const Mdisk = 3.5e10;    // M☉
const Mbulge = 1.2e10;   // M☉
const CONV_RHO = (1.989e30 * 5.61e26) / (3.086e21)**3;

// Gaia-era rotation data
const OBS_R = [4,5,6,7,8,9,10,11,12,14,16,18,20,22,24,27.3];
const OBS_V = [220,228,232,231,230,229,228,227,226,224,222,219,215,208,200,173];
const OBS_ERR = [10,8,7,7,6,6,6,6,7,7,8,9,10,11,13,17];

// Baryonic velocity placeholder
function vBaryonic(R) {
  // In the full implementation, this uses Freeman's disk formula
  // with modified Bessel functions I0, I1, K0, K1.
  const vb2 = G * Mbulge / Math.max(R, 0.2);
  return Math.sqrt(Math.max(0, vb2));
}

// BeeTheory dark density
function rhoDark(r, ell, lam) {
  const N = 40;
  const dRp = 30.0 / N;
  let sum = 0;

  for (let i = 0; i < N; i++) {
    const Rp = (i + 0.5) * dRp;
    const Sig = Sig0 * Math.exp(-Rp / Rd);
    const K = (2 * Math.PI * ell / r)
      * Math.sinh(Math.min(r / ell, 30))
      * Math.exp(-Math.min((r + Rp) / ell, 30));

    sum += Sig * Rp * K * dRp;
  }

  return (lam / ell) * sum;
}

// Enclosed dark mass
function enclosedDarkMass(R, ell, lam) {
  const Nr = 30;
  const dr = R / Nr;
  let M = 0;

  for (let j = 0; j < Nr; j++) {
    const rj = (j + 0.5) * dr;
    M += 4 * Math.PI * rj * rj * rhoDark(rj, ell, lam) * dr;
  }

  return M;
}

// Dark and total velocity
function vDM(R, ell, lam) {
  return Math.sqrt(Math.max(0, G * enclosedDarkMass(R, ell, lam) / R));
}

function vTotal(R, ell, lam) {
  const vb = vBaryonic(R);
  const vd = vDM(R, ell, lam);
  return Math.sqrt(vb * vb + vd * vd);
}

// Chi-squared
function chiSquared(ell, lam) {
  let s = 0;

  for (let i = 0; i < OBS_R.length; i++) {
    const dv = (vTotal(OBS_R[i], ell, lam) - OBS_V[i]) / OBS_ERR[i];
    s += dv * dv;
  }

  return s / (OBS_R.length - 2);
}

A complete version should include accurate implementations of the modified Bessel functions I₀, I₁, K₀, and K₁ for the Freeman disk velocity.

References

• Abramowitz, M., Stegun, I. A. — Handbook of Mathematical Functions, Dover, 1972.
• Bovy, J., Rix, H.-W. — A Direct Dynamical Measurement of the Milky Way’s Disk Surface Density Profile, ApJ 779, 115, 2013.
• Freeman, K. C. — On the disks of spiral and S0 galaxies, ApJ 160, 811, 1970.
• McMillan, P. J. — The mass distribution and gravitational potential of the Milky Way, MNRAS 465, 76, 2017.
• Ou, X., Eilers, A.-C., Necib, L., Frebel, A. — The dark matter profile of the Milky Way inferred from its circular velocity curve, MNRAS 528, 693, 2024.
• Pato, M., Iocco, F., Bertone, G. — Dynamical constraints on the dark matter distribution in the Milky Way, JCAP 12, 001, 2015.
• Portail, M. et al. — Dynamical modelling of the Galactic bulge and bar, MNRAS 465, 1621, 2017.