A Christmas Star (Minus the Star, Plus a Moon Nobody Asked For)

Summary

‘Tis the season to stare at light curves. While most people were unwrapping presents on December 25th, I was staring at synthetic flux drop-offs and debugging a limb-darkening model. The result is a browser-based simulator for exoplanet transit photometry, including binary eclipses and exomoon scenarios. It does not detect any real exomoons. It does, however, correctly model why detecting them is comically hard.

Source: sebastianspicker/exoplanet-exomoon-simulation

Background

The gift of transits

When a planet crosses in front of its host star, the observed stellar flux drops by a fraction proportional to the ratio of their projected areas:

For Jupiter transiting the Sun, \( \delta \approx 1\% \). For Earth, about 84 ppm. For an exomoon orbiting a Jupiter-sized planet — well, unwrap that calculation yourself:

A Moon-sized exomoon around a Sun-like star contributes roughly 7 ppm of flux variation. The Kepler space telescope’s photometric precision was on the order of 20–30 ppm per 6-hour cadence for bright targets. Ho ho hold on — that signal is buried.

Why stars are not uniformly bright (and why that ruins everything)

A star is not a flat disk of uniform intensity. It is darker at the limb than at the centre — an effect called limb darkening — because the line of sight through the stellar atmosphere is shallower at the edges, sampling cooler, less emissive layers. The quadratic limb-darkening law is:

where \( \mu = \cos\theta \) is the cosine of the angle from disk centre, and \( u_1, u_2 \) are stellar-type-dependent coefficients. This matters for transit modelling because the depth of the light curve dip changes as the planet traverses from limb to centre to limb — the transit is not a flat-bottomed box, it is a rounded trough. Fitting it incorrectly biases \( R_p / R_\star \) and, more critically for exomoon searches, generates false residuals that look suspiciously like a secondary dip.

Exomoon detection: the indirect approach

No exomoon has been unambiguously confirmed as of the time of writing this post (Christmas Day, 2025 — yes, really). The most promising indirect signatures are:

Transit Timing Variations (TTV). The planet–moon system orbits their common barycentre. This causes the planet’s transit to arrive slightly early or late relative to a pure Keplerian ephemeris. The timing offset scales as:

where \( m_m / m_p \) is the moon-to-planet mass ratio, \( a_m \) is the moon’s semi-major axis around the planet, and \( v_p \) is the planet’s orbital velocity. For an Earth-mass moon at 10 planetary radii around a Jupiter-mass planet at 1 AU, this is on the order of minutes — measurable, in principle, with long baselines.

Transit Duration Variations (TDV). The same barycentre wobble modulates the planet’s velocity along the line of sight at transit ingress, changing transit duration. TDV and TTV are 90° out of phase, which lets you solve for both moon mass and orbital radius given enough transits.

Neither signal is clean in practice. Stellar activity, instrument systematics, and other planets in the system all contribute correlated noise at similar timescales. The residuals of the best exomoon candidate to date — Kepler-1625b-i (Teachey & Kipping, 2018) — remain contested. Season’s readings: disputed.

The Simulation

What it actually does

The simulator is a TypeScript application (Vite build, runs in-browser) built around a deterministic, SI-unit physics core. The main pipeline:

UI parameters
  → V4 normalisation
  → runtime creation  (realtime | reference)
  → orbital integrator (Kepler + N-body)
  → geometry & photometry
  → flux decomposition
  → canvas render + plots

Two runtime modes:

• Realtime — fast integrator, interactive rendering, good for exploration
• Reference — high-fidelity integrator, deterministic export, good for sanity-checking against known systems

The photometry layer computes quadratic limb-darkened transit flux, handles binary eclipse geometry (for eclipsing binary configurations), and exposes hooks for phase curves and instrument noise.

The diagnostics layer is the part I find most useful: energy conservation checks across the integration, radial velocity time series, astrometry, and transit timing outputs. If your N-body integrator is drifting, the energy plot tells you immediately.

The repo ships a real-systems.snapshot.json with versioned data from the NASA Exoplanet Archive — so you can load, e.g., TRAPPIST-1 or HD 209458 as a starting configuration.

What it deliberately does not do

The relativistic corrections are approximations. This is not a GR integrator. For the systems it is designed for (short-period planets around Sun-like stars), the relativistic perihelion precession is tiny — Mercury’s 43 arcseconds per century is the canonical example and that is already a demanding target — but for millisecond pulsars or extremely compact binaries, do not trust it.

The atmospheric module exposes hooks but is not a radiative-transfer solver. If you want realistic transmission spectra, point yourself at something like petitRADTRANS and use this for the orbital geometry only.

Discussion

The simulation is educational in intent — hence the built-in didactic mode (black-box exploration → hypothesis → reveal → A/B comparison → rubric scoring). But the physics is not dumbed down: the limb darkening is real, the N-body integrator tracks multi-body gravitational interactions, and the TTV outputs are computed from first principles rather than parameterised fits.

The thing I kept running into while building this is how much of exomoon detection reduces to a residuals-hunting problem. You fit the best planet-only model you can, examine the timing and duration residuals, and look for a coherent signal. The simulator lets you inject a synthetic exomoon of specified mass and orbital radius, generate synthetic light curves with configurable noise, and see what the residuals look like — which is exactly the kind of intuition-building exercise that is tedious to set up from scratch with, say, a raw BATMAN lightcurve model and a custom integrator.

Limitations worth being honest about. The performance budget is real: some effects are profile-gated to keep the interactive mode responsive, which means the reference mode exists specifically for cases where you want the full physics at the cost of speed. For a publication-quality simulation you would want a dedicated N-body code (REBOUND is the obvious choice), not a browser runtime. This is a tool for understanding the problem, not for writing papers about it — which, fitting for a Christmas project, is exactly what I have time for right now.

References

• Teachey, A. & Kipping, D. M. (2018). Evidence for a large exomoon orbiting Kepler-1625b. Science Advances, 4(10). https://arxiv.org/abs/1810.02362

• Kipping, D. M. (2009). Transit timing effects due to an exomoon. MNRAS, 392(1), 181–189. https://arxiv.org/abs/0810.2243

• Mandel, K. & Agol, E. (2002). Analytic light curves for planetary transit searches. ApJL, 580, L171. https://arxiv.org/abs/astro-ph/0210099

• Claret, A. (2000). A new non-linear limb-darkening law for LTE stellar atmosphere models. A&A, 363, 1081–1190.

Merry Christmas. If you came here expecting warmth and cheer, I offer instead a synthetic light curve with a 7 ppm exomoon signal buried in 30 ppm of photon noise. Practically the same thing.

For the physical version of this — a lamp, a ball, and a smartphone measuring real transit light curves in a classroom — see Hunting Exoplanets with Your Phone. For context on where those experiments came from, see The Lab Goes Home.

Changelog

• 2026-03-05: Corrected the description of the limb-darkening variable from “$\mu = \cos\theta$ is the angle from disk centre” to “$\mu = \cos\theta$ is the cosine of the angle from disk centre.” $\theta$ is the angle; $\mu$ is its cosine.
• 2026-03-05: Corrected Claret (2000) page range from 1081–1090 to 1081–1190. The paper contains extensive tables of limb-darkening coefficients spanning 109 pages.