Light Curves of Binary Stars

Stars that hide each other

Algol, the second-brightest star in Perseus, loses about two thirds of its light for a few hours every two days and twenty-one hours, and then recovers. In 1783 John Goodricke, a deaf amateur astronomer of eighteen, timed these dimmings precisely and suggested that a darker body passing in front of the star could cause them [1]. He was right. Algol is a pair of stars whose orbit we see almost edge-on, and spectra taken in 1889 showed the bright star moving back and forth in step with the eclipses. Binaries of this kind are more than a curiosity. The masses of stars can be weighed directly only through their orbits, and eclipsing binaries are the main source of the stellar masses and radii on which the theory of stellar structure is tested [3].

The orbit

Two stars of masses $M_1$ and $M_2$ move on ellipses about their common centre of mass, with the same period and the same eccentricity $e$. Their relative orbit has semi-major axis $a$, and Kepler's third law ties it to the period. In units of astronomical units, years and solar masses it reads

and the simulation takes $a$ and the masses as inputs and derives the period, so the two can never disagree. Each star stays on the far side of the centre of mass from the other, on an ellipse scaled by the other star's share of the mass,

Where a star is at a given time follows from Kepler's equation. The mean anomaly $\mathcal M=2\pi(t-T_0)/P$ grows uniformly in time from the moment $T_0$ of closest approach, the eccentric anomaly $E$ solves

and the true anomaly, the actual angle swept around the focus, follows from

To see the orbit from Earth it is turned by the argument of periastron $\omega$ within its own plane and then tilted by the inclination $i$, where $i=90^\circ$ means we see it exactly edge-on. The component $Z$ along our line of sight decides which star is in front,

The light curve

The simulation treats each star as a uniform disc whose surface brightness grows as the fourth power of its temperature, $S\propto T^4$, and leaves out limb darkening and the distortion of stars that nearly touch. When the discs overlap on the sky by an area $A_\cap$, the star in front hides that much of the one behind it, and the total flux is

where $\delta_1=1$ when star 1 is the one behind and $\delta_2=1$ when star 2 is. The curve shows $F$ divided by the uneclipsed total, against orbital phase.

Even this simple model shows how much a light curve gives away, which is what Henry Norris Russell turned into a method in 1912 [2]. For a circular orbit the same area is hidden in both eclipses, so the ratio of the two dips is the ratio of the two surface brightnesses, and with it the ratio of the temperatures,

The deeper dip is always the one in which the hotter star is hidden. The duration of the eclipses measures the sizes. From first to last contact the stars travel a distance $2(R_1+R_2)$ relative to each other, and a flat bottom appears when the smaller star passes entirely in front of or behind the larger one, so the widths of the dips give the two radii in units of the separation. Whether there are eclipses at all depends on the tilt. For a circular orbit the discs touch on the sky only if

which for the default eclipsing binary means an inclination within about twelve degrees of edge-on. An eccentric orbit leaves its own fingerprint. The stars move faster near periastron, so the two eclipses are no longer half an orbit apart, and to first order in the eccentricity the second one is displaced in phase by $(2/\pi)\,e\cos\omega$. The eclipse closer to periastron is also the shorter one.

Inside the code

Kepler's equation has no solution in closed form, so it is solved numerically at every point of the curve. Newton's method converges in a few steps, but for very eccentric orbits it can shoot past the root. The solver keeps a bracket that always contains the solution and falls back to bisection whenever a Newton step would leave it, which makes it reliable up to $e=0.95$, a fact the tests check at every phase.

// Kepler's equation M = E - e sin E has no closed-form solution.
// Newton's method converges fast, but for e close to 1 it can overshoot,
// so every step is kept inside a bracket [lo, hi] that always holds the root.
let lo = m - e, hi = m + e;
let E = e < 0.8 ? m : (m >= 0 ? Math.min(Math.PI, hi) : Math.max(-Math.PI, lo));
for (let i = 0; i < 60; i++) {
  const f = E - e * Math.sin(E) - m;
  if (Math.abs(f) < 1e-12) break;
  if (f > 0) hi = E; else lo = E;
  let next = E - f / (1 - e * Math.cos(E));     // Newton step
  if (!(next > lo && next < hi)) next = 0.5 * (lo + hi);   // fall back to bisection
  E = next;
}

The orbit is then turned and tilted into the sky, and each star placed on its own ellipse.

