Long-Run Stability Search Workbook (reusable)

A general-purpose analysis workbook for any completed ExperimentRunner trial with stability detection enabled - not specific to one run. To reuse on a different run, only the RUN_DIR cell below needs to change - everything else derives from that trial's own metadata.txt, header, and stability_candidates.csv (this workbook handles both the older 2-column crowding-only CSV format and the newer 4-column crowding+proximity+accel format automatically).

Checks four things, in order of how likely they are to fool you:

1. Candidate count trend - the raw stability-detection signal(s), plus an explicit EMA-coldstart-warmup check (a real artifact found while building this: candidates can spike right at step 0 simply because the EMA variance hasn't accumulated any real history yet, not because anything is genuinely stable).
2. RMS radius - is the swarm actually settling, or just becoming smoothly chaotic while still dispersing? (long_stability_search's first run found candidate count growing to 67-75% of all particles while RMS radius grew 66x - a likely false positive from exactly this confusion.)
3. Charge segregation - are the two charge types spatially separating ("one colour prevailing" in a local neighborhood), tested directly rather than eyeballed from the viewer.
4. Bulk motion coherence - are particles actually moving in a common direction (a real collective-motion signal), or does it just look that way from one viewing angle?

None of this replaces looking at the visualization yourself - it's meant to give the visual impression a quantitative gut-check before trusting it.

import sys
import csv
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

# vvv THE ONLY CELL THAT SHOULD NEED TO CHANGE TO REUSE THIS ON ANOTHER RUN vvv
RUN_DIR = Path("~/Documents/sim/long_stability_search_v2/long_stability_search_v2_n200k").expanduser()
# ^^^ ------------------------------------------------------------------- ^^^

assert RUN_DIR.exists(), f"{RUN_DIR} does not exist"
header = bio.read_header(RUN_DIR / "simulation_data.bin")
N_TOTAL = header.num_posons + header.num_negons
metadata_text = (RUN_DIR / "metadata.txt").read_text()
print(metadata_text)

Scenario Name: long_stability_search_v2_n200k
Initial Type Model: Experiment Trial (long_stability_search_v2)
Aggregated Posons: 100000
Aggregated Negons: 100000
Logical Time Steps: 3800
Capture Interval: 100
Integration Flag: 0
Sub-steps Per Frame: 4
Use Octree: false
Octree Theta: 0.5
Stability Detection Enabled: true
Stability Window: 15
Stability Variance Threshold (crowding): 0.05
Stability Variance Threshold (proximity): 0.05
Stability Variance Threshold (acceleration): 0.05
Elapsed Seconds: 11362.92958425

-- Cluster ID Boundaries --
Cluster 1 (Cluster 1) : IDs 0 to 199999

1. Candidate count trend (and the EMA coldstart-warmup artifact)

Loads whichever CSV schema is present (2-column: step,candidateCount from before proximity/acceleration existed; 4-column: step,crowdingCandidates, proximityCandidates,accelCandidates).

candidates_path = RUN_DIR / "stability_candidates.csv"
with open(candidates_path) as f:
    reader = csv.DictReader(f)
    fieldnames = reader.fieldnames
    rows = list(reader)

signal_columns = [c for c in fieldnames if c != "step"]
steps = np.array([int(r["step"]) for r in rows])
signal_series = {col: np.array([int(r[col]) for r in rows]) for col in signal_columns}

print(f"Signals present: {signal_columns}")
for col in signal_columns:
    counts = signal_series[col]
    print(f"\n--- {col} ---")
    for pct in [0, 1, 5, 25, 50, 75, 100]:
        idx = min(int(len(steps) * pct / 100), len(steps) - 1)
        print(f"  step {steps[idx]:>5} ({pct:>3}% through run): {counts[idx]:>7} ({100*counts[idx]/N_TOTAL:.1f}% of particles)")

fig, ax = plt.subplots(figsize=(10, 4.5))
for col in signal_columns:
    ax.plot(steps, 100 * signal_series[col] / N_TOTAL, label=col, alpha=0.85)
ax.set_xlabel("raw step")
ax.set_ylabel("% of particles flagged")
ax.set_title("Stability candidate signals vs. raw step")
ax.legend()
plt.tight_layout()
plt.show()

Signals present: ['crowdingCandidates', 'proximityCandidates', 'accelCandidates']

--- crowdingCandidates ---
  step     0 (  0% through run):       0 (0.0% of particles)
  step    38 (  1% through run):       0 (0.0% of particles)
  step   190 (  5% through run):     114 (0.1% of particles)
  step   950 ( 25% through run):    9187 (4.6% of particles)
  step  1900 ( 50% through run):   29823 (14.9% of particles)
  step  2850 ( 75% through run):   57718 (28.9% of particles)
  step  3799 (100% through run):   78344 (39.2% of particles)

--- proximityCandidates ---
  step     0 (  0% through run):   98771 (49.4% of particles)
  step    38 (  1% through run):      67 (0.0% of particles)
  step   190 (  5% through run):    8183 (4.1% of particles)
  step   950 ( 25% through run):   93019 (46.5% of particles)
  step  1900 ( 50% through run):  133781 (66.9% of particles)
  step  2850 ( 75% through run):  155649 (77.8% of particles)
  step  3799 (100% through run):  168101 (84.1% of particles)

--- accelCandidates ---
  step     0 (  0% through run):  110748 (55.4% of particles)
  step    38 (  1% through run):    2718 (1.4% of particles)
  step   190 (  5% through run):   22462 (11.2% of particles)
  step   950 ( 25% through run):  124713 (62.4% of particles)
  step  1900 ( 50% through run):  180860 (90.4% of particles)
  step  2850 ( 75% through run):  194879 (97.4% of particles)
  step  3799 (100% through run):  198236 (99.1% of particles)

Coldstart-warmup check: each signal's EMA variance is zero-initialized, so a particle can look spuriously "stable" for the first few steps simply because its EMA hasn't accumulated any real variability yet - not because anything is genuinely stable. Watch for a spike in the first ~window steps (see Stability Window in the metadata above) followed by a drop - that shape is the coldstart artifact, not a finding. Only trends that persist well past the warmup window are worth trusting.

window_guess = 15
for m in metadata_text.splitlines():
    if m.startswith("Stability Window:"):
        window_guess = int(m.split(":")[1].strip())

warmup_end = min(window_guess * 3, len(steps) - 1)
print(f"Stability window = {window_guess}; treating the first ~{warmup_end} steps as coldstart/warmup.")
for col in signal_columns:
    counts = signal_series[col]
    print(f"{col}: step 0 = {counts[0]} ({100*counts[0]/N_TOTAL:.1f}%), "
          f"step {warmup_end} = {counts[warmup_end]} ({100*counts[warmup_end]/N_TOTAL:.1f}%), "
          f"final = {counts[-1]} ({100*counts[-1]/N_TOTAL:.1f}%)")

Stability window = 15; treating the first ~45 steps as coldstart/warmup.
crowdingCandidates: step 0 = 0 (0.0%), step 45 = 0 (0.0%), final = 78344 (39.2%)
proximityCandidates: step 0 = 98771 (49.4%), step 45 = 217 (0.1%), final = 168101 (84.1%)
accelCandidates: step 0 = 110748 (55.4%), step 45 = 3536 (1.8%), final = 198236 (99.1%)

Limitation worth knowing: stability_candidates.csv records per-step counts for each signal, not which particles were flagged - so this workbook cannot test whether requiring multiple signals to agree simultaneously (shown in a synthetic concept check to eliminate false positives - see but_stability_detection_feature in project memory) would actually help on this real run. That would need per-particle candidate flags exported from the app, which isn't done yet.

2. RMS radius: is the swarm actually settling?

The population-level sanity check. A genuinely high candidate count should correspond to a swarm that's stopped expanding - if RMS radius is still growing while candidates climb, that's the same false-positive pattern found in the first long_stability_search run.

bin_path = RUN_DIR / "simulation_data.bin"
frame_indices, all_positions, all_charges = bio.load_all_frames(bin_path, header=header)
rms = np.array([bio.rms_radius(all_positions[i]) for i in range(len(frame_indices))])

for i in range(0, len(frame_indices), max(1, len(frame_indices) // 10)):
    print(f"baked frame {frame_indices[i]:>5}: RMS radius = {rms[i]:.2f}")
print(f"baked frame {frame_indices[-1]:>5}: RMS radius = {rms[-1]:.2f}  (last)")
print(f"\nRMS radius change over the run: {rms[0]:.2f} -> {rms[-1]:.2f} ({rms[-1]/rms[0]:.2f}x)")

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(frame_indices, rms)
ax.set_xlabel("raw step (baked frames)")
ax.set_ylabel("RMS radius")
ax.set_title("Is the swarm settling (flat) or still dispersing (still rising)?")
plt.tight_layout()
plt.show()

baked frame     0: RMS radius = 62.46
baked frame     3: RMS radius = 101.52
baked frame     6: RMS radius = 160.93
baked frame     9: RMS radius = 238.15
baked frame    12: RMS radius = 330.85
baked frame    15: RMS radius = 437.39
baked frame    18: RMS radius = 555.84
baked frame    21: RMS radius = 684.81
baked frame    24: RMS radius = 823.06
baked frame    27: RMS radius = 969.40
baked frame    30: RMS radius = 1122.71
baked frame    33: RMS radius = 1281.89
baked frame    36: RMS radius = 1445.96
baked frame    38: RMS radius = 1557.06  (last)

RMS radius change over the run: 62.46 -> 1557.06 (24.93x)

3. Charge segregation ("one colour prevailing")

Two independent checks: whether the two charge types' centroids are physically separating (large-scale segregation into two clouds), and whether each particle's local neighborhood is becoming more charge-homogeneous (same-charge neighbors becoming more common than opposite-charge ones - a more local, "clumping by type" signal that a distant centroid-separation check alone could miss).

# NOTE: all_charges is a SINGLE [N] array, not indexed per frame - charge
# never changes for a given particle across a run (see load_all_frames'
# docstring). Index it directly, not as all_charges[i].
pos_mask = all_charges > 0
neg_mask = ~pos_mask

def charge_split_centroids(positions):
    pos_centroid = positions[pos_mask].mean(axis=0) if pos_mask.any() else np.zeros(3)
    neg_centroid = positions[~pos_mask].mean(axis=0) if neg_mask.any() else np.zeros(3)
    return pos_centroid, neg_centroid

centroid_separation = []
pos_subgroup_rms = []
neg_subgroup_rms = []
overall_centroid_drift = []  # distance of the swarm's own centroid from the origin
for i in range(len(frame_indices)):
    pos_c, neg_c = charge_split_centroids(all_positions[i])
    centroid_separation.append(np.linalg.norm(pos_c - neg_c))
    pos_subgroup_rms.append(bio.rms_radius(all_positions[i][pos_mask]))
    neg_subgroup_rms.append(bio.rms_radius(all_positions[i][neg_mask]))
    overall_centroid_drift.append(np.linalg.norm(all_positions[i].mean(axis=0)))
centroid_separation = np.array(centroid_separation)
pos_subgroup_rms = np.array(pos_subgroup_rms)
neg_subgroup_rms = np.array(neg_subgroup_rms)
overall_centroid_drift = np.array(overall_centroid_drift)

print("Centroid separation between poson-cloud and negon-cloud (relative to overall RMS radius),")
print("plus each subgroup's own internal spread and how far the WHOLE swarm has drifted from the origin")
print("(a big centroid-separation number is not meaningful on its own if the whole swarm is also just")
print("translating a long way from where it started - it needs to be compared to the subgroups' own spread):")
for i in range(0, len(frame_indices), max(1, len(frame_indices) // 10)):
    print(f"  baked frame {frame_indices[i]:>5}: separation={centroid_separation[i]:.2f}  "
          f"pos_subgroup_rms={pos_subgroup_rms[i]:.2f}  neg_subgroup_rms={neg_subgroup_rms[i]:.2f}  "
          f"whole_swarm_drift_from_origin={overall_centroid_drift[i]:.2f}  "
          f"separation/pos_subgroup_rms={centroid_separation[i]/max(pos_subgroup_rms[i],1e-9):.4f}")

# Local charge homogeneity: sample a subset of particles each frame (full n^2
# neighbor search isn't needed for a diagnostic check) and look at the charge
# of each sampled particle's single nearest neighbor.
rng = np.random.default_rng(0)
SAMPLE_SIZE = min(2000, all_positions[0].shape[0])
same_charge_neighbor_fraction = []
for i in range(len(frame_indices)):
    pos_i = all_positions[i]
    n = pos_i.shape[0]
    sample_idx = rng.choice(n, size=SAMPLE_SIZE, replace=False)
    same_count = 0
    # Nearest neighbor via octree for speed at this n.
    root = bio.build_octree(pos_i, all_charges, max_leaf_size=8)
    for idx in sample_idx:
        leaf_node = bio._find_leaf_node(pos_i[idx], root)
        leaf_indices = leaf_node.leaf_indices
        if len(leaf_indices) <= 1:
            continue
        best_j, best_d = None, np.inf
        for j in leaf_indices:
            if j == idx:
                continue
            d = np.linalg.norm(pos_i[j] - pos_i[idx])
            if d < best_d:
                best_d, best_j = d, j
        if best_j is not None and all_charges[best_j] == all_charges[idx]:
            same_count += 1
    same_charge_neighbor_fraction.append(same_count / SAMPLE_SIZE)
same_charge_neighbor_fraction = np.array(same_charge_neighbor_fraction)

print("\nFraction of sampled particles whose nearest neighbor shares their charge (0.5 = no segregation, charges randomly mixed locally; >0.5 = same-charge clumping):")
for i in range(0, len(frame_indices), max(1, len(frame_indices) // 10)):
    print(f"  baked frame {frame_indices[i]:>5}: {same_charge_neighbor_fraction[i]:.3f}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(frame_indices, centroid_separation / np.maximum(pos_subgroup_rms, 1e-9))
axes[0].axhline(0, color="gray", linestyle=":", alpha=0.5)
axes[0].set_title("Poson/negon centroid separation, relative to each subgroup's own spread")
axes[0].set_xlabel("raw step"); axes[0].set_ylabel("separation / subgroup RMS radius")

axes[1].plot(frame_indices, same_charge_neighbor_fraction)
axes[1].axhline(0.5, color="gray", linestyle=":", alpha=0.5, label="no segregation (0.5)")
axes[1].set_title("Local same-charge nearest-neighbor fraction")
axes[1].set_xlabel("raw step"); axes[1].set_ylabel("fraction")
axes[1].legend()
plt.tight_layout()
plt.show()

Centroid separation between poson-cloud and negon-cloud (relative to overall RMS radius),
plus each subgroup's own internal spread and how far the WHOLE swarm has drifted from the origin
(a big centroid-separation number is not meaningful on its own if the whole swarm is also just
translating a long way from where it started - it needs to be compared to the subgroups' own spread):
  baked frame     0: separation=0.07  pos_subgroup_rms=62.46  neg_subgroup_rms=62.46  whole_swarm_drift_from_origin=0.17  separation/pos_subgroup_rms=0.0012
  baked frame     3: separation=0.00  pos_subgroup_rms=101.52  neg_subgroup_rms=101.52  whole_swarm_drift_from_origin=0.28  separation/pos_subgroup_rms=0.0000
  baked frame     6: separation=0.00  pos_subgroup_rms=160.92  neg_subgroup_rms=160.93  whole_swarm_drift_from_origin=0.46  separation/pos_subgroup_rms=0.0000
  baked frame     9: separation=0.01  pos_subgroup_rms=238.14  neg_subgroup_rms=238.15  whole_swarm_drift_from_origin=0.44  separation/pos_subgroup_rms=0.0000
  baked frame    12: separation=0.01  pos_subgroup_rms=330.84  neg_subgroup_rms=330.86  whole_swarm_drift_from_origin=0.80  separation/pos_subgroup_rms=0.0000
  baked frame    15: separation=0.01  pos_subgroup_rms=437.39  neg_subgroup_rms=437.39  whole_swarm_drift_from_origin=1.20  separation/pos_subgroup_rms=0.0000
  baked frame    18: separation=0.02  pos_subgroup_rms=555.83  neg_subgroup_rms=555.84  whole_swarm_drift_from_origin=1.60  separation/pos_subgroup_rms=0.0000
  baked frame    21: separation=0.01  pos_subgroup_rms=684.81  neg_subgroup_rms=684.80  whole_swarm_drift_from_origin=2.04  separation/pos_subgroup_rms=0.0000
  baked frame    24: separation=0.02  pos_subgroup_rms=823.06  neg_subgroup_rms=823.06  whole_swarm_drift_from_origin=2.38  separation/pos_subgroup_rms=0.0000
  baked frame    27: separation=0.02  pos_subgroup_rms=969.39  neg_subgroup_rms=969.41  whole_swarm_drift_from_origin=2.83  separation/pos_subgroup_rms=0.0000
  baked frame    30: separation=0.02  pos_subgroup_rms=1122.71  neg_subgroup_rms=1122.71  whole_swarm_drift_from_origin=3.08  separation/pos_subgroup_rms=0.0000
  baked frame    33: separation=0.01  pos_subgroup_rms=1281.88  neg_subgroup_rms=1281.89  whole_swarm_drift_from_origin=3.33  separation/pos_subgroup_rms=0.0000
  baked frame    36: separation=0.02  pos_subgroup_rms=1445.96  neg_subgroup_rms=1445.95  whole_swarm_drift_from_origin=3.82  separation/pos_subgroup_rms=0.0000

Fraction of sampled particles whose nearest neighbor shares their charge (0.5 = no segregation, charges randomly mixed locally; >0.5 = same-charge clumping):
  baked frame     0: 0.281
  baked frame     3: 0.201
  baked frame     6: 0.164
  baked frame     9: 0.174
  baked frame    12: 0.134
  baked frame    15: 0.133
  baked frame    18: 0.118
  baked frame    21: 0.112
  baked frame    24: 0.086
  baked frame    27: 0.070
  baked frame    30: 0.068
  baked frame    33: 0.065
  baked frame    36: 0.063

4. Bulk motion coherence ("moving in similar directions")

HistoricParticle (the baked-frame format) only stores position + charge, not velocity - so this estimates each particle's effective velocity from the position difference between consecutive baked frames, divided by the number of raw steps between them (captureInterval). This is coarser than the true per-substep velocity and will miss fine-grained motion, but a real bulk/collective drift - which is what "all moving in similar directions" visually suggests - should still show up clearly at this resolution.

The metric is a standard polarization order parameter from collective- motion literature: |mean(v)| / mean(|v|), ranging from 0 (velocities point in random directions, no net drift) to 1 (every particle moving in exactly the same direction).

capture_interval = 1
for m in metadata_text.splitlines():
    if m.startswith("Capture Interval:"):
        capture_interval = int(m.split(":")[1].strip())

polarization = []
for i in range(1, len(frame_indices)):
    dt_steps = frame_indices[i] - frame_indices[i - 1]
    effective_vel = (all_positions[i] - all_positions[i - 1]) / max(dt_steps, 1)
    speeds = np.linalg.norm(effective_vel, axis=1)
    mean_vel_mag = np.linalg.norm(effective_vel.mean(axis=0))
    mean_speed = speeds.mean()
    polarization.append(mean_vel_mag / max(mean_speed, 1e-12))
polarization = np.array(polarization)

print("Polarization order parameter (0 = random directions, 1 = fully aligned):")
for i in range(0, len(polarization), max(1, len(polarization) // 10)):
    print(f"  baked frame {frame_indices[i+1]:>5}: {polarization[i]:.4f}")

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(frame_indices[1:], polarization)
ax.axhline(0, color="gray", linestyle=":", alpha=0.5)
ax.set_ylim(-0.05, 1.05)
ax.set_xlabel("raw step")
ax.set_ylabel("polarization |<v>| / <|v|>")
ax.set_title("Bulk motion coherence over the run")
plt.tight_layout()
plt.show()

Polarization order parameter (0 = random directions, 1 = fully aligned):
  baked frame     1: 0.0082
  baked frame     4: 0.0032
  baked frame     7: 0.0021
  baked frame    10: 0.0042
  baked frame    13: 0.0022
  baked frame    16: 0.0036
  baked frame    19: 0.0037
  baked frame    22: 0.0033
  baked frame    25: 0.0018
  baked frame    28: 0.0015
  baked frame    31: 0.0026
  baked frame    34: 0.0014
  baked frame    37: 0.0029

4b. Radial coherence - a different kind of "moving in similar directions"

Global polarization near zero rules out "the whole swarm drifting off in one direction" - but a self-similarly expanding swarm (consistent with the RMS radius growth above) would show near-zero polarization too, while still looking, from any single zoomed-in viewing angle, like "everything is moving the same way": in a small local patch of a radially-expanding cloud, nearby particles' outward-pointing velocity vectors really are nearly parallel, even though the population's velocities point in every direction overall when you average over the whole sphere.

This checks that directly: for each particle, the cosine similarity between its (effective) velocity direction and the direction from the swarm's own centroid to that particle. +1 = every particle moving straight outward (radial expansion), -1 = every particle moving straight inward (collapse), 0 = no radial relationship.

radial_coherence = []
for i in range(1, len(frame_indices)):
    dt_steps = frame_indices[i] - frame_indices[i - 1]
    effective_vel = (all_positions[i] - all_positions[i - 1]) / max(dt_steps, 1)
    centroid = all_positions[i].mean(axis=0)
    outward_dir = all_positions[i] - centroid
    outward_norm = np.linalg.norm(outward_dir, axis=1, keepdims=True)
    vel_norm = np.linalg.norm(effective_vel, axis=1, keepdims=True)
    valid = (outward_norm[:, 0] > 1e-9) & (vel_norm[:, 0] > 1e-9)
    cos_sim = np.sum(
        (effective_vel[valid] / vel_norm[valid]) * (outward_dir[valid] / outward_norm[valid]), axis=1
    )
    radial_coherence.append(cos_sim.mean())
radial_coherence = np.array(radial_coherence)

print("Radial coherence (+1 = expanding outward together, -1 = collapsing inward together, 0 = no radial pattern):")
for i in range(0, len(radial_coherence), max(1, len(radial_coherence) // 10)):
    print(f"  baked frame {frame_indices[i+1]:>5}: {radial_coherence[i]:.4f}")

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(frame_indices[1:], radial_coherence)
ax.axhline(0, color="gray", linestyle=":", alpha=0.5)
ax.set_ylim(-1.05, 1.05)
ax.set_xlabel("raw step")
ax.set_ylabel("radial coherence")
ax.set_title("Is the population expanding outward together (even though global polarization is ~0)?")
plt.tight_layout()
plt.show()

Radial coherence (+1 = expanding outward together, -1 = collapsing inward together, 0 = no radial pattern):
  baked frame     1: 0.3761
  baked frame     4: 0.4013
  baked frame     7: 0.4227
  baked frame    10: 0.4517
  baked frame    13: 0.4907
  baked frame    16: 0.5310
  baked frame    19: 0.5773
  baked frame    22: 0.6201
  baked frame    25: 0.6611
  baked frame    28: 0.6969
  baked frame    31: 0.7271
  baked frame    34: 0.7501
  baked frame    37: 0.7703

Findings for long_stability_search_v2_n200k (2026-07-18)

• Candidate signals show a clear coldstart artifact, then real (but likely false-positive) late-run growth: proximity/accel start absurdly high at step 0 (49%/55% of all 200,000 particles) purely from EMA zero- initialization, crash down by step ~45 as the EMA warms up (0.1%/1.8%), then climb for the rest of the run to 84% (proximity) and 99.1% (acceleration) of all particles by the end. Crowding follows the same late-run-growth shape (0% -> 39%) without the coldstart spike (its EMA input isn't zero at step 0 the same way). Proximity and especially acceleration are more prone to this than crowding was in the first run - acceleration flagging 99% of the swarm means it's not discriminating anything by the end.
• RMS radius grew 24.9x (62.5 -> 1557.1) - the swarm is still dispersing throughout, not settling. Combined with the above, this is the same false-positive pattern as long_stability_search v1, now confirmed on a run with all three signals: none of crowding/proximity/acceleration distinguishes "genuinely stable" from "aging into a smooth, coasting, still-expanding regime" - exactly the physical confound flagged after v1 (net force ~1/r weakens naturally as the swarm's own scale grows).
• No large-scale charge segregation: poson-cloud and negon-cloud centroids stay essentially coincident throughout (separation is <0.15% of either subgroup's own spread, every frame checked) - the two charge types are not separating into distinct regions.
• But a real, striking local pairing signal: the fraction of particles whose nearest neighbor shares their charge drops from 28.1% at step 0 to just 6.3% by the end - i.e. by the end of the run, over 93% of nearest- neighbor pairs are opposite-charge. This is physically sensible (same charge repels, opposite attracts) and a genuinely interesting signal - it looks like real dipole-like pairing is happening, distinct from (and arguably more meaningful than) any of the three z-score stability signals. This is plausibly what "one colour prevailing" looked like in the viewer: not regional segregation, but individual same-charge particles being pushed apart while opposite-charge pairs bond tightly, which can read as "clumps" depending on rendering/zoom.
• Global velocity polarization is near zero (~0.002-0.008) throughout - the swarm is NOT drifting off in one shared direction. But radial coherence climbs steadily from 0.38 to 0.77 - particles are increasingly moving radially outward from the swarm's own center together. This reconciles the visual "moving in similar directions" impression: it's a self-similar radial expansion becoming cleaner/more coherent over time (consistent with the RMS radius growth), not bulk translation and not bound/orbital structure.

Overall read: this run looks like a swarm settling into an increasingly clean, self-similar ballistic expansion (radial coherence rising, RMS radius rising smoothly) rather than forming genuine bound composite structures - except for the local charge-pairing signal, which is worth a closer, dedicated look (e.g. actually finding and tracking individual opposite-charge pairs' separation over time, to see whether any of them stay bound rather than also drifting apart with the general expansion).

Summary (reusable guidance)

The point of this workbook is that these signals (candidate trend past warmup, RMS radius trend, charge segregation, motion polarization, radial coherence) together give a much more reliable read than eyeballing the 3D viewer, especially since the viewer only shows one projection/zoom level at a time and can't show trends over the whole run at once. When re-running on a different trial, re-derive a findings summary like the one above from that trial's actual numbers rather than assuming these same conclusions carry over.