feat(accrete): parametric giant migration model behind -r (Phase 1)

Implements migrate_giant() per research/modern/11-giant-migration.md, replacing
the crude uniform-random relocation in dist_planetary_masses (accrete.cpp). Behind
allow_planet_migration (-r), default OFF -> the default golden path is byte-
identical (verify.sh --determinism ref=294147B, unchanged from pre-migration;
goldens clean; ctest 23/23).

Model: inward drift a_final = a_form*exp(-dt/tau) to an inner trap (~0.05 AU disk
edge), tau = viscous Type II (a^1.5) above the ~0.5 M_Jup gap mass (Crida 2006) or
Tanaka 2002 Type I (M^-1 a^(25/14)) below, with the empirical rate reduction in
MIGRATION_EFFICIENCY; hot orbits (perihelion < ~0.05 AU) tidally circularized; only
bodies > 0.1 M_Jup migrate (terrestrials untouched -> inner architecture preserved).

Determinism: a FIXED two-draw order from the per-system RandomContext and an
ANALYTIC perihelion clamp (no rejection loops). Verified -r byte-identical across
-T1/-T2/-T4/-T8.

Calibration (MIGRATION_EFFICIENCY=1e-4, 600-seed harness --migrate):
  hot-Jupiter occurrence (a<0.1 AU)   0.5%   (matches transit rate Howard12/Fressin13)
  giant a-bins: <0.1:3 0.1-1:35 1-3:39 3-10:504 >10:325  (cold peak preserved)
  mass p50/p90  0.71 / 343 M_earth   (no inflation)   density-sane 0.946
  architecture (-r): peas 0.28, hill 19, sigma_e 0.052, min period-ratio 1.012 (>1)
The 1e-4 efficiency is a calibrated value (not derived): StarGen's snapshot Type
I/II model + the steep isothermal Tanaka rate over-migrate massively, so a strong
suppression is needed to match the observed ~0.5% HJ rate -- consistent with the
known Type-I over-efficiency problem (Ida & Lin 2008a reduce by >=10x; the snapshot
model needs more).

Also adds --migrate to validate_population.py (passes -r) to measure it.

KNOWN LIMITATIONS (documented follow-ups, NOT migration bugs):
- cold:hot ratio is ~169 (target ~10-20) ONLY because StarGen over-produces cold
  giants (~84% of systems vs observed ~10-20%) -- a giant-FREQUENCY issue distinct
  from migration (giant occurrence / metallicity).
- giant eccentricity dichotomy not achieved (<e>~0 hot AND cold): the post-accretion
  Cresswell-Nelson gas damping circularizes migration's broad e. Reconciling needs
  exempting cold giants from the damping, which changes the default path/goldens --
  deferred to its own PR.
- hot Jupiters are mostly 0.3-0.5 M_Jup (Type I cores); the snapshot model does not
  carry a >1 M_Jup core inward (it is Type II by final mass) -- a model limitation.

Stays opt-in (-r) per the plan's phased rollout; default-ON is a later, separately-
gated decision (needs the cold-frequency + eccentricity follow-ups first).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: sirus20x6

accrete.cpp