// From orbital phase to positions on the sky (phase 0 is periastron).
const E = solveKepler(2 * Math.PI * phase, e);
const f = trueAnomaly(E, e);
const r = a * (1 - e * Math.cos(E));
const xo = r * Math.cos(f), yo = r * Math.sin(f);
const xr = xo * Math.cos(w) - yo * Math.sin(w);      // turn by omega
const yr = xo * Math.sin(w) + yo * Math.cos(w);
const X = xr, Y = yr * Math.cos(inc), Z = yr * Math.sin(inc);   // tilt by i
// each star sits on its own ellipse about the centre of mass
star1 = { x: -X * m2 / (m1 + m2), y: -Y * m2 / (m1 + m2), z: -Z * m2 / (m1 + m2) };
star2 = { x:  X * m1 / (m1 + m2), y:  Y * m1 / (m1 + m2), z:  Z * m1 / (m1 + m2) };

The overlap of two discs is the classic lens formula, two circular segments with the limiting cases handled first.

// Area shared by two discs of radii R and r whose centres are d apart:
// two circular segments, the classic "lens" formula.
function overlapArea(R, r, d) {
  if (d >= R + r) return 0;                                    // apart
  if (d <= Math.abs(R - r)) return Math.PI * Math.min(R, r) ** 2;   // one inside the other
  const A = Math.acos((d*d + R*R - r*r) / (2*d*R));
  const B = Math.acos((d*d + r*r - R*R) / (2*d*r));
  return R*R*A + r*r*B - 0.5 * Math.sqrt((-d+R+r) * (d+R-r) * (d-R+r) * (d+R+r));
}

The flux then follows directly, with the star nearer to us hiding part of the other.

// Uniform discs with surface brightness S ∝ T^4. The star nearer to us
// (larger z) hides the overlap area of the one behind it.
const S1 = (T1 / 5778) ** 4, S2 = (T2 / 5778) ** 4;
const overlap = overlapArea(R1, R2, skyDistance);
const hidden1 = star2.z > star1.z ? overlap : 0;
const hidden2 = star2.z > star1.z ? 0 : overlap;
const flux = (S1 * (Math.PI*R1*R1 - hidden1) + S2 * (Math.PI*R2*R2 - hidden2))
           / (S1 * Math.PI*R1*R1 + S2 * Math.PI*R2*R2);

The presets

The eclipsing binary is a pair of Sun-like stars, 1.1 and 0.9 solar masses, on a circular orbit 0.05 astronomical units wide, seen almost edge-on. Its period is 2.89 days, and its unequal temperatures give one deep and one shallow dip. The wide binary is tilted too far for eclipses, and its light stays flat. The eccentric orbit, with $e=0.65$, shows the displaced and unequal eclipses described above, and the near-contact pair shows wide eclipses from stars almost touching. The last preset is Algol itself, with the masses, separation and inclination measured by interferometry [4]: a hot star of 3.17 solar masses and a cooler, larger companion of 0.70 solar masses, 0.062 astronomical units apart. Kepler's third law then gives a period of 2.866 days, against the observed 2.867, and even the simple disc model makes the main eclipse remove 72 percent of the light, close to what Goodricke saw. Algol also hides a paradox. Its less massive star is the more evolved one, which stars of the same age should not allow. The explanation is that the companion was once the heavier star and has since passed much of its mass to its partner.

Try it yourself

Lower the inclination slowly from the edge-on view and watch the eclipses grow shallower and then vanish, close to the angle given by the condition above. Make one star hotter and see which dip deepens. Give the orbit some eccentricity and turn $\omega$, and the second eclipse slides away from phase one half. The sliders for the separation, the radii and the eccentricity mark in red the values at which the two stars would overlap at their closest approach, and a warning appears if you enter that region.

References

1. [1] J. Goodricke, A series of observations on, and a discovery of, the period of the variation of the light of the bright star in the head of Medusa, called Algol, Phil. Trans. R. Soc. Lond. 73, 474 (1783).
2. [2] H. N. Russell, On the determination of the orbital elements of eclipsing variable stars, Astrophys. J. 35, 315 (1912).
3. [3] R. W. Hilditch, An Introduction to Close Binary Stars (Cambridge University Press, 2001).
4. [4] F. Baron et al., Imaging the Algol triple system in the H band with the CHARA interferometer, Astrophys. J. 752, 20 (2012).