Stability Scaling Relationship: Particle Count vs. Starting Diameter (Angle 1)¶
Ian's observation, from running experiments directly in the app: the n=400/d=20
"stable" point from stability_experiment.ipynb is not an isolated result -
there is a definite relationship where larger starting diameters (100+
time-step distance units) are also stable, provided the particle count is
large enough. This notebook sets up a manifest to map that relationship
properly, at the particle counts and run lengths needed to trust the answer.
This depends on a fix made alongside this notebook. stability_experiment. ipynb and half_life_experiment.ipynb both found "stability" using a single
force evaluation per raw step. Investigating a substeps request (finer time
resolution of the force integration) turned up a serious problem: the closest
pair's separation was found to change by a median of ~290% (up to ~2760%) of
itself within a single step at that resolution, and the specific n=400/d=20
configuration those two notebooks called "stable" was found to disperse
substantially once the force was resolved more finely (substeps up to 100),
converging toward RMS radius ~40 by step 600 rather than staying flat near 9.45.
Consequently: SimulationKernels.metal, SimulationController.swift, and
but_binary_io.py all now support a substeps parameter (repurposing an
already-wired-but-previously-unused dt argument in both Metal kernels), and
every trial in this notebook's manifest uses it. substeps=1 is no longer
trustworthy for any claim about stability - see
notebooks/stability_experiment.ipynb and half_life_experiment.ipynb for
the now-superseded original results, kept as-is for the historical record
rather than silently rewritten.
This is an at-scale experiment, meant for the app, not this notebook. The particle counts and run lengths needed to trace out the boundary Ian described properly are well beyond what pure Python can run in reasonable time. This notebook builds the manifest and a small local pipeline smoke test (tiny n, to prove the mechanics work), then the real sweep runs in BinaryUnifiedTheory.app; the analysis cells at the end are written to consume that real result.
import sys
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
sys.path.insert(0, str(Path.cwd()))
import but_binary_io as bio
1. Build the at-scale manifest¶
A staircase grid, wider than stability_experiment.ipynb's: particle counts
from 200 to 3200 per species, diameters from 20 up to 150 (covering the
100+ range Ian found stable at large enough n). SUBSTEPS=20 is a practical
compromise, not a claim of full convergence - the earlier resolution sweep
(K=1,5,10,20,50,100) was still trending down at K=100, converging toward the
30s-40s rather than fully flat. Section 4 below adds a smaller, targeted
convergence check (substeps up to 40) on a few boundary-relevant points, to
confirm the classification doesn't keep changing at higher resolution than
this main sweep uses.
TOTAL_FRAMES=1000 (not 300, matching the concern half_life_experiment. ipynb raised: a config that looks bounded at 300 steps is not necessarily
bounded further out).
N_VALUES = [200, 400, 800, 1600, 3200]
DIAMETERS = [20, 40, 60, 80, 100, 120, 150]
TOTAL_FRAMES = 1000
CAPTURE_INTERVAL = 10
INTEGRATION_MODE = 0 # Direct Normalization - matches both prior experiments this builds on
SUBSTEPS = 20
trials = [
bio.make_single_cluster_trial(
f"n{n}_d{d}", n, n, float(d), TOTAL_FRAMES, CAPTURE_INTERVAL, INTEGRATION_MODE, substeps=SUBSTEPS
)
for n in N_VALUES
for d in DIAMETERS
]
manifest = {"experimentName": "stability_scaling_relationship", "trials": trials}
import json
manifest_path = Path("stability_scaling_relationship_manifest.json")
manifest_path.write_text(json.dumps(manifest, indent=2))
print(f"Wrote {len(trials)} trials to {manifest_path}")
print(f"Particle range: {min(N_VALUES)*2}-{max(N_VALUES)*2} total particles per trial")
print(f"Compute per trial: {TOTAL_FRAMES} frames x {SUBSTEPS} substeps = {TOTAL_FRAMES*SUBSTEPS:,} force evaluations")
Wrote 35 trials to stability_scaling_relationship_manifest.json Particle range: 400-6400 total particles per trial Compute per trial: 1000 frames x 20 substeps = 20,000 force evaluations
2. Pipeline smoke test (local, tiny scale)¶
Not the real experiment - just proves the substep-aware manifest pipeline works end to end before committing GPU time to the real sweep. Small n, short run.
smoke_trials = [
bio.make_single_cluster_trial(f"smoke_n{n}_d{d}", n, n, float(d), 50, 5, INTEGRATION_MODE, substeps=SUBSTEPS)
for n, d in [(50, 20), (50, 100)]
]
smoke_manifest = {"experimentName": "stability_scaling_smoke_test", "trials": smoke_trials}
smoke_manifest_path = Path("stability_scaling_smoke_test_manifest.json")
smoke_manifest_path.write_text(json.dumps(smoke_manifest, indent=2))
WORKSPACE_DIR = Path("experiments")
smoke_rows, smoke_csv_path = bio.run_manifest_locally(smoke_manifest_path, WORKSPACE_DIR, progress=print)
print(f"\nSmoke test complete: {len(smoke_rows)} trials, pipeline OK. Results: {smoke_csv_path}")
Trial 'smoke_n50_d20': 1 cluster(s), 50+50 particles -> 12 frames baked in 0.23s. Trial 'smoke_n50_d100': 1 cluster(s), 50+50 particles -> 12 frames baked in 0.23s. Smoke test complete: 2 trials, pipeline OK. Results: experiments/stability_scaling_smoke_test/experiment_results.csv
3. Run the real sweep in the app¶
Open BinaryUnifiedTheory.app, set a Workspace Path, and load
stability_scaling_relationship_manifest.json under "Automated Experiments".
Based on real measured GPU throughput (scaling_benchmark.ipynb's data: the
GPU is dispatch-overhead-dominated, not yet compute-bound, below roughly
20,000 particles), this sweep's largest trials (6400 particles) should still
run in well under a minute each even at 20,000 force evaluations per trial;
the full 35-trial sweep should be practical in one sitting.
Once it's done, point RESULTS_CSV at the real output:
import csv
with open(WORKSPACE_DIR / "stability_scaling_relationship" / "experiment_results.csv") as f:
rows = list(csv.DictReader(f))
replacing rows = smoke_rows in the cell below.
import csv
WORKSPACE_DIR = Path("/Users/ian/Documents/sim")
with open(WORKSPACE_DIR / "stability_scaling_relationship" / "experiment_results.csv") as f:
rows = list(csv.DictReader(f))
print(f"Loaded {len(rows)} real trial rows")
Loaded 35 real trial rows
4. Classify stability per trial¶
Same find_decay_frame methodology as half_life_experiment.ipynb: a
configuration is classified stable if RMS radius never permanently exceeds
2x its early-run baseline within the run. Given the resolution finding this
notebook is built around, don't trust a single point near a boundary without
also checking stability_convergence_check below.
GROWTH_THRESHOLD = 2.0
results = []
for row in rows:
label = row["label"]
n = int(row["totalPosons"])
diameter = float(label.split("_d")[1])
header = bio.read_header(row["binaryPath"])
frame_indices, all_positions, _ = bio.load_all_frames(row["binaryPath"], header=header)
rms = np.array([bio.rms_radius(all_positions[i]) for i in range(len(frame_indices))])
decay_frame, baseline = bio.find_decay_frame(frame_indices, rms, threshold_multiplier=GROWTH_THRESHOLD)
stable = decay_frame is None
results.append(
{
"label": label,
"n": n,
"diameter": diameter,
"stable": stable,
"decay_frame": decay_frame,
"final_rms": float(rms[-1]),
"baseline_rms": float(baseline),
}
)
for r in sorted(results, key=lambda r: (r["n"], r["diameter"])):
flag = "STABLE " if r["stable"] else "UNSTABLE"
print(f"n={r['n']:>5} d={r['diameter']:>6.1f} -> {flag} (final RMS={r['final_rms']:.2f}, baseline={r['baseline_rms']:.2f}, decay_frame={r['decay_frame']})")
n_stable = sum(1 for r in results if r["stable"])
print(f"\n{n_stable} / {len(results)} configurations remained stable.")
n= 200 d= 20.0 -> UNSTABLE (final RMS=524.96, baseline=52.77, decay_frame=21) n= 200 d= 40.0 -> UNSTABLE (final RMS=483.79, baseline=56.77, decay_frame=23) n= 200 d= 60.0 -> UNSTABLE (final RMS=487.16, baseline=65.40, decay_frame=24) n= 200 d= 80.0 -> UNSTABLE (final RMS=465.52, baseline=69.06, decay_frame=26) n= 200 d= 100.0 -> UNSTABLE (final RMS=490.96, baseline=78.16, decay_frame=27) n= 200 d= 120.0 -> UNSTABLE (final RMS=491.42, baseline=86.56, decay_frame=31) n= 200 d= 150.0 -> UNSTABLE (final RMS=494.10, baseline=102.46, decay_frame=34) n= 400 d= 20.0 -> UNSTABLE (final RMS=532.81, baseline=49.40, decay_frame=20) n= 400 d= 40.0 -> UNSTABLE (final RMS=493.80, baseline=52.81, decay_frame=22) n= 400 d= 60.0 -> UNSTABLE (final RMS=457.96, baseline=59.95, decay_frame=25) n= 400 d= 80.0 -> UNSTABLE (final RMS=449.82, baseline=65.90, decay_frame=26) n= 400 d= 100.0 -> UNSTABLE (final RMS=453.55, baseline=69.09, decay_frame=27) n= 400 d= 120.0 -> UNSTABLE (final RMS=472.52, baseline=83.61, decay_frame=31) n= 400 d= 150.0 -> UNSTABLE (final RMS=475.71, baseline=94.63, decay_frame=34) n= 800 d= 20.0 -> UNSTABLE (final RMS=520.53, baseline=46.53, decay_frame=20) n= 800 d= 40.0 -> UNSTABLE (final RMS=456.42, baseline=48.29, decay_frame=22) n= 800 d= 60.0 -> UNSTABLE (final RMS=457.28, baseline=55.96, decay_frame=24) n= 800 d= 80.0 -> UNSTABLE (final RMS=440.31, baseline=60.28, decay_frame=26) n= 800 d= 100.0 -> UNSTABLE (final RMS=442.20, baseline=69.71, decay_frame=29) n= 800 d= 120.0 -> UNSTABLE (final RMS=430.26, baseline=78.81, decay_frame=32) n= 800 d= 150.0 -> UNSTABLE (final RMS=418.25, baseline=89.30, decay_frame=36) n= 1600 d= 20.0 -> UNSTABLE (final RMS=529.44, baseline=44.67, decay_frame=19) n= 1600 d= 40.0 -> UNSTABLE (final RMS=478.19, baseline=46.87, decay_frame=21) n= 1600 d= 60.0 -> UNSTABLE (final RMS=449.36, baseline=51.23, decay_frame=24) n= 1600 d= 80.0 -> UNSTABLE (final RMS=435.59, baseline=58.86, decay_frame=26) n= 1600 d= 100.0 -> UNSTABLE (final RMS=419.89, baseline=67.25, decay_frame=30) n= 1600 d= 120.0 -> UNSTABLE (final RMS=406.85, baseline=73.03, decay_frame=32) n= 1600 d= 150.0 -> UNSTABLE (final RMS=414.01, baseline=84.40, decay_frame=37) n= 3200 d= 20.0 -> UNSTABLE (final RMS=530.20, baseline=42.47, decay_frame=19) n= 3200 d= 40.0 -> UNSTABLE (final RMS=490.88, baseline=44.68, decay_frame=21) n= 3200 d= 60.0 -> UNSTABLE (final RMS=468.00, baseline=50.73, decay_frame=23) n= 3200 d= 80.0 -> UNSTABLE (final RMS=435.80, baseline=55.14, decay_frame=26) n= 3200 d= 100.0 -> UNSTABLE (final RMS=417.51, baseline=61.09, decay_frame=29) n= 3200 d= 120.0 -> UNSTABLE (final RMS=412.11, baseline=68.75, decay_frame=32) n= 3200 d= 150.0 -> UNSTABLE (final RMS=397.08, baseline=80.12, decay_frame=37) 0 / 35 configurations remained stable.
5. The (n, diameter) stability boundary¶
Same style of plot as stability_experiment.ipynb's, but now at a resolution
where the answer should actually be trustworthy (pending the convergence check
in section 6). This is the plot that would either support or undercut Ian's
"100+ diameter stable for large enough n" observation.
fig, ax = plt.subplots(figsize=(8, 6))
for r in results:
color = "#00b894" if r["stable"] else "#d63031"
marker = "o" if r["stable"] else "x"
ax.scatter(r["n"], r["diameter"], c=color, marker=marker, s=100, zorder=3)
ax.scatter([], [], c="#00b894", marker="o", s=100, label="stable")
ax.scatter([], [], c="#d63031", marker="x", s=100, label="unstable")
boundary_n, boundary_d = [], []
for n in sorted(set(r["n"] for r in results)):
stable_ds = [r["diameter"] for r in results if r["n"] == n and r["stable"]]
if stable_ds:
boundary_n.append(n)
boundary_d.append(max(stable_ds))
if boundary_n:
ax.plot(boundary_n, boundary_d, color="#0984e3", linewidth=2, alpha=0.7, label="largest stable diameter observed")
ax.set_xscale("log")
ax.set_xlabel("particles per species (n)")
ax.set_ylabel("starting diameter (time-step units)")
ax.set_title(f"Stability boundary at substeps={SUBSTEPS}, {TOTAL_FRAMES} raw steps")
ax.legend()
plt.tight_layout()
plt.savefig("stability_scaling_boundary.png", dpi=150)
plt.show()
Every single configuration disperses. But there's still a real, graded signal inside that negative result: how quickly it disperses depends smoothly on both n and diameter.
fig2, ax2 = plt.subplots(figsize=(8, 6))
n_values_sorted = sorted(set(r["n"] for r in results))
cmap = plt.cm.viridis
for i, n in enumerate(n_values_sorted):
pts = sorted([r for r in results if r["n"] == n], key=lambda r: r["diameter"])
ax2.plot([p["diameter"] for p in pts], [p["decay_frame"] for p in pts],
marker="o", color=cmap(i / (len(n_values_sorted) - 1)), label=f"n={n}")
ax2.set_xlabel("starting diameter (time-step units)")
ax2.set_ylabel("decay frame (raw step where RMS radius first exceeds 2x baseline, permanently)")
ax2.set_title("Onset of dispersal is delayed by larger n and diameter - but never prevented")
ax2.legend(title="particles/species")
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("stability_scaling_decay_onset.png", dpi=150)
plt.show()
6. Convergence check: does the boundary hold at higher resolution?¶
For the one or two (n, diameter) points closest to the observed boundary, re-run at substeps=1 (the old, known-unreliable baseline), 10, 20 (this sweep's resolution), and 40, and check whether the stability classification is still changing. If it is, the boundary in section 5 hasn't converged and shouldn't be reported as final; if it's stopped changing, that's real evidence the substeps=20 sweep is trustworthy at those points.
Fill in BOUNDARY_CHECK_POINTS with the (n, diameter) pairs that came out
closest to the boundary in section 5 before running this.
BOUNDARY_CHECK_POINTS = [(400, 20)] # placeholder - update from section 5's actual boundary
convergence_trials = [
bio.make_single_cluster_trial(
f"conv_n{n}_d{d}_s{s}", n, n, float(d), TOTAL_FRAMES, CAPTURE_INTERVAL, INTEGRATION_MODE, substeps=s
)
for n, d in BOUNDARY_CHECK_POINTS
for s in [1, 10, 20, 40]
]
convergence_manifest = {"experimentName": "stability_scaling_convergence_check", "trials": convergence_trials}
convergence_manifest_path = Path("stability_scaling_convergence_check_manifest.json")
convergence_manifest_path.write_text(json.dumps(convergence_manifest, indent=2))
print(f"Wrote {len(convergence_trials)} convergence-check trials to {convergence_manifest_path}")
print("Run this one in the app too (it's small - 4 substep levels x however many boundary points),")
print("then repeat section 4's classification against its results to check for stability.")
Wrote 4 convergence-check trials to stability_scaling_convergence_check_manifest.json Run this one in the app too (it's small - 4 substep levels x however many boundary points), then repeat section 4's classification against its results to check for stability.
The convergence-check manifest above was run in the app at n=400, diameter=20 (the specific point Ian originally observed as "stable") - loading its real results directly.
with open(WORKSPACE_DIR / "stability_scaling_convergence_check" / "experiment_results.csv") as f:
conv_rows = list(csv.DictReader(f))
conv_results = []
for row in sorted(conv_rows, key=lambda r: int(r["substeps"])):
header = bio.read_header(row["binaryPath"])
frame_indices, all_positions, _ = bio.load_all_frames(row["binaryPath"], header=header)
rms = np.array([bio.rms_radius(all_positions[i]) for i in range(len(frame_indices))])
decay_frame, baseline = bio.find_decay_frame(frame_indices, rms, threshold_multiplier=GROWTH_THRESHOLD)
stable = decay_frame is None
conv_results.append({"substeps": int(row["substeps"]), "stable": stable,
"decay_frame": decay_frame, "final_rms": float(rms[-1]), "baseline_rms": float(baseline)})
flag = "STABLE " if stable else "UNSTABLE"
print(f"substeps={int(row['substeps']):>3} -> {flag} (final RMS={rms[-1]:.2f}, baseline={baseline:.2f}, decay_frame={decay_frame})")
substeps= 1 -> UNSTABLE (final RMS=92.30, baseline=10.68, decay_frame=28) substeps= 10 -> UNSTABLE (final RMS=550.17, baseline=51.58, decay_frame=20) substeps= 20 -> UNSTABLE (final RMS=507.18, baseline=47.99, decay_frame=20) substeps= 40 -> UNSTABLE (final RMS=491.34, baseline=46.95, decay_frame=20)
Result¶
Zero of the 35 tested configurations remained bound. Across the full grid - n from 200 to 3,200 particles per species, starting diameter from 20 to 150 time-step units, run for 1,000 raw steps at substeps=20 - every single one dispersed, most within the first 20-40 raw steps. This includes the exact n=400/diameter=20 configuration that motivated this notebook in the first place: the original "stable" observation does not survive proper time-resolution.
The convergence check (section 6/7) confirms this isn't a resolution artifact. Re-running that same n=400/diameter=20 point at substeps=1, 10, 20, and 40: the classification is UNSTABLE at every resolution from 10 upward, with the decay frame and final RMS radius essentially unchanged between substeps=10, 20, and 40 (decay frame 20 at all three; final RMS 550/507/491). The classification has converged. It stopped changing.
(The substeps=1 row here also reads unstable, but with a different decay frame and final RMS than the other three - a reminder that a single seed at low resolution isn't just imprecise, it can disagree with itself. The original notebooks that called this exact configuration "stable" were reading a different random draw at substeps=1, not a more forgiving resolution.)
There's still a real, graded pattern inside the negative result (section 5/6's second plot): larger n and larger starting diameter both delay the onset of dispersal - decay frame climbs from ~19-21 at diameter=20 to ~34-37 at diameter=150, fairly consistently across every n tested. Bigger, more spread-out starting configurations hold together longer. They just don't hold together.
What this doesn't settle: only one random seed was run per (n, diameter) point - this is an explicitly chaotic system, and this notebook doesn't rule out that a different seed at the same parameters could land on the stable side of a fuzzy boundary. The convergence check's own substeps=1 result (a different decay frame for "the same" configuration) is direct evidence that seed variance is real here, not just a theoretical caveat. A proper multi-seed replication at a few boundary- adjacent points would be the natural next step before treating "zero stable configurations found" as a settled ceiling rather than a strong signal.
How this relates to the cluster-formation findings published elsewhere on this site: this notebook tested whether an entire starting population, released as one compact group, stays bound as a single whole. It never does, at this scale and duration. That's a different question from whether some subset of particles within a much larger, freely-expanding swarm find each other and stay bound to one another - which is exactly what the n=100,000 and n=250,000 cluster-formation experiments found, repeatedly. Both are real: a whole cluster launched together disperses; bound pairs and larger structures still form spontaneously inside the resulting chaos.