Skip to content

Your first SB2 disentangling

disentangle_sb2 separates the two stellar spectra of a double-lined spectroscopic binary from observations taken at several orbital phases. This walk-through builds a tiny synthetic SB2 dataset where the answer is known, so you can sanity-check the algorithm before pointing it at real data.

1. Synthesise two stars in orbit

Two flat-continuum stars with a single Gaussian absorption line each - one at 4860 Å (primary, slow rotator) and one at 4862 Å (secondary, faster). We observe them at six orbital phases, with the two radial velocities anti-phased.

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

C = 299792.458  # km/s

def _gaussian_line(wave, centre, depth, fwhm):
    sigma = fwhm / 2.3548
    return depth * np.exp(-0.5 * ((wave - centre) / sigma) ** 2)

wave = np.linspace(4830, 4890, 2000)
v1_list = [-50, -30, +0, +30, +50, +20]    # primary velocities (km/s)
v2_list = [+50, +30,  0, -30, -50, -20]    # secondary, anti-phased

spectra = []
for v1, v2 in zip(v1_list, v2_list):
    centre1 = 4860.0 * (1 + v1 / C)
    centre2 = 4862.0 * (1 + v2 / C)
    flux = (
        np.ones_like(wave)
        - _gaussian_line(wave, centre1, depth=0.40, fwhm=2.5)
        - _gaussian_line(wave, centre2, depth=0.30, fwhm=1.5)
    )
    flux += np.random.default_rng(int(v1 + 200)).normal(0, 0.003, wave.size)
    spectra.append(Spectrum1D(wave, flux))

ctx = WorkContext(spectra=spectra)

2. Disentangle

from spectro_kernel import run_algorithm

run_algorithm(
    "disentangle_sb2",
    ctx,
    {
        "v1_kms": v1_list,
        "v2_kms": v2_list,
        "n_iter": 80,
        "n_grid": 8000,
    },
)

primary = ctx.extras["primary_spectrum"]
secondary = ctx.extras["secondary_spectrum"]

3. Inspect the recovered spectra

primary is the rest-frame spectrum of star 1 - it should now show only the 4860 Å line, with the 4862 Å line gone (it has been assigned to secondary).

def _line_center(spec, search_lo, search_hi):
    region = spec.select_range(search_lo, search_hi)
    return float(region.wavelength[int(region.flux.argmin())])

print("recovered P:", _line_center(primary,   4855, 4865), "Å  (expected 4860)")
print("recovered S:", _line_center(secondary, 4858, 4866), "Å  (expected 4862)")

You should see both centres within a fraction of an Ångström of the truth.

Caveats

  • disentangle_sb2 is the simplified wavelength-domain method (Simon & Sturm 1994 spirit). It assumes the input velocities are known - feed it values from your orbital solution. If you don't have one yet, run measure_radial_velocity on the strongest line of each epoch first.
  • The accuracy improves with more epochs at well-separated phases, and with spectra that have decent SNR (>30) in the wavelength range you care about.

Going further

  • Run it on a real SB2 (e.g. one of the Spectroscopic Binaries Catalog references). Be honest about your velocity uncertainties - they propagate into the disentangled spectra.
  • Pair with fit_voigt_line once the components are recovered: most SB2 components have pressure-broadened wings.
  • Use cross_correlate_rv on the disentangled primary against a synthetic template to refine the orbital solution iteratively.