@@ -747,6 +747,76 @@ void accrete::coalesce_planetesimals(long double a, long double e, long double m
  * @return planet* 
  */
 
+// --- Giant migration (research/modern/11-giant-migration.md; -r only) ----------
+// Parametric inward disk migration of a giant / giant-core. Deterministic: a FIXED
+// number of draws (two) in a FIXED order from the per-system RandomContext, and an
+// ANALYTIC perihelion clamp (no rejection loops), so serial == parallel and the
+// stream is robust to threshold tweaks. Modifies a and e in place. Runs only under
+// allow_planet_migration (-r), so the default golden path never executes it.
+static void migrate_giant(long double &a, long double &e,
+                          long double total_mass_solar, long double stell_mass_ratio,
+                          long double min_distance, RandomContext *rng) {
+  const long double a_form = a;
+  const long double m_p_earth = total_mass_solar * SUN_MASS_IN_EARTH_MASSES;
+
+  // Fixed draw order (always two): available-time fraction, then ecc quantile.
+  const long double u_t = rng->randDouble(0.0L, 1.0L);
+  const long double u_e = rng->randDouble(0.0L, 1.0L);
+
+  // Inner trap: disk inner edge / magnetospheric cavity, never inside min_distance.
+  long double a_stop = MIGRATION_INNER_EDGE_AU;
+  if (a_stop < min_distance) {
+    a_stop = min_distance;
+  }
+
+  // Migration timescale (yr): viscous Type II above the gap-opening mass, faster
+  // Tanaka Type I below it; efficiency folds in the empirical rate reduction.
+  const long double gap_mass_earth =
+      MIGRATION_GAP_MASS_MJUP * JUPITER_MASS * stell_mass_ratio;
+  long double tau_yr;
+  if (m_p_earth >= gap_mass_earth) {
+    tau_yr = MIGRATION_TYPE2_NORM_YR * std::pow(a_form, 1.5L) / sqrt(stell_mass_ratio);
+  } else {
+    tau_yr = MIGRATION_TYPE1_NORM_YR / m_p_earth * std::pow(a_form / 10.0L, 25.0L / 14.0L);
+  }
+  tau_yr /= MIGRATION_EFFICIENCY;
+  if (tau_yr <= 0.0L) {
+    return;
+  }
+
+  // Available migration time = fraction of disk lifetime (early formers get the
+  // full budget). The spread carves the hot/warm/cold distribution.
+  const long double dt = u_t * DISK_LIFETIME_YEARS;
+  long double a_new = a_form * std::exp(-dt / tau_yr);
+  if (a_new < a_stop) {
+    a_new = a_stop;
+  }
+  if (a_new > a_form) {
+    a_new = a_form;  // inward only
+  }
+
+  // Eccentricity: broad Rayleigh for cold/warm giants; tidal circularization for
+  // hot orbits whose perihelion crosses the circularization locus.
+  long double e_broad = GIANT_ECC_SIGMA * sqrt(-2.0L * std::log(1.0L - u_e));
+  if (e_broad > GIANT_ECC_MAX) {
+    e_broad = GIANT_ECC_MAX;
+  }
+  long double e_new = e_broad;
+  if (a_new * (1.0L - e_broad) < MIGRATION_CIRC_AU) {
+    e_new = 0.0L;  // tides circularize hot Jupiters
+  }
+  // Analytic clamp: keep perihelion >= min_distance without a rejection loop.
+  if (a_new > 0.0L && a_new * (1.0L - e_new) < min_distance) {
+    e_new = 1.0L - (min_distance / a_new);
+    if (e_new < 0.0L) {
+      e_new = 0.0L;
+    }
+  }
+
+  a = a_new;
+  e = e_new;
+}
+
 auto accrete::dist_planetary_masses(sun &the_sun, long double inner_dust,
                            long double outer_dust,
                            long double outer_planet_limit,
@@ -902,45 +972,15 @@ auto accrete::dist_planetary_masses(sun &the_sun, long double inner_dust,
           else {
             min_distance = 0.015;
           }
-          long double new_a = 0;
-          long double new_e = 0;
-          new_a = random_number(min_distance, a);
-          new_e = random_eccentricity(ecc_coef);
-          for (int i = 0;(calcPerihelion(new_a, new_e) < min_distance) && (i < 1000); i++) {
-            new_a = random_number(min_distance, a);
-            new_e = random_eccentricity(ecc_coef);
-          }
-          if (total_mass >= ((0.7 * JUPITER_MASS) / SUN_MASS_IN_EARTH_MASSES)) {
-            if (random_ctx->randInt(14) == 0 && min_distance < 0.2) {
-              new_a = random_number(min_distance, 0.2);
-              new_e = random_eccentricity(ecc_coef);
-              for (int i = 0;(calcPerihelion(new_a, new_e) < min_distance) && (i < 1000); i++) {
-                new_a = random_number(min_distance, 0.2);
-                new_e = random_eccentricity(ecc_coef);
-              }
-            }
-            if (random_ctx->randInt(10) == 0 &&
-                calcPerihelion(new_a, 0.1) > min_distance) {
-              while (new_e < 0.1) {
-                new_e = random_eccentricity(0.25);
-                for (int i = 0;(calcPerihelion(new_a, new_e) < min_distance) && (i < 1000); i++) {
-                  new_a = random_number(min_distance, a);
-                  new_e = random_eccentricity(ecc_coef);
-                }
-              }
-            }
-          } else if (total_mass >= (2.0 / SUN_MASS_IN_EARTH_MASSES) &&
-                     total_mass <= (10.0 / SUN_MASS_IN_EARTH_MASSES) &&
-                     a >= (4.0 * sqrt(stell_luminosity_ratio)) &&
-                     (dust_mass / total_mass) >= 0.75) {
-            if (random_ctx->randInt(4) == 0) {
-              new_a = random_number(min_distance, a / 2.0);
-              new_e = random_eccentricity(ecc_coef);
-              for (int i = 0;calcPerihelion(new_a, new_e) < planet_inner_bound && i < 1000;i++) {
-                new_a = random_number(min_distance, a / 2.0);
-                new_e = random_eccentricity(ecc_coef);
-              }
-            }
+          long double new_a = a;
+          long double new_e = e;
+          // Only giants / giant-cores migrate; terrestrials stay at their
+          // formation orbit (protects the inner peas-in-a-pod architecture).
+          // See research/modern/11-giant-migration.md and migrate_giant() above.
+          if ((total_mass * SUN_MASS_IN_EARTH_MASSES) >=
+              (MIGRATION_MIN_MASS_MJUP * JUPITER_MASS)) {
+            migrate_giant(new_a, new_e, total_mass, stell_mass_ratio, min_distance,
+                          random_ctx);
           }
           a = new_a;
           e = new_e;

const.h

@@ -186,6 +186,27 @@ constexpr double DISK_LIFETIME_YEARS       = 3.0E6;
 constexpr double DISK_SIGMA0_SOLAR_PER_AU2 = 1.0E-6;
 constexpr double DISK_DAMP_BRACKET_FLOOR   = 0.10;
 
+/* Giant migration (research/modern/11-giant-migration.md; behind -r, default off).
+ * A parametric inward drift a_final = a_form*exp(-dt/tau) to an inner trap, with
+ * tidal circularization of hot orbits. tau follows Tanaka 2002 Type I
+ * (tau ~ M_p^-1 a^(25/14)) below the gap-opening mass and a viscous Type II
+ * (tau ~ a^1.5) above it (Crida 2006 q_c~5e-4 -> ~0.5 M_Jup); the empirical
+ * efficiency reduction (Ida & Lin 2008a) is folded into MIGRATION_EFFICIENCY. The
+ * inner trap is the disk inner edge / magnetospheric cavity (~0.05 AU; Lin,
+ * Bodenheimer & Richardson 1996) and the tidal-circularization locus is ~0.04-0.05
+ * AU (Wang 2023). Cold/warm giants keep a broad Rayleigh eccentricity
+ * (<e>~0.27; Kipping 2013). These are CALIBRATION knobs tuned against the
+ * Phase-0 population baseline (scripts/validate_population.py --migrate). */
+constexpr double MIGRATION_MIN_MASS_MJUP = 0.1;     // only bodies > 0.1 M_Jup migrate
+constexpr double MIGRATION_GAP_MASS_MJUP = 0.5;     // Type II gap-opening threshold (per M_star)
+constexpr double MIGRATION_TYPE1_NORM_YR = 1.1E5;   // Tanaka 2002 Type I normalization (yr)
+constexpr double MIGRATION_TYPE2_NORM_YR = 7.0E5;   // Type II viscous normalization at 1 AU (yr)
+constexpr double MIGRATION_EFFICIENCY    = 0.0001;     // overall rate multiplier (primary calibration knob)
+constexpr double MIGRATION_INNER_EDGE_AU = 0.05;    // inner trap a_stop (disk edge / cavity)
+constexpr double MIGRATION_CIRC_AU       = 0.05;    // tidal circularization locus (hot -> e=0)
+constexpr double GIANT_ECC_SIGMA         = 0.215;   // Rayleigh sigma -> <e>~0.27 for cold giants
+constexpr double GIANT_ECC_MAX           = 0.80;    // eccentricity cap
+
 // SI-unit physical constants used by the enviro.cpp acceleration helpers. These
 // are the EXACT values those functions previously redefined locally, kept
 // verbatim so generated output stays byte-identical. They intentionally differ

scripts/validate_population.py

@@ -43,7 +43,7 @@ def repo_root() -> str:
     return os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
 
 
-def generate(bin_path: str, work: str, seed: int, mass: float) -> dict | None:
+def generate(bin_path: str, work: str, seed: int, mass: float, migrate: bool = False) -> dict | None:
     """Run stargen for one seed and return the parsed system, or None on failure."""
     html = os.path.join(work, "html")
     os.makedirs(html, exist_ok=True)
@@ -56,8 +56,11 @@ def generate(bin_path: str, work: str, seed: int, mass: float) -> dict | None:
         os.remove(out)
     except FileNotFoundError:
         pass
+    argv = [bin_path, f"-s{seed}", f"-m{mass}", "-JS", "-o", "sys"]
+    if migrate:
+        argv.append("-r")  # enable giant migration (research/modern/11-giant-migration.md)
     rc = subprocess.run(
-        [bin_path, f"-s{seed}", f"-m{mass}", "-JS", "-o", "sys"],
+        argv,
         cwd=work, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL,
     )
     if rc.returncode != 0 or not os.path.exists(out):
@@ -311,6 +314,7 @@ def main() -> int:
     ap.add_argument("--mass", type=float, default=1.0)
     ap.add_argument("--json", action="store_true", help="emit metrics as JSON")
     ap.add_argument("--check", action="store_true", help="exit nonzero on gross-sanity violation")
+    ap.add_argument("--migrate", action="store_true", help="enable giant migration (passes -r)")
     args = ap.parse_args()
 
     if not os.path.exists(args.bin):
@@ -323,7 +327,7 @@ def main() -> int:
     systems = []
     with tempfile.TemporaryDirectory(dir=os.environ["TMPDIR"]) as work:
         for i in range(args.seeds):
-            s = generate(args.bin, work, args.seed_start + i, args.mass)
+            s = generate(args.bin, work, args.seed_start + i, args.mass, args.migrate)
             if s is not None:
                 systems.append(s)