A model for eclipsing stars.
Prologue.
This is not a scientific paper. It is my take on eclipsing stars using a simple model with little assumptions. It is about a model with as few (adjustable) parameters as I can get away with. While still adhering to elementary physics.
1. Introduction.
Eclipsing binaries are double stars that pass in front of each other while orbiting. Here we are only concerned with stars that are so close together, and far away from us, that they manifest themselves as one point of light (star). These stars show periodic behaviour as the light diminishes when one star is behind the other.
On the face of it, the light curve is simple: add the light when both stars are seen, otherwise subtract the occluded parts. And indeed such systems exist. When the stars are really close to each other - of the order of a few stellar radii -, other effects start to emerge. The stars heat each other, making a hot hemisphere where the stars face each other. Depending on the time, the hot spots are seen in phases. Like the phases of the Moon or Venus. Also at close distances, tidal attaction will distort the spherical shape of the stars into prolate spheroids. The spheroids present a surface, larger than the undisturbed star when looked at sideways, and smaller when looked at head on.
In this note, the star at the center of the (calculated) orbit will be designated as star 1; the orbiting star is star 2.
2. Stellar Orbit in 3 dimensions.
We start with a stellar orbit model in 3 dimensions. As both stars are merged into one point, we set the semi-major axis to 1 and the rotation of the line of nodes to 0. Those parameters of a stellar orbit model have no effect in an eclipsing star system.
When the orbit has grinded down to a completely circular one, 2 more parameters of the steller orbit model vanish: the eccentricity is 0 by definition for circles, and the longitude from north disappear by want of a periastron. We set the value for the longitude to 0. This latter model has only 3 parameters left: the period, the phase, and the inclination.
The stellar orbit model yields 3 coordinates (x,y,z). The +x coordinate is pointing to the east, -y is pointing north and z is towards the observer. We can transform the rectangular coordinates (x,y,z) to spherical coordinates (ρ,φ,θ). Now ρ is the spatial distance between the stars, φ is the rotational angle in the sky-plane, measured from north, and θ is the angle measured from the xy-plane, up (to the observer) is positive and down is negative. The apparent distance (in the xy-plane) between the stars is indicated as d.
We minimally need 4 more parameters: the radii of both stars and their luminosities. The radii are given in fractions of the semi-major axis, i.e. r1,r2 < 1. The luminosities need to be additive, so they cannot be magnitudes. We need to translate to the magnitude, m, to luminosty, L, and scale them to some nice number range.
 |
(1) |
|---|
In total, these simple eclipsing star models have 6 or 9 adjustable parameters.
3. Constraints
There are 2 obvious constraints for the parameters. Firstly the stars should not crash or rupture each other and secondly they should actually show eclipses.
Constraints are implemented as an extra log prior, added to the log likelihood. It is knowledge we have beforehand, but it can not be implemented with simple, mutually independent probability distributions, acting on the parameters. It is the relation between the parameters that define the prior.
The no-crash constraint entails that the sum of the stellar radii must be smaller than the distance between the stars, or better even the Roche limits must be smaller. The distance between the stars depends on the eccentricity, e, of the system. The closest approach, the periastron distance, is (1 - e) times the semi-major axis, that we have defined as being 1. Here we have a simple prior that is 0 when a crash occurs and 1 otherwise.
The second constraint, actual eclipses, demands that the inclination must be near enough to perpendicular to the sky plane. The constraint is a function of the stellar radii, the eccentricity, the inclination and the angle between periastron and north (longitude). From these parameters we can calculate the nearest approach as projected on the sky plane. We define an extra prior which accepts the orbit setting when eclipses are possible (Pr = 1), rejects them when there can not be any eclipses (Pr = 0) and a linear relation in the in-between zone.
However, some variable double star systems, the so-called heartbeat stars, are not always eclipsing. The heating and distortion effects are the cause of the variability.
4. Overlap.
When the projected distance, d, is larger than the sum or the radii, there is no overlap. When the difference between the radii is smaller than the distance there is complete overlap. In these cases the overlap area, Ao, is defined as
 |
(2) |
|---|
Am is the area of the smaller star.
The overlap of 2 stars approaching each other is given by the pink area in figure 1. The area can be calculated as the sum of the two sectors starting in A and B and subtracting the 2 triangles ACD and BCD.
| Figure 1. Area of two overlapping circles. |
The sector A is found as: 2 β r12. The triangle ACD is equal to AE * ED = r12 cos( β ) sin( β )
For sector B and triangle BCD hold similar equations so that the combined overlap area, Ao, equals
 |
(3) |
|---|
The angle α and β are found using the cosine rule.
From the overlap area, we define 2 visibility functions, one for each star. The sum of these visibilities, multiplied with the luminosities, yield the light curve of the double star.
 |
