Binary Star Systems — two-body orbits, Roche lobes, and simulating them in Three.js

1. What is a binary star system?

A binary star system consists of two stars bound together by their mutual gravity, orbiting a common centre of mass (the barycentre). Surveys of nearby stars suggest that roughly 50–85% of solar-type stars have at least one stellar companion, depending on spectral type — massive O and B stars are almost always multiple, while low-mass red dwarfs are more often single. Our own Sun is unusual in being (as far as we know) a solitary star.

Binaries range enormously in scale. Some pairs orbit with periods of hours, nearly touching; others take hundreds of thousands of years to complete one orbit, separated by a light-year or more. The famous star Alpha Centauri is actually a triple system: Alpha Centauri A and B form a close binary with an 80-year period, while the faint red dwarf Proxima Centauri orbits the pair at a distance of roughly 13 000 AU.

2. The two-body problem and reduction to one body

Two point masses m₁ and m₂ attract each other with Newtonian gravity. Naively this looks like it needs two separate equations of motion, but a classic trick — reducing to the relative coordinate — turns it into a single one-body problem, exactly like a planet orbiting the Sun.

In other words: solve the ordinary Kepler two-body orbit for the separation vector between the stars, using the combined mass M = m₁ + m₂. Then split that relative orbit between the two actual stars in inverse proportion to their masses — the more massive star traces a smaller ellipse closer to the barycentre.

3. Kepler's third law for binaries

