Reverse N-Body Problem and Astrometric Parallax¶
Status: Foundational reference Relates to: [[TTV Sensitivity and N-Body Math]], [[TTV Reverse N-Body Inference]], NEO distance measurement Last derived: 2026-03-20
Part I: The Reverse N-Body Problem in Full¶
This document expands on [[TTV Sensitivity and N-Body Math]] by treating the inverse problem in full mathematical detail, covering the degeneracy structure, MCMC implementation strategy, and Bayesian model order selection.
1. State Space and Parameterization¶
1.1 Orbital Elements¶
A Keplerian orbit is fully specified by six elements. For the TTV inverse problem, we use:
$$\boldsymbol{\theta}i = {P_i, t, e_i \cos\omega_i, e_i \sin\omega_i, i_i, \Omega_i, m_i/M_\star}$$
where: - $P_i$: period - $t_{0,i}$: time of inferior conjunction (transit epoch) - $e_i \cos\omega_i$, $e_i \sin\omega_i$: eccentricity vector components (ensures uniform prior in $e$, avoids singularity at $e=0$) - $i_i$: orbital inclination (note: TTV constraints become degenerate for near-90° inclinations) - $\Omega_i$: longitude of ascending node - $m_i/M_\star$: mass ratio
For co-planar orbits: $i_i = i_b$ for all planets, $\Omega_i - \Omega_b$ is a single free parameter (relative nodal longitude, often set to 0 for simplicity in 2D problems).
For the minimal 2-planet case (transiting planet $b$ + non-transiting perturber $c$), with co-planar orbits and known $M_\star$, the free parameters are:
$$\boldsymbol{\theta} = {P_b, t_{0,b}, e_b\cos\omega_b, e_b\sin\omega_b, P_c, t_{0,c}, e_c\cos\omega_c, e_c\sin\omega_c, \mu_c}$$
where $\mu_c = m_c/M_\star$ is the perturber mass ratio. 9 free parameters.
Note: $m_b/M_\star$ does not strongly affect the observed TTV because the TTV is driven by the perturber's acceleration of planet $b$, which scales as $m_c$, not $m_b$ (to first order). However, $m_b$ affects the perturber's orbit through the back-reaction; this term enters at second order in $m_b/M_\star$ and is typically small.
1.2 Initial Conditions for the Integrator¶
TTVFast and similar integrators work in Cartesian initial conditions $(x, y, z, \dot{x}, \dot{y}, \dot{z})$, converted from orbital elements at a reference epoch $t_\text{ref}$:
$$\mathbf{r}_i = a_i \left[(1 - e_i^2)/(1 + e_i \cos f_i)\right] \hat{\mathbf{r}}(\Omega_i, \omega_i, i_i, f_i)$$
where $f_i$ is the true anomaly at $t_\text{ref}$, computed from the mean anomaly $M_i = 2\pi (t_\text{ref} - t_{0,i})/P_i$ via Kepler's equation: $$M_i = E_i - e_i \sin E_i, \quad \tan(f_i/2) = \sqrt{\frac{1+e_i}{1-e_i}} \tan(E_i/2)$$
This conversion must be done carefully — errors in Kepler's equation solving propagate directly into the transit time predictions.
2. The Forward Model: TTVFast¶
2.1 Algorithm Summary¶
TTVFast (Deck et al. 2014) uses a Wisdom-Holman symplectic integrator with a fixed timestep $\Delta t$. [Source: Deck et al. (2014), ApJ 787, 132 — TTVFast algorithm and validation] The key feature is that transit times are found by a Hermite interpolation between timesteps bracketing the transit (when the $z$-coordinate of planet $b$ changes sign with $\dot{z} < 0$).
The transit time at transit $n$ is found by solving: $$z(t_n) = 0, \quad \dot{z}(t_n) < 0$$
using Newton's method on the interpolated trajectory. The precision of transit times is: $$\sigma_{t_n}^\text{TTVFast} \approx \left(\frac{\Delta t}{2\pi}\right)^4 \cdot P_\text{short}$$
For $\Delta t = P_\text{short}/10$, the integration error is $\sim (0.1/2\pi)^4 \times P = 2.5\times10^{-6} P$. For $P = 5$ d: $\sim 1$ s error. Acceptable.
Recommended TTVFast timestep: $\Delta t = 0.05 \times \min(P_b, P_c)$
2.2 Runtime Scaling¶
For $N_\text{transits}$ transit epochs over a baseline $T_\text{base}$: - Integration time: $\propto T_\text{base} / \Delta t = T_\text{base} \times f_s / P_\text{short}$ - For $T_\text{base} = 2$ yr, $P_\text{short} = 5$ d, $f_s = 20$: $N_\text{steps} = 730 \times 20 = 14600$ timesteps - Each timestep: $O(n^2)$ force evaluations, $n = 3$ bodies: constant factor ~9 - TTVFast on modern hardware: $\sim 0.5$–2 ms per likelihood call
For MCMC with $10^6$ evaluations: $10^6 \times 2\,\text{ms} = 2000\,\text{s} \approx 30\,\text{min}$ per chain. This is tractable on a single CPU core; multiple chains run in parallel.
3. The Posterior Geometry and Its Challenges¶
3.1 Known Multimodalities¶
Period ratio ambiguity: The TTV super-period $P_\text{TTV}$ constrains the combination $|j/P_c - (j-k)/P_b|$ but not which $j:(j-k)$ resonance is responsible. Multiple period ratios can produce the same super-period, leading to isolated posterior modes.
Phase ambiguity: For a sinusoidal TTV signal with amplitude $A$ and phase $\phi$, two solutions related by $\phi \to \pi - \phi$ can give the same amplitude but different perturber positions. If only part of the TTV cycle is observed (baseline $< P_\text{TTV}$), the posterior is bimodal in phase.
The $m_c$-$e_c$ degeneracy: Higher $e_c$ changes the resonant forcing function $f(j,k,e_b,e_c)$. A lower-mass perturber on an eccentric orbit can produce the same TTV as a higher-mass perturber on a circular orbit. This degeneracy is broken when higher-order (short-period) TTV terms are included — "chopping" signals at $P_c$ that scale differently with $e_c$ vs $m_c$.
3.2 Prior Recommendations¶
Use the following priors (all in $\log$ unless noted):
| Parameter | Prior | Rationale |
|---|---|---|
| $P_b$ | Gaussian centered on known period, $\sigma = 3\times$ measured uncertainty | Informative; known from initial fit |
| $t_{0,b}$ | Gaussian centered on known epoch, same approach | Informative |
| $e_b\cos\omega_b$, $e_b\sin\omega_b$ | Uniform in $[-0.3, 0.3]^2$ | Small eccentricity prior; hot Jupiters circularized |
| $\ln P_c$ | Uniform in $[\ln P_b/3, \ln 10 P_b]$ | Wide; no strong period prior |
| $t_{0,c}$ | Uniform over one period $P_c$ | Phase-uniform |
| $e_c\cos\omega_c$, $e_c\sin\omega_c$ | Uniform in $[-0.5, 0.5]^2$ | Moderate eccentricity allowed |
| $\ln \mu_c$ | Uniform in $[\ln 10^{-7}, \ln 10^{-2}]$ | Earth to Jupiter mass |
The log-uniform prior on $\mu_c$ is physically motivated: we have no a priori preference for a specific mass range, and equal probability per decade of mass is the maximum entropy prior on a scale parameter.
3.3 Identifiability Conditions¶
The 9-parameter model is identifiable (in principle) when: 1. At least 10 transit times are observed ($M \geq 10$) 2. The observation baseline exceeds $P_\text{TTV}/2$ 3. The TTV signal is detected at $> 3\sigma$ (otherwise the posterior is essentially the prior) 4. The "chopping" signal ($P_c$-period variations in the timing residuals after subtracting the sinusoidal TTV) is detectable — this is needed to break the $m_c$-$e_c$ degeneracy
Condition 4 requires: - Baseline $\gg P_c$ (many cycles of $P_c$) - Photometric precision good enough that per-transit $\sigma_{t_m} < P_c \times (m_c/M_\star)$ (i.e., the chopping amplitude is detectable)
4. MCMC Implementation¶
4.1 Sampling Algorithm¶
The multimodal posterior with isolated modes makes standard Metropolis-Hastings inefficient. Recommended algorithms:
Parallel tempering (PT-MCMC):
- Run $N_T$ chains at temperatures $T_1 = 1 < T_2 < \ldots < T_{N_T}$
- Chain $k$ explores the "tempered" posterior $p(\theta)^{1/T_k}$, which is broader (higher T = more exploration)
- Periodically swap configurations between adjacent temperature chains
- The $T=1$ chain yields the correct posterior samples
- Implementation: ptemcee (Python) [Source: standard parallel tempering MCMC; see Vousden et al. (2016) for ptemcee implementation]
Nested sampling:
- Does not get stuck in local modes
- Simultaneously samples the posterior and computes the evidence $Z = \int \mathcal{L} \pi \, d\theta$
- Evidence is needed for model comparison (Section 5)
- Implementation: dynesty or MultiNest (via PyMultiNest) [Source: Speagle (2020), MNRAS 493, 3132 — dynesty; Feroz et al. (2019), OJAp 2, 10 — MultiNest]
- Recommended for TTV analysis because of the need for both posterior samples and evidence
Hamiltonian Monte Carlo (HMC): - Very efficient in parameter space if the posterior is smooth and unimodal - Not recommended for multimodal TTV posteriors without significant posterior pre-conditioning
4.2 Proposal Strategy for MCMC¶
For emcee-based sampling (affine-invariant ensemble sampler): - Walkers: $N_w = \max(100, 10 k_\text{params}) = 90$ walkers for 9 parameters - Initialization: Start walkers near the best-fit from a preliminary optimization (Nelder-Mead or scipy.optimize) with small Gaussian scatter - Burn-in: $\sim 500$–1000 steps (monitor autocorrelation time $\tau_\text{ac}$) - Production: $50 \tau_\text{ac}$ steps for well-converged chains - Convergence check: Gelman-Rubin $\hat{R} < 1.01$ for all parameters across 4 independent chains
4.3 Computational Budget Estimate¶
For TTVFast at 2 ms/call, 9 parameters, target $5\times10^4$ independent samples after thinning:
- Auto-correlation time $\tau_\text{ac} \approx 50$–200 steps (typical for multimodal)
- Total steps needed: $50 \tau_\text{ac} \times N_w = 50 \times 100 \times 90 = 450{,}000$
- Total likelihood calls: $450{,}000$
- Time: $450{,}000 \times 2\,\text{ms} = 900\,\text{s} = 15\,\text{min}$ on 1 core
- With 8 chains in parallel (4 temperatures, 2 replicas each): ~2 min
TTV inference is computationally tractable on a standard laptop in $< 30$ minutes for a 2-planet system. BOINC / cloud compute is only needed for 3+ planet systems or very long baselines requiring many timesteps.
[NOVEL] The specific compute budget analysis (2 ms/call × ~450,000 calls ≈ 15 min per 2-planet TTV inference on a single CPU) and the resulting conclusion that BOINC is unnecessary for Stage 1, is original to OpenAstro.
5. Bayesian Model Comparison: Selecting the Number of Planets¶
5.1 The Evidence Integral¶
For model $\mathcal{M}_n$ ($n$ planets total in the system, with 1 transiting and $n-1$ perturbers), the Bayesian evidence is:
$$Z_n = P(\text{data} | \mathcal{M}_n) = \int \mathcal{L}(\boldsymbol{\theta} | \text{data}) \pi(\boldsymbol{\theta} | \mathcal{M}_n) \, d\boldsymbol{\theta}$$
The Bayes factor comparing $n=3$ to $n=2$: $$B_{32} = \frac{Z_3}{Z_2}$$
Jeffreys' scale interpretation: [Source: Jeffreys (1961), Theory of Probability; Kass & Raftery (1995), JASA 90, 773] - $\ln B_{32} < 1$: weak evidence, models indistinguishable - $1 < \ln B_{32} < 3$: moderate evidence for $n=3$ - $3 < \ln B_{32} < 5$: strong evidence for $n=3$ - $\ln B_{32} > 5$: very strong evidence for $n=3$
5.2 Complexity Penalization¶
The evidence naturally penalizes overfitting (Occam's razor). Adding 5 parameters (${P_c', t_{0,c}', e_c'\cos\omega_c', e_c'\sin\omega_c', \mu_c'}$) to go from $n=2$ to $n=3$ costs:
$$\ln \text{Occam factor} \approx -\frac{k_\text{new}}{2} \ln\left(\frac{\pi_0}{\sigma_\text{posterior}}\right) \approx -\frac{5}{2} \ln\left(\frac{\Delta P}{\sigma_P}\right)$$
where $\Delta P$ is the prior range and $\sigma_P$ is the posterior width. For $\Delta P / \sigma_P \sim 100$ (prior 100 times wider than posterior for each new parameter):
$$\ln \text{Occam factor} \approx -\frac{5}{2} \ln 100 = -11.5$$
So the data must improve $\ln \mathcal{L}$ by at least 11.5 (equivalent to $\chi^2$ reduction of 23) to justify the 5 extra parameters. This corresponds to $\sim 4.8\sigma$ improvement in the maximum likelihood fit.
In practice: Fit $n=2$, compute $\chi^2_\nu$ and residuals. If $\chi^2_\nu > 1.5$ and residuals show periodic structure at a period not present in the $n=2$ model, run nested sampling for both $n=2$ and $n=3$, compute Bayes factor. Report the preferred model.
5.3 Dealing with False Positives¶
A TTV signal can be mimicked by: - Stellar variability (star spots modulating the transit depth asymmetrically) - Instrumental systematics (flat-field variations correlated with field rotation) - Contaminating eclipsing binaries (background EB in the photometric aperture produces a periodic eclipse that shifts the measured midtime)
Validation tests: 1. Check that the TTV period does not correlate with stellar rotation period (from autocorrelation of out-of-transit flux) 2. Check that the TTV signal is consistent across different wavelengths (star spot signals are wavelength-dependent; gravitational TTVs are not) 3. Check that the apparent TTV is not reproduced by varying the comparison star ensemble
6. Degeneracy Breakdown: When Can You Distinguish One Perturber from Two?¶
6.1 The Superposition Principle (First Order)¶
At first order in mass ratios, the TTV signal from multiple perturbers superimposes linearly:
$$\delta t_n = \sum_{c} \delta t_n^{(c)}$$
where $\delta t_n^{(c)}$ is the TTV from perturber $c$ alone. Each perturber near a distinct MMR produces a TTV with a distinct super-period $P_\text{TTV}^{(c)}$. If the super-periods are well-separated, the signals are orthogonal and individually recoverable.
The sum is distinguishable from a single perturber when: - $|P_\text{TTV}^{(1)} - P_\text{TTV}^{(2)}| / P_\text{TTV}^{(1)} \gtrsim 1$ (non-degenerate super-periods) - The observation baseline covers both super-periods - The SNR per component exceeds 3
It is indistinguishable when: - Both perturbers are near the same MMR (nearly the same super-period) - The total TTV signal looks like a single sinusoid — no harmonic content to distinguish 2 from 1 - In this case, the marginalized posterior on "number of planets" is uninformative
6.2 Second-Order Effects that Break Degeneracies¶
Beyond first order in mass: - Resonant angle libration vs circulation: If the planets are in true resonance ($\phi$ librates), the TTV is bounded and quasi-periodic. If they are merely near-resonance ($\phi$ circulates), the TTV evolves and is technically a chirp. Observable over long baselines. - TTVs of perturber $c$: If $c$ also transits, its TTV is independently constrained, providing mass information on $b$ and breaking the $m_b/m_c$ ratio degeneracy. - Radial velocities: If RV are available, $m_b \sin i_b$ is measured, and combined with TTV $m_b/m_c$, individual masses are determined.
Part II: Astrometric Parallax for NEO Distance Measurement¶
7. Baseline and Distance Relationship¶
For a Near-Earth Object (NEO) at distance $\Delta$ (AU), observed simultaneously from two sites separated by baseline $B$ on Earth, the parallax angle is:
$$\varpi = \frac{B}{\Delta} \quad \text{(in AU, if } B \text{ in AU)}$$
The angular parallax: $$\alpha_\text{parallax} = \frac{B}{\Delta} \text{ radians} = \frac{B_\text{km}}{1.496 \times 10^8 \text{ km}} \cdot \frac{1}{\Delta_\text{AU}} \text{ radians}$$
For $B = D_\oplus / 2 = 6371$ km (half the Earth's diameter — the maximum baseline from two antipodal sites), $\Delta = 0.1$ AU (a close NEO):
$$\alpha_\text{parallax} = \frac{6371}{1.496 \times 10^8 \times 0.1} = 4.26 \times 10^{-4}\,\text{rad} = 87.8\,\text{arcsec}$$
That is an enormous parallax. For $\Delta = 1$ AU (a more typical NEO): ~8.8 arcsec.
NEO parallax is trivially detectable from ground-based astrometry for objects within 1 AU, given arcsecond-level astrometry. [Source: standard parallax measurement; see Bernstein & Khushalani (2000), AJ 120, 3323]
7.1 Astrometric Precision Required¶
To measure distance to $P\%$ relative precision, the parallax must be measured to $P\%$: $$\sigma_\alpha < P\% \times \alpha_\text{parallax}$$
For 10% distance precision on a 0.1 AU NEO: $\sigma_\alpha < 8.8$ arcsec. Trivially achievable (any astrometric measurement with plate-solving to 1-arcsec precision is adequate).
For 1% distance precision: $\sigma_\alpha < 0.88$ arcsec. Achievable with standard astrometry.
For 0.1% distance precision: $\sigma_\alpha < 0.088$ arcsec = 88 mas. Requires careful PSF-fit astrometry using Gaia reference stars. Achievable with a 20 cm telescope and Gaia-calibrated plate solutions (per-frame astrometric precision ~50–100 mas).
7.2 Timing the Simultaneous Observations¶
For the parallax to be measured cleanly, the two observations must be simultaneous (or within a time short enough that the NEO's motion is negligible). For a NEO at 0.1 AU with typical angular velocity:
$$\dot{\theta}\text{NEO} \approx \frac{v\text{NEO}}{D} \approx \frac{30\,\text{km/s}}{0.1 \times 1.496 \times 10^8\,\text{km}} \approx 2 \times 10^{-6}\,\text{rad/s} = 0.41\,\text{arcsec/s}$$
For a 1% distance precision requirement ($\sigma_\alpha = 0.88$ arcsec), the maximum time delay between the two observations is: $$\Delta t_\text{max} = \frac{\sigma_\alpha}{\dot{\theta}} = \frac{0.88}{0.41} \approx 2.1\,\text{s}$$
For simultaneous parallax of a close NEO, the observations must be within a few seconds of each other. This is achievable with NTP synchronization at $\sim$0.1 s level.
7.3 The Science Case¶
Amateur network parallax for NEOs is scientifically useful for: 1. Newly discovered NEOs lacking follow-up orbits — a rapid distance measurement helps confirm/refine the orbit 2. Close-approach events where precise distance determines whether a subsequent close-approach threatens impact 3. For validation of OpenAstro astrometry pipeline — comparing measured parallax to orbital prediction verifies the astrometric precision
The practical procedure: 1. Both sites observe the NEO simultaneously over ~10 minutes 2. Extract positions relative to Gaia reference stars in each image (plate solve) 3. The NEO position in each image differs by $\alpha_\text{parallax} = B / \Delta$ 4. Invert for $\Delta$
This requires absolutely precise knowledge of the two site coordinates and the site baseline vector. Use GPS for both sites; the baseline is then known to meter precision, adequate for arcsecond parallax measurements.
[NOVEL] The proposal to use simultaneous astrometric observations from the OpenAstro distributed network for NEO parallax distance measurement — as both a science product and a pipeline validation tool — is original to OpenAstro.
8. Error Budget for NEO Parallax¶
The uncertainty in the parallax measurement: $$\sigma_\varpi = \sqrt{\sigma_{\alpha_1}^2 + \sigma_{\alpha_2}^2}$$
For matched astrometric precision $\sigma_\alpha$ at each site: $$\sigma_\varpi = \sqrt{2} \sigma_\alpha$$
The uncertainty in distance: $$\sigma_\Delta = \frac{B}{\varpi^2} \sigma_\varpi = \frac{B}{\varpi^2} \sqrt{2} \sigma_\alpha = \frac{\Delta^2}{\sqrt{2} B} \sigma_\alpha \cdot 2$$
Wait, more cleanly: since $\Delta = B/\varpi$: $$\frac{\sigma_\Delta}{\Delta} = \frac{\sigma_\varpi}{\varpi} = \frac{\sqrt{2} \sigma_\alpha}{\alpha_\text{parallax}} = \frac{\sqrt{2} \sigma_\alpha \Delta}{B}$$
For $\sigma_\alpha = 0.1$ arcsec = $5\times10^{-7}$ rad, $B = 5000$ km = $3.3\times10^{-5}$ AU, $\Delta = 0.1$ AU:
$$\frac{\sigma_\Delta}{\Delta} = \frac{\sqrt{2} \times 5\times10^{-7}}{3.3\times10^{-5} / 0.1} = \frac{\sqrt{2} \times 5\times10^{-7}}{3.3\times10^{-4}} = \frac{7.07\times10^{-7}}{3.3\times10^{-4}} = 0.21\% $$
0.2% distance precision with 100 mas astrometry and a 5000 km baseline. Excellent.
9. Summary Tables¶
9.1 TTV Inverse Problem Summary¶
| Quantity | Value | Assumption |
|---|---|---|
| Free parameters (2-planet, coplanar) | 9 | Circular orbits: 5 |
| Required transit sample | $\geq 10$ | For meaningful posterior |
| Required baseline | $\geq P_\text{TTV}/2$ | At least half TTV cycle |
| Compute per MCMC call (TTVFast) | ~2 ms | 3-body, 2-yr baseline |
| Total MCMC time (single CPU) | ~15–30 min | $10^6$ evaluations |
| Minimum detectable TTV (5 mmag, 20 transits) | ~1 min | Hot Jupiter |
| Minimum detectable perturber mass (near 2:1) | ~15 $M_\oplus$ | $ |
| Degeneracy broken by | Chopping + radial velocities |
9.2 NEO Parallax Summary¶
| Scenario | $\Delta$ (AU) | Baseline (km) | Required $\sigma_\alpha$ | Distance precision |
|---|---|---|---|---|
| Close NEO, rough | 0.05 | 3000 | 10 arcsec | 1% |
| Typical NEO | 0.5 | 5000 | 1 arcsec | 1% |
| Distant NEO | 2 | 10000 | 0.5 arcsec | 1% |
| Close NEO, precise | 0.05 | 3000 | 0.1 arcsec | 0.01% |
10. References¶
- Deck et al. (2014), ApJ, 787, 132 — TTVFast algorithm, validation
- Lithwick, Xie & Wu (2012), ApJ, 761, 122 — analytic TTV, chopping signals
- Fabrycky et al. (2012), ApJ, 750, 114 — TTV masses for Kepler multi-planet systems
- Ford et al. (2012), ApJ, 750, 113 — Bayesian analysis of Kepler TTVs
- Buchner (2023), Statistics and Computing — nested sampling review (dynesty)
- Feroz et al. (2019), Open J. Astrophys. — MultiNest algorithm
- Bernstein & Khushalani (2000), AJ, 120, 3323 — NEO orbit fitting from astrometry
- Pravec et al. (2019), Icarus — NEO astrometry with amateur telescopes