(4) |
|---|
The Kronecker δ returns 1 when the condition is true, 0 otherwise. It determines which star is in front of the other.
Limb darkening or variations in temperature over the surface of the stars (like sun spots) are not taken into account.
| Figure 2. Eclipsing binary star with a variety of additional settings. |
The light curve of a simple eclipsing star can be seen in figure 2, the black line.
5. Spot Illumination.
When two stars are real close, they heat each other up, on the near sides only. The heating is proportional to the temperature of the other star times it surface, and inversely proportional to the true distance squared. The equation was taken from Dzygunenko and Tvardovskyi or PDF. It gives the extra temperture on star 2, caused by star 1
 |
(5) |
|---|
The hot spot is always facing the other star, so the observer sees
it in phases, with the larger contribution when z < 0, i.e. star 2
is almost behind star 1.
The phase modulates the spot with
cos2( 2 ( θ - 90 ) ).
The proportionality in eq.5 defines another adjustable parameter, fs, which combined with the fact that the luminosity is proportional to T4, yields the result, to be added to L2
 |
(6) |
|---|
There is also a heating of star 1 by star 2, which is the same, mutatis
mutandis, as eq.6, except that the phase is reversed, as star 1 needs to
be behind star 2 for visible full illumination:
phase = ( 1 + z / ρ ).
Even though the hotter star produces a much larger effect on the other than reverse, we implemented both as we don't know in advance which star will be the hotter of the two.
Similar to limb darkening and other temperature variations, we consider the star surface as uniform when doing the overlap calculations. On the other hand, when the stars are eclipsing, the rear one is in "full moon", while the near one is in "new moon". So no terminators in sight.
The effects of spot illumination can be seen in figure 2, the green line.
6. Tidal Distortion.
Another effect present in close binaries is tides due to gravitational differences on the near and far side of a star. The tidal forces follows from Newton's law of gravity:
 |
(7) |
|---|
A similar formula holds for star 2.
The result of this gravitational pull is a distortion of the spherical star into a prolate (elongated) spheroid. We assume that the stars are tidally locked or that the distortion effects act immediately resulting in an elongation directed to the other star.
The ellipticity of a uniform fluid sphere of radius, r2, is given, in first approximation, by equations 1.468 of teaching site --admittedly not a very good reference, but it is all I could find--
 |
(8) |
|---|
In eq.8 we have the ratio of the masses of the stars, m1 / m2. When we normalize the sum of the stellar masses to 1, we can have the (relative) mass of star 1, m1, as extra parameter, from which the (relative) mass of the second star follows as m2 = 1 - m1.
Second order approximations would turn the stars into egg-shapes where the pointy parts are directed to each other, eventually resulting th mass transfer from one star to the other. However that would not be the stable situation we would like to consider.
In the analogous equation for the ellipticity of star 1, the same ratio appears, but now in the inverse. So with one stellar mass as parameter, we can calculate the tidal distortion for both stars.
Assuming that the total volume of the star is preserved, during tidal distortion, we have that the cube of the (nominal) radius, r, equals the semimajor axis, a, times the semiminor axis, b, squared. Together with the ellipticity, ε, which connectd a and b, via b = a ( 1 - ε ), we find
 |
(9) |
|---|
Obviously the ellipticity has to be kept strictly within the range [0,1]. Equation 8, does not automatically guarantee that. We need another constraint on the combination of parameters that yields ε.
Projecting a prolate spheroid, yields an ellipse with the same minor axis as the spheroid, and an apparent major axis, c, varying between both axis. The size of the apparent major axis depends on the aspect angle θ.
 |