For an elliptical relative orbit with semi-major axis a (the sum of the two stars' individual semi-major axes, a = a₁ + a₂), the orbital period follows directly from Newton's generalisation of Kepler's third law:

This single equation is why binaries are so valuable: if you can measure the orbital period P and the physical separation a (from parallax and angular separation), you immediately get the total mass of the system. Combined with the mass ratio from the centre-of-mass motion, both individual masses fall out.

// Kepler's third law solved for total mass, SI-friendly units
const G = 6.674e-11;           // m^3 kg^-1 s^-2

function totalMassFromOrbit(semiMajorAxisM, periodSeconds) {
  return (4 * Math.PI * Math.PI * Math.pow(semiMajorAxisM, 3))
       / (G * periodSeconds * periodSeconds);
}

// Example: Alpha Centauri A+B, a ≈ 23.5 AU, P ≈ 79.9 years
const AU = 1.496e11;
const YEAR = 3.156e7;
const Msun = 1.989e30;
const M = totalMassFromOrbit(23.5 * AU, 79.9 * YEAR);
console.log(M / Msun);  // ≈ 2.0 solar masses combined

4. Centre of mass and mass ratio

Because both stars orbit the shared barycentre, their individual semi-major axes are inversely proportional to their masses. A heavy primary star barely moves; a light companion swings through a much wider orbit:

This is exactly how astronomers detect exoplanets and even black holes indirectly — a visible star "wobbling" around an invisible companion's shared centre of mass reveals the companion's existence through the same equations that describe two visible stars orbiting each other. The Doppler radial-velocity method for finding exoplanets is a direct descendant of binary-star orbit analysis.

5. Roche lobes and mass transfer

Real stars are not points — they have finite size and, in a close binary, can even change shape under their companion's gravity. The Roche lobe is the teardrop-shaped equipotential surface around each star inside which material is gravitationally bound to that star. If a star expands (say, as it evolves off the main sequence and becomes a giant) until it fills its Roche lobe, gas can flow through the inner Lagrange point L1 onto the companion — a process called Roche lobe overflow.

Eggleton's widely used approximation gives the effective radius of the Roche lobe (as a fraction of the orbital separation a) purely as a function of mass ratio q = m₁/m₂:

// Eggleton (1983) Roche lobe radius approximation
function rocheLobeRadius(q, separation) {
  const q23 = Math.pow(q, 2 / 3);
  const q13 = Math.pow(q, 1 / 3);
  const ratio = (0.49 * q23) / (0.6 * q23 + Math.log(1 + q13));
  return ratio * separation;
}

Mass transfer through Roche lobe overflow drives some of the most spectacular phenomena in astrophysics: novae (thermonuclear explosions on an accreting white dwarf's surface), X-ray binaries (accretion onto a neutron star or black hole, heating infalling gas to millions of kelvin), and the Algol paradox — where the less massive star in the Algol system is, counter-intuitively, the more evolved one, because it used to be the heavier star before it dumped most of its mass onto its companion.

6. Classifying binaries: visual, spectroscopic, eclipsing

Astronomers detect and classify binary systems by how they can observe the orbital motion:

• Visual binaries: wide enough apart that a telescope resolves both stars as separate points of light, and their relative motion can be tracked over years (Alpha Centauri A/B is a visual binary).
• Spectroscopic binaries: too close to resolve visually, but the orbital motion periodically Doppler-shifts each star's spectral lines — blueshifted when approaching us, redshifted when receding. Most close binaries are discovered this way.
• Eclipsing binaries: the orbital plane is oriented nearly edge-on to Earth, so each star periodically passes in front of the other, producing a periodic dip in the combined brightness — a light curve. Algol (from Arabic "the ghoul") was the first eclipsing binary recognised, its regular dimming once thought supernatural.
• Astrometric binaries: only one star is visible, but its position wobbles on the sky as it orbits an unseen companion — this is how Sirius B was predicted mathematically decades before it was first observed telescopically in 1862.

7. Simulating an eclipsing light curve

A simplified light curve model treats each star as a uniform disk and computes the overlapping area between the two disks as seen in projection. The combined flux drops whenever one disk overlaps the other:

// Circle-circle overlap area (for simple limb-darkening-free eclipses)
function circleOverlapArea(r1, r2, d) {
  if (d >= r1 + r2) return 0;                 // no overlap
  if (d <= Math.abs(r1 - r2)) {
    const rMin = Math.min(r1, r2);
    return Math.PI * rMin * rMin;             // full eclipse
  }
  const d1 = (d*d - r2*r2 + r1*r1) / (2*d);
  const d2 = d - d1;
  const a1 = r1*r1 * Math.acos(d1 / r1) - d1 * Math.sqrt(r1*r1 - d1*d1);
  const a2 = r2*r2 * Math.acos(d2 / r2) - d2 * Math.sqrt(r2*r2 - d2*d2);
  return a1 + a2;
}

// Combined flux, in units of the two stars' individual luminosities
function combinedFlux(L1, L2, r1, r2, projectedSeparation, primaryInFront) {
  const overlap = circleOverlapArea(r1, r2, projectedSeparation);
  const eclipsedFraction = overlap / (Math.PI * (primaryInFront ? r2 : r1) ** 2);
  const hiddenLum = primaryInFront ? L2 : L1;
  return L1 + L2 - hiddenLum * eclipsedFraction;
}

Two distinct dips appear per orbital cycle: a primary minimum when the fainter, usually cooler star passes in front of (or behind) the brighter one, and a shallower secondary minimum half an orbit later. The relative depth and width of the two dips let astronomers back out both stars' radii and temperatures — this is exactly how the TESS and Kepler space telescopes measure properties of eclipsing binaries (and, with the same technique, transiting exoplanets).

8. Three.js: rendering two orbiting stars

To animate a binary visually, integrate the reduced two-body equation of motion (Section 2) with a simple symplectic integrator — velocity Verlet keeps the orbit stable over many revolutions, unlike naive Euler integration which spirals the orbit outward.

import * as THREE from 'three';

const G = 1.0;          // scene units, not SI
const m1 = 1.6;        // solar masses (primary)
const m2 = 0.9;        // solar masses (companion)
const M  = m1 + m2;

// Relative separation vector state (star2 - star1)
let rel = new THREE.Vector3(6, 0, 0);
let vel = new THREE.Vector3(0, 0, 0.65);  // initial relative velocity

function acceleration(r) {
  const dist = r.length();
  const factor = -G * M / (dist * dist * dist);
  return r.clone().multiplyScalar(factor);
}

function stepVelocityVerlet(dt) {
  const a0 = acceleration(rel);
  rel.add(vel.clone().multiplyScalar(dt).add(
       a0.clone().multiplyScalar(0.5 * dt * dt)));
  const a1 = acceleration(rel);
  vel.add(a0.clone().add(a1).multiplyScalar(0.5 * dt));
}

// Each frame: advance the relative orbit, then split by mass ratio
function animate() {
  requestAnimationFrame(animate);
  stepVelocityVerlet(0.01);

  star1Mesh.position.copy(rel.clone().multiplyScalar(-m2 / M));
  star2Mesh.position.copy(rel.clone().multiplyScalar( m1 / M));

  renderer.render(scene, camera);
}

Trail rendering

A fading orbital trail behind each star makes the ellipse shape immediately visible. Push each frame's position into a fixed-size ring buffer and update a THREE.Line geometry:

const TRAIL_LEN = 400;
const trailPositions = new Float32Array(TRAIL_LEN * 3);
const trailGeo = new THREE.BufferGeometry();
trailGeo.setAttribute('position',
  new THREE.BufferAttribute(trailPositions, 3));
const trail = new THREE.Line(trailGeo,
  new THREE.LineBasicMaterial({ color: 0x60a5fa, transparent: true, opacity: 0.5 }));
scene.add(trail);

let trailIndex = 0;
function pushTrail(pos) {
  trailPositions[trailIndex*3]   = pos.x;
  trailPositions[trailIndex*3+1] = pos.y;
  trailPositions[trailIndex*3+2] = pos.z;
  trailIndex = (trailIndex + 1) % TRAIL_LEN;
  trailGeo.attributes.position.needsUpdate = true;
}

9. Stability, chaos, and circumbinary planets

Two stars alone in an elliptical orbit are perfectly stable forever (in the idealised Newtonian two-body problem — no energy loss, no perturbations). Add a third body — another star or a planet — and the system generally becomes chaotic, as in the classic three-body problem. Planets can still orbit stably around a binary pair if they are either:

• Circumbinary: orbiting far outside both stars, at roughly a_planet > 3–4 × a_binary for stability — like the fictional Tatooine, and the real Kepler-16b.
• Circumstellar (S-type): orbiting tightly around just one star, well inside the gravitational influence of the companion — stable if the companion is far enough away.

Orbits in between these two stable zones are dynamically unstable and get ejected or collide within a few thousand orbits — a direct, testable consequence of the chaotic three-body dynamics that a simple N-body integrator (see the N-body article) can demonstrate by scanning initial planet distances and checking which survive after a long integration.

10. Extensions and improvements

• Eccentric orbits: set the initial relative velocity below or above the circular-orbit value to get elliptical trajectories, and watch the orbital speed change with distance (Kepler's second law — equal areas in equal times).
• Contact binaries: shrink the separation until both stars fill and overlap their Roche lobes, forming a shared, dumbbell-shaped photosphere — render this with a single deformed mesh (metaball-style blending) instead of two spheres.
• Accretion disk: add a flattened particle disk around the mass-gaining star, using the same InstancedMesh particle technique from the spiral galaxy simulation, with an additive-blended bloom pass for a hot inner disk.
• Full three-body scan: numerically integrate a circumbinary test planet at varying initial distances and plot which survive vs. which are ejected — a compact, visually compelling demonstration of orbital chaos.
• Real light curve fitting: extend the light curve model in Section 7 with limb darkening (stars appear dimmer near their edge) for a closer match to real Kepler/TESS data.