Skip to content

Exoplanet transit RV curve

A planetary transit drags the host star's spectral lines back and forth by the Rossiter-McLaughlin amplitude (a few hundred m/s for a Hot Jupiter). With a handful of in-transit spectra plus a template, cross_correlate_rv reconstructs the RV curve in two lines of code.

This walk-through synthesises six exposures across a transit, recovers the velocities with cross-correlation, fits a sine, and reports the amplitude.

1. Synthesise an in-transit time series

Three Gaussian absorption lines, all shifted together at every epoch by an RV following a simple sine model peaking at mid-transit.

import numpy as np
from spectro_kernel.types import Spectrum1D, WorkContext

C = 299792.458  # km/s
true_amplitude_kms = 0.40        # Rossiter-McLaughlin peak
n_epochs = 6
phases = np.linspace(-0.5, 0.5, n_epochs)         # transit centred on phase 0
v_true = true_amplitude_kms * np.sin(np.pi * phases)

REST = (4861.3, 5172.0, 5890.0)                   # Hβ, Mg b, Na D - vacuum-ish

def _gauss(x, c, depth, fwhm):
    sigma = fwhm / 2.3548
    return depth * np.exp(-0.5 * ((x - c) / sigma) ** 2)

wave = np.linspace(4800, 6000, 6000)
spectra = []
for v in v_true:
    flux = np.ones_like(wave)
    for centre in REST:
        flux -= _gauss(wave, centre * (1 + v / C), depth=0.35, fwhm=0.6)
    flux += np.random.default_rng(int(v * 10_000)).normal(0, 0.002, wave.size)
    spectra.append(Spectrum1D(wave, flux))

template = Spectrum1D(
    wave,
    1.0 - sum(_gauss(wave, c, depth=0.35, fwhm=0.6) for c in REST),
)

2. Cross-correlate every epoch against the template

cross_correlate_rv takes the template via ctx.extras['template_spectrum'] and returns the velocity that maximises the CCF.

from spectro_kernel import run_algorithm

recovered = []
for sp in spectra:
    ctx_i = WorkContext(spectrum=sp)
    ctx_i.extras["template_spectrum"] = template
    run_algorithm(
        "cross_correlate_rv",
        ctx_i,
        {"velocity_window_kms": 5.0, "velocity_step_kms": 0.05},
    )
    recovered.append(ctx_i.metrics["radial_velocity_kms"])

for p, vt, vr in zip(phases, v_true, recovered):
    print(f"phase {p:+.2f}  true {vt:+.2f} km/s  recovered {vr:+.2f} km/s")

The recovered values should sit within tens of m/s of the truth - limited mostly by the noise injected above.

3. Fit a sine

We know the model: v(phase) = A * sin(π * phase). Drop into a one-parameter least-squares fit and you have your Rossiter-McLaughlin amplitude.

from scipy.optimize import curve_fit

def model(phase, amplitude):
    return amplitude * np.sin(np.pi * phase)

popt, pcov = curve_fit(model, phases, recovered, p0=[0.3])
amplitude_fit = float(popt[0])
amplitude_err = float(np.sqrt(np.diag(pcov))[0])
print(f"Rossiter-McLaughlin amplitude: {amplitude_fit:.3f} ± {amplitude_err:.3f} km/s "
      f"(injected: {true_amplitude_kms:.3f} km/s)")

4. Render the velocity curve

A standard line plot via the visualisation algorithms or matplotlib - your choice. Quickest is to drop the time series into ctx.light_curves and re-use the periodogram tooling for an out-of-transit baseline.

from spectro_kernel.types import LightCurve

ctx = WorkContext()
ctx.light_curves["rv"] = LightCurve(
    time=phases, flux=np.asarray(recovered), time_unit="phase"
)
# pair with phase_fold + lomb_scargle on a full-orbit dataset.

What you just did

  • Synthesised a Hot-Jupiter style transit.
  • Recovered the in-transit RV with cross_correlate_rv (no per-line bookkeeping, no Doppler maths in user code).
  • Fit a sine to the recovered velocities and extracted the Rossiter-McLaughlin amplitude with realistic uncertainties.

Plug a real CCD frame into the full reduction walkthrough to get your Spectrum1D from raw FITS, then this same notebook is your transit analysis. Five new pages of arXiv per night, with one shared kernel.

Going further

  • Replace the synthetic template by a stellar atlas (Vega, the Sun, BT-Settl)
  • read_fits handles them; cross_correlate_rv does not care where the template comes from.
  • Switch cross_correlate_rv for measure_radial_velocity to compare a CCF estimate with a per-line Gaussian fit. The disagreement diagnoses bad lines.
  • Drop the recovered RVs into a Markov chain Monte-Carlo orbit fitter (out of scope here - but the recovered values are clean enough to feed into emcee directly).