(10) |
|---|
In figure 3, we display 6 positions of the secundary star in the binary system. The observer is at the top. The variables a1, b1, and c1 are defined as the true semi-major axis, the semi-minor axis and the apparent semi-major axis of the spheroid, respectively. As the spheroid is rotationally symmetric along the long axis, the apparent semi-minor axis is the same as the true one.
| Figure 3. Gravitational tides. The shape of star 2 changes considerably in its elliptic orbit. On star 1, shape changes are hardly noticeable. |
Assuming, as we did before, that the surface temperature is the same everywhere, we see the luminosity increase when looking sideways at the prolate spheroid and decrease when looking head-on. The luminosity changes proportional to the apparent surface area: c1 * b1 and c2 * b2, respectively.
The effects of tidal distortion can be seen in figure 2, the red and blue lines.
7. Symmetry.
As we are seeing only one dot of light which sums the contribution of both stars, the orbit that fit the light curve, is not unique. Even after we fixed the line of nodes to pointing north and the semi-major axis to 1, there are still several completely identical solutions. We can find the other solutions by mirroring of the main axes.
| Figure 4. Orbital symmetries. |
In figure 4 we display 4 panels. The first one, panel a, shows the starting point, from which we mirror to get to the other panels. In each panel we see the light curve at the top. The numbered arrows refer to the positions of star 2 in the insets low-left and low-right. The insets give two views of the orbit; to the left as eclipsing view and to the right in a sideways view, with the observer further to the right. The effects of spot illumination and tidal distortion is only shown in star 2. It would be confusing to show it also for star 1.
Panel b shows the orbit mirrored in the y-z plane, most clearly in the left inset. We see star 2 pass on the other side over star 1. This can be achieved by changing the inclination into ( 180 - inclination )
Panel c shows the orbit mirrored in the x-z plane, most clearly in the right inset. The stars are running anti-clockwise. This mirroring comes about when both the inclination and the longitude of the periastron are increased by 180 degrees. The orbit is flipped with respect to the sky plane, due to the increase in inclination, and the orbit itself is moved forward by increasing the longitude by half a period.
Panel d shows the orbit mirrored in the x-y plane. What was back is now up front and vice versa. To keep the same light curve we also have to flip the stars. To achieve this we increase the inclination by 180 degrees and exchange the radii and luminosities of the stars, and change the mass parameter into 1 - m1.
The mirrorings are summarized in the table below.
| mirror | panel | incl | long | radius | lumen | mass_1 |
|---|---|---|---|---|---|---|
| y z | b | π-i | ||||
| x z | c | +π | +π | |||
| x y | d | +π | 1<=>2 | 1<=>2 | 1 - m1 |
The 3 mirrorings on fundamental planes, can be combined into 8 parameter sets that all produce the same light curve.
We could make the choice here to allow all these solutions and see where the final ends. However we already have a 10 dimensional parameter space where the solution must be found in a tiny area, with in some dimensions almost no gradient leading to it. If e.g. the period is off by a very small fraction, it is just as bad as when is is off by a large factor.
8. Practicalities.
Constricting the parameter search space as much as possible is a must. First and foremost we need to know the period.
There are 2 problems here. Firstly, we have irregularly spaced data. This prohibits everyones favorite method, the FFT. And secondly, a light curve does not resemble (co)sines in any way. That makes the Lomb-Scargle method much less effective.
However, inspired by the idea that Lomb-Scargle boils essentially down to the fitting of (co)sines to a linear series to frequencies, we replaced the (co)sines for splines. The splines have less of a problem to follow the intrincacies of a light curve, when it has enough knots. Assuming that the eclipse time is about 10 % of the total time, we would need about 20 knots in a periodic configuration. So we at least catch the eclipses at one knot location. The fitting is more complicated than with (co)sines, but it is still a linear problem and easily doable in one (quasi) matrix inversion.
Generally, the eclipses have steep slopes, prohibiting gradients along which to slide toward the true minimum. We need a geometrically spaced grid of about 1000 points per octave, to find a hint of a possible minimum. Subsequently we need to do a downhill search to check and visually inspect whether the true period is found.
Given a period and a splines fit, we can find from the distance between the minima what the minimal value of the eccentricity has to be. The periastron is somewhere between the minima, where the distance is shortest. The luminosities can be found from the depths of the minima and the sign of the inclination from the order of the minima since periastron.
All these values are of course initial guesses, to be fed into the priors for the parameters.
Only with these priors we can hope to find a believable solution.
9. Implementation.
The material in sections 1 to 7, are implemented in EclipsingStarModel which calculates how the light intensity changes as a function of time, and given a set of parameters.
It starts off with a StellarOrbitModel using dummy values for the missing parameters, to produce a 3-dimensional stellar orbit in x, y and z. The distance and orientation of the stars in the sky plane is given by x and y, and z determines which star is in front of the other.
The method tidalDistortion determines how much the stellar dimensions vary in their orbit around each other. The apparent semi-major axes determine when the stars start to overlap.
The method overlap calculates the amount of overlap when the stars pass in front of each other. It assumes that the stars are circular; not exactly true but near enough. At eclipse times, we are looking at the stars almost head-on, when they are truely circular. From the fraction of overlap and the orientation of the stars, we find the visibleFraction.
spotIllumination is a method that calculates the heating of the stars on each other. This results in a brightening of the hemispheres on the near side. It alters the nominal luminosities of the stars. On top of this, the luminosities also change due to change in apparent surface as a result of the tidal distortion.
The luminosities multiplied with the visible fractions sum into the light curve. This is done in the eponymous method lightCurve.
For all these methods there are methods to calculate the derivative to time and partial derivatives to each of the parameters. They were obtained from the Online Derivative Calculator. --Thanks, thanks thanks. A very usefull site.--
Section 8 is implemented in PeriodicScout, still a somewhat experimental class which should only be trusted after careful inspection of the results. Even then, they only provide starting points (priors) for the optimalisation of the final parameters.