ThreeBody_Stability: reproducible circumstellar stability limits

This repository contains a reproducible, Python 3 workflow for querying and rebuilding the circumstellar stability-limit tables from Quarles et al. (2020), The Astronomical Journal, 159, 80.

The core quantity is the critical semimajor-axis ratio

$$ a_c = a_p / a_\mathrm{bin} $$

where $a_p$ is the planet semimajor axis and $a_\mathrm{bin}$ is the binary semimajor axis. The tables included here give conservative stability limits as a function of binary mass ratio $\mu$, binary eccentricity $e_\mathrm{bin}$, and initial planet inclination $i_p$.

The results from all the simulations are available on Zenodo.org (http://doi.org/10.5281/zenodo.3579202) as a series of compressed tar archives.

What changed relative to the original public repo

The original repository exposed a useful lookup script, but the workflow was hard to reproduce end-to-end. This version separates the project into:

• a small installable Python package, threebody_stability;
• versioned input data tables under src/threebody_stability/data/;
• command-line tools under scripts/;
• extraction code for rebuilding a_crit_Incl[...].txt from raw simulation archives;
• plotting, fitting, and Monte Carlo utilities;
• tests and a minimal GitHub Actions workflow.

The lookup no longer uses the removed/deprecated scipy.interpolate.interp2d. It builds a regular grid from the tables and uses deterministic bilinear interpolation.

Install

git clone https://github.com/saturnaxis/ThreeBody_Stability.git
cd ThreeBody_Stability
python -m venv .venv
source .venv/bin/activate  # Windows: .venv\Scripts\activate
python -m pip install -e .[dev]

A Conda environment is also provided:

conda env create -f environment.yml
conda activate threebody-stability
python -m pip install -e .[dev]

Tutorial notebooks

The refreshed repository includes a short notebook series under notebooks/:

• 00_quick_start.ipynb — basic package and command-line lookup.
• 01_get_ac_interpolation.ipynb — detailed explanation of how the modern get_ac interpolation works.
• 02_plot_stability_maps.ipynb — visual checks of the packaged stability maps.
• 03_raw_archives_and_derived_products.ipynb — raw-archive extraction, fitted surfaces, and seeded Monte Carlo products.

Run them after installing the notebook extras:

python -m pip install -e .[dev,notebooks]
jupyter lab notebooks/

Notebook cells avoid printing absolute local paths, so rerun outputs are cleaner on different machines.

The most important notebook for users who want to understand the lookup calculation is notebooks/01_get_ac_interpolation.ipynb. It loads an a_crit_Incl[...] table, builds the interpolation grid, shows the four surrounding grid points, and reproduces the bilinear interpolation by hand before comparing against the package result.

Quick lookup

Use either the backwards-compatible script:

python scripts/get_ac.py 0.2 0.3 20.0 0

or the installed console command:

threebody-stability get-ac --mu 0.2 --ebin 0.3 --abin 20.0 --inclination 0

Both commands report $a_c$ and the corresponding physical limit $a_c a_\mathrm{bin}$ in AU.

Python API

from threebody_stability import critical_ratio, critical_semimajor_axis

ratio = critical_ratio(mu=0.2, ebin=0.3, inclination=0)
ac_au = critical_semimajor_axis(mu=0.2, ebin=0.3, abin=20.0, inclination=0)
print(ratio, ac_au)

Reproduce the packaged lookup tables

The bundled tables are:

• a_crit_Incl[0].txt
• a_crit_Incl[30].txt
• a_crit_Incl[45].txt
• a_crit_Incl[180].txt

To rebuild a table from a raw Zenodo archive that contains the MaxEcc_Incl[...].tar.gz files:

python scripts/extract_critical_limits.py \
  --archive raw/MaxEcc_Incl[0].tar.gz \
  --inclination 0 \
  --output src/threebody_stability/data/a_crit_Incl[0].txt

Then test the package:

pytest

Reproduce figures and derived products

Plot all four stability maps:

python scripts/plot_stability_maps.py --output figures/stability_maps.png

The plotting helpers use the paper-style Figure 7 color ordering: gist_ncar_r with fixed limits from 0 to 0.8 for $a_c/a_\mathrm{bin}$.

Fit a Holman & Wiegert-style polynomial to a table:

python scripts/fit_stability_surface.py --inclination 0

Estimate the binary-population stability fraction with a fixed random seed:

python scripts/monte_carlo_binary_pdf.py --n-samples 100000 --seed 42

Input-domain notes

The interpolation tables cover:

• $0.001 \le \mu \le 0.999$;
• $0 \le e_\mathrm{bin} \le 0.8$;
• $i_p \in {0^\circ, 30^\circ, 45^\circ, 180^\circ}$.

Values outside that domain are rejected by default. Use --clip only when you intentionally want nearest-boundary behavior.

Reproducibility boundaries

This repository makes the published stability limits reproducible from the tabulated data and from the raw simulation-output archives. Re-running the original hundreds of millions of N-body integrations requires the external modified Mercury6 workflow used in the paper and substantial compute time. The Python code here is organized so that once raw MaxEcc_Incl[...].tar.gz archives exist, the table extraction, interpolation, fitting, plotting, and Monte Carlo post-processing are reproducible.

Citation

If you use this code or the included stability limits, please cite:

@article{Quarles2020,
  author = {{Quarles}, B. and {Li}, G. and {Kostov}, V. and {Haghighipour}, N.},
  title = {Orbital Stability of Circumstellar Planets in Binary Systems},
  journal = {The Astronomical Journal},
  year = 2020,
  volume = {159},
  number = {3},
  eid = {80},
  pages = {80},
  doi = {10.3847/1538-3881/ab64fa},
  archivePrefix = {arXiv},
  eprint = {1912.11019},
  primaryClass = {astro-ph.EP}
}