Analysis of Eclipsing Binary data - Mini-Project

Diego Granziol
25/01/2014
Abstract
I describe my analysis of the eclipsing binary data, taken by Radhakrishnan & Sarma and by Tomkin on the eclipsing binary system R Canis
Majoris. I determine the orbital period as 1.139480 0.0000122 days and
use this to fit the phase radial velocity data for both stars. The orbits
are found to be nearly circular which allows a deduction of the radii of
both stars as 1.12 R (Solar Radii) and 0.8943 R and their masses as
1.002092105 M (Solar Masses) and 0.156081863 M. I also estimate the
effective black body surface temperatures of the two stars as 7393.076K
and 5189.508K. The measurements on the second star suggest it is no
longer on the main sequence and must have had a far greater mass in
the past. A possible scenario for the evolution of the binary system is
discussed in the main text.

Introduction

The future evolution of a star is largely determined by its mass. Although the
situation is changing now, historically stellar masses could only be determined
directly for objects in binary systems and inferred for single stars by comparing
their spectra with those of binaries. The study of binaries, and eclipsing binaries
in particular, remains an important way of determining stellar masses and of
providing insights into the later stages of stellar evolution [ref.Hilditch 2001] The
purpose of the work descibed here is to verify physical and geometric characteristics of the eclipsing binary system R Canis Majoris using previously reported
radial veloctity and light curves [Ref. Radhakrishnan & Sarma, Tomkin]. The
purpose of this experiment is to infer fundamental properties about the eclipsing binary R Canis Majoris from the light and radial velocity curves. We will
then determine paramaters such as the orbital period, the radii, the mass, the
effective temperature and the extent of distortion from a spherical shape for
both stars. I list the experimental objectives here
Determine the orbital period, that is the time taken by both stars to orbit
their centre-of-mass and produce phased light curves for both the Johnson-B
and Johnson-V filters.
Phase radial velocity data for both stars using the derived period.

 Fit sinusoids to both phased radial velocity curves and assess whether or
not the assumption of circular orbits is reasonable.
For circular orbits whose planes lie in the line-of-sight, fitting sinusoids
to the phased radial velocity curves gives orbital speeds; use these to estimate
masses of both stars and their orbital radii.
Deduce the radii of both stars from eclipse timings.
Estimate effective temperatures of both stars from (B-V) colours.
Compare masses, radii and effective temperatures of both stars with those
of the Sun and deduce whether or not these are Main Sequence stars. The
reality of any conclusion reached will depend on the precision with which the
mass-ratio is known.
One of the two stars will be identified as having evolved beyond the Main
Sequence; provide a scenario by which this may have happened.
A detailed fit to the light and radial velocity curve data is needed to
establish the extent to which each star is distorted from its spherical shape; this
shows the extent to which mass-loss and mass- exchange had a role in the earlier
evolution of both stars.
Compare residuals in fits to light curves with those expected, given the
precision of the data, and deduce whether or not there is further information
which could (in principle) be extracted.

1.1
1.1.1

Period Determination
Introduction to the Experiment

The first estimate of the period is achieved by the method of inspection. We are
provided with values of magnitudes of stellar flux from the JohnsonB-band
light

F1
curves in the file B.dat. Armed with the formula x1 x2 = 2.5 log F2 (Where
x1 & x2 refer to magnitudes and F1 & F2 refer to the fluxes) we immediately
 
1
realize that the smaller the value of F1 the greater the value of 2.5 log F
F2 .
Thus in order to identify the minima, i.e smallest values of F1 we need to
identify the greatest magnitudes of x1 x2 . By inspection we find the Times of
the Minima (heliocentric Julian Date of observation) and in the adjacent column
we include the difference in Time between supposed minimas.

1.1.2

Table

Times of Minima (s)

2444606.2928
2444607.4235
2444647.1883
2444648.3210
2444664.2299
2444672.1742
2444998.1999
2444999.3302
2445015.2392
1.1.3

Difference in Time between last Minima

N/A
1.1307
39.7648
1.1327
14.7748
7.9493
326.0207
1.1303
15.909

Calculation

We begin by assuming that our smallest difference in time corresponds to the

period. We then refine our value by finding the closest integer value of n to
which our approximate value divides into a larger Time difference. Then divide
this Time difference by that value of n. We do this step twice and the calculation
is outlined below. Using 1.1303 days as our first estimate 14.7748
1.1303 = 13.072
13. We now can gain a superior estimate of 14.7748
= 1.136523077 repeating
13
39.7648
this procedure for the largest time difference of 39.7648 yields 1.136523077
=
39.7648
34.988 35 thus an even more improved estimate of 35 = 1.136137143.
This is carried out one final time with our largest time difference of 326.0207
326.0207
326.0207
= 1.13596
1.136137143 = 286.9664103 287 giving us a final estimate of
287
to 5dp. The unwary analyst may of course be tempted to skip the first two steps
and proceed immediately to the final step. We note however that 326.0207
1.1303 =
288.4373175 which gives a refined value of 1.132016319 notably different to our
calculate value.
1.1.4

Error analysis

This conclusion can be verified by a simple error calculation. Let us assume

that there exists an exact interval time which corresponds to a complete orbit.
Let us call this time Tp and let us say that the smallest time between which we
register two minima can be written as Tp +e1 , where e1 is the error which we can
model as a symmetric distribution centered around 0. For a larger time between
minima we can write this as nTp + e2 , where n corresponds to the exact number
of revolutions of the secondary around the primary planet. We approximate n
nT +e
as Tpp+e12 = n n, where the delta represents the deviation from the actual
value of n. Now, to remove any ambiguity in the integer value of n, we need
the value of | | to be less than a 1/2. Thus we have
<e2 >+n2 <e21 >2n<e1 e2 >
e22 +n2 e21 2ne1 e2
2
1
< 2 >= T2 p +<e1 2 >+2T
= eT2pne
+e1 =
(Tp +e1 )2
p <e1 >
<e2 >+n2 <e2 >

T2p +<e1 2 >1 41

Where we have assumed that e1 and e2 are independent.

4 < e22 > +(4n2 1) < e21 > Tp2

T 2 +<e2 >4<e2 >

1
2
n2 p 4<e
2>
1
n 100.3427349 100
Where to get an estimate of the max n that would work with this method,
we assumed that our final value calculated was the True value and that our
deviations in the larger times were the root of the mean square error (i.e one STD
DEV away from the zero mean). Given the tiny values of < e21 > & < e22 >the

T2

formula can be well approximated by 4<ep2 > which gives a value of 100.3498233
1
for n (i.e essentially unchanged), much lower than the value of 288 given by the
above naive approach. It essentially comes down to the fact that our repeatedly
refined calculation shows that our best estimates indicate that for the longest
time between minimas 287 orbits as opposed to 288 in the simple calculation
have been completed.

1.2
1.2.1

Period Refinement
Introduction

The next stage was to refine the period time by the string-length method
proposed by Dworetsky (Ref[4]) using the script strphase as described in detail
in section 4.2 of Ref[1]. As proposed in the script we started by selecting the
narrow period of Tp 0.0001 days, increasing in steps of 0.000005, i.e ranging
the value from 1.13586 to 1.13606. The recorded string-length which we wish
to minimize is recorded for the values of interest (i.e surrounding the minimum)
below
1.2.2

String lengths to be minimized and associated periods

String Length
9.694
9.327
9.028
9.142
1.2.3

Period (Days)
1.135935
1.135940
1.135945
1.135950

Further Refinement of the Period

We thus refine our period even further by centering around 1.135945, allowing
strphase to vary in each direction by 0.000001. We get a minimum string length
value of 8.9370.205 for a Tp = 1.1359454, where the value of .205 represents one
STD DEV for the string value, thus we find the values of Tp which give a string
length of one Standard Deviation on either side of our calculated minimum and
take the average. Hence we try to find two values of Tp with a STD DEV of
9.347 and take the midpoint. Giving us
= 1.1359475 days
Timproved = 1.135939+1.135956
2
Timproved = 1.1359475 0.0000088 days

1.2.4

Using the V light Curve to improve the estimate further

We repeat the same procedure with the V data. Finding that we have a minimum at 1.13594443 days which gives a string length of 8.427 .204 and by
repeating the same method of taking the mid point between the values two
standard deviations away we get a value of
= 1.1359485 0.0000085 days
TV lightcurve = 1.135940+1.135957
2
1.2.5

Final Period Result

Averaging the two values and summing the errors in quadrature we get
Tf inal = 1.139480 0.0000122 days
1.2.6

B and V phase light plots

We display the plots of the B and V light curves using Dipso

1.2.7

B light curve (1)

Where the magnitude of the B light curve data is plotted on y against the phase
(measured in multiples of Pi radians) is plotted on x. This plot has been zoomed
in around radians to make for more accurate measurements for calculations
used in Section 2.

1.2.8

B light curve (2)

More general overview of the B Light curve, with phase measured in radians
plotted on the x axis and the magnitude plotted on the y axis.

1.2.9

V light Curve (1)

Identical Plot as in 1.2.7 except this time using Johnson V-band data.

1.2.10

V light curve (2)

Identical Plot as in 1.2.8 except using Johnson V-band data.

1.3
1.3.1

Radial Velocity Curve Fitting

Introduction

We fit functions of the form v = s + a sin(2 + p) (i.e assuming perfectly

circular orbits) to the observed radial velocity curves for both stars using dipso.
For circular orbits we expect the phase difference p to be an integer multiple of
(including 0). Completing the procedure we get two sets of results, one for
the primary star (the larger one) and one for the secondary.
1.3.2

Radial Velocity data

Phase in units of radians

1.3.3

Approximate values of Radial Velocities and Phase offsets

From the graph above we see that for the secondary the difference between
maximum and minimum is around 375 so we get an approximate amplitude
of 187.5 which we use as our first estimate for the secondary radial velocity.
Repeating the the procedure for the Primary we get a radial velocity of about
30. We find that it seems to be centered around roughly -32 so we take that as
our systemic velocity first guess. Whilst the program is not so sensitive to phase
as it is to velocity estimates, we note that if we trace a sinusoidal curve onto
the secondary, that it roughly reaches the 0 point (i.e -32 in our approximation)
because of a design feature in the program not dealing well with an input of 0 we
give it a small non zero value (i.e 0.1). Repeating the procedure for the primary,
we note that tracing a sinusoidal curve through the points and following through
we get an inverted sine curve. Thus we take as a starting point 3.14 as our phase
offset for the primary, the justification for which is provided below.
sin(x + ) = sin(x) cos() + cos(x) sin() = sin(x)
1.3.4

Results from non linear least squares fit

We use the program bnryrvfit to carry out a non-linear least squares fit. We load
the input file binary.dat with the already preset template. We use our previously
obtained estimates for the Systemic Velocity, Radial velocity (for both primary
and secondary), the Phase offset (for both primary and secondary) and the
Period in days (Tp as calculated earlier). Included in the results are also the
10

graphical plots obtained from the fit program. For completeness, and in line
with the method above, we also quote the range which encompasses a 95%
confidence interval around the result. Initially we plot the radial phase velocity
observations for both stars. The secondary is shown in white and the primary
in red. (Fig 1.3.2)
1.3.5

Secondary star results

Systemic Velocity = 32.2270938 kms1

Radial Velocity = 185.390930 0.659569 kms1
95% confidence interval = 183.978858 186.003002 kms1
Phase Offset = 0.0864732908 0.020116 radians
95% confidence interval = 0.0387876481 0.134158933 radians

11

1.3.6

Primary star results

Systemic Velocity = 32.2270938 kms1

Radial Velocity = 28.8757506 0.54489 kms1
95% confidence interval = 27.7390724 30.0124289 kms1
Phase Offset = 3.24332685 0.012718 radians
95% confidence interval = 3.21679608 3.26985762 radians
1.3.7

Comments on Results

It is interesting to note that for the primary star that, even with the 95% confidence interval the phase offset does not encompass the value (3.14..) We
are forced to conclude that the orbit is clearly not circular, but obviously the
ellipse does not have too large an eccentricity (or the value would be much more
distorted from ) so the approximation will still be of use in further calculations. The systemic velocity is not quoted with a standard deviation as it is not
particularly fundamental to our calculations. It will be allowed to vary within
the fitting programs to achieve the best fit and then lowest chi-squared value. It
is also interesting to note that a moments thought would have revealed that for
a circular orbit, the two radial velocities would have always been out of phase
with each other. Assuming we can treat the binary stars as a closed system,
quite clearly by Newtons 3rd law the Centre of mass experiences no net forces
in this closed system and hence cannot move in space with anything other than
12

a uniform velocity from our reference frame. In that frame of reference, given
that there are no net forces, the centre of mass cannot move. Thus by whatever
angle the first star has moved through, the other will have had to move through
that same angle, but in the opposite direction (i.e out of phase).

2
2.1

Light Curve Synthesis

Estimating component masses and orbital radii

v2
1
We find that the mass ratio r = m
m 2 = v1 =
we find the variance given by
2
2 (r)
2 (r)
2 (a1 )
+ a(a2 2 ) = (6.420298214)
2 =
r2 =
a21
2
(r) = 0.1232862899
r = 6.42 0.123

2.1.1

a2
a1

185.39093
28.8875706

(0.659569)2
(185.390930)2

= 6.420298214 and
(.54489)2
(28.8757506)2

Derivation of mass ratio

Assuming circular motion for both stars

m v2
m v2
Gm1 m2
F = r11 1 = r22 2 = (r
2
1 +r2 )
where r1 and r2 specify the distance from stars 1 and 2 to the C.O.M
1 = 2 = (or centre of mass would move)
v = r
m1 r1 2 = m2 r2 2
m1
r2
v2
m2 = r1 = v1
2.1.2

Calculation of the inter-stellar radii

Given that we know the radial velocities (calculated in section 1.3.6) and we
have an estimate for the period of the orbit we can estimate radii with the
v 1 Tp
simple relationship r1 = v1 = 2
and similarly for r2 yielding the results
(1.135948243600)(28.8757806
r1 =
) = 451050956m = 0.648527614 R}
2
r2 = 2895881646m = 4.163740685 R}
r1 + r2 = 4.812268299 R}
Where r1 &r2 signify the distances between stars 1 and 2 to the centre of
mass of the binary system.
2.1.3

Error analysis of the inter-stellar radii

Applying basic calculus of 2 variables and summing in quadrature (i.e assuming

independence
between the two variables) gives us a total error e of

v 2 (T )2 +T 2 (v)2

(where v is either v1 or v2 depending and the sigmas

e=
2
refer to the errors in v and T respectively)
we thus get on inserting the correct values of of v, T, (T ) and (v)
e(r1 ) = 1.8870162%
e(r2 ) = 0.3557736%
13

2.1.4

Calculation of the stellar masses

Manipulating the basic circular motion algebra and remembering the definition
of r as the mass ratio/radius ratio and r1/r2 as the distance from star 1/2 from
the centre of mass, we get expressions for the masses of both stars
r 3 2 r(1+r)2
= 1.993261405 1030 kg = 1.002092105 M }
M1 = 1 G
r 3 2 (1+r 1 )2

M2 = 2 Gr
= 3.104624332 1029 kg = 0.156081863 M }
M1 + M2 = 1.158173968 M }
2.1.5

Error analysis of the stellar masses

The calculation is analogous to that above except for having 3 variables. Adding
in quadrature we get an error for M1 and M2 which is best expressed in percentages
e(M1 ) = 5.29842%
e(M2 ) = 7.65355%
2.1.6

Estimating Effective Radii

From the images 1.2.9 and 1.2.10 we make an estimate of the time in which the
secondary star is blocking the incoming light from the primary and the estimate
the time in which it does this and additionally drops outside of the field of view,
this will be clarified by the algebra below. Where we assume that the stars can
be modeled as disks and that the eclipses are short with respect to the orbital
durations and thus we can model the trajectories as linear.
2(r1 +r2 )

0.2091246306
= 3266.611692
2
v2 +v1 = long = = (
)
1.135948243600

2(r1 r2 )
v2 +v1

= short =
r1 = 1.12 R}
r2 = 0.8943 R}

2.2
2.2.1

0.02338661688
2
( 1.135948243600
)

= 365.308457

Estimating Temperatures
Primary Star Preamble

We measure the B and V magnitude from figures 1.2.8 and 1.2.9 and we calculate
B-V.
V (magnitude) = 0.31
B(M agnitude) = 0.54
B V = 0.23
2.2.2

Determination of Temperature

As given by the regression coefficients of Lynas and Gray in the astronomical

literature [Ref.2] we calculate the Temperature of the star
T = A + Bp + Cp2
where p = B V

14

T = 8515 + (5270)(0.23) + 1440(0.23)2

T1 = 7393.076K
2.2.3

Secondary Star Preamble

Unlike the primary star, assuming a synchronized binary system, the secondary
star will have one hemisphere permanently directed towards the primary. Thus
the observer would see the hot hemisphere immediately before and after the
secondary eclipse (as explained in section 6.2 of the manuscript). Thus the
secondary effective temperature requires a slightly more involved calculation
which is outlined in Appendix B [Ref.1]. Following the calculation and using
the values inferred from the figures 1.2.8 and 1.2.9 for the secondary eclipse we
get
100.4V2 + 100.4V1 = 100.4V12
100.4V2 + 100.4(0.305) = 100.4(0.54)
V2 = 3.02798623
100.4B2 = 100.4B12 100.4B1
B2 = 3.481377662
B2 V2 = 0.453391432
2.2.4

Determination of Temperature

Using a slightly different set of coefficients as we know the mass of the secondary
is much smaller than that of the sun but it still has a roughly solar radius, it
must be on the Red Giant branch, using tabulated values from Lynas and Gray
[Ref.2].
T = A1 + B 1 p + C 1 p 2
T = 4842 + 1413(0.453391432) 1426(0.453391432)2
T2 = 5190K

3
3.0.5

Fitting Light and Radial Velocity Curves with

Nightfall
Editing the V.dat file

Following the instructions laid out in section 6.4 and 6.5 respectively we insert
the following information into the beginning of the V.dat file
#B V
This tells the program that this is the Johnson-V band data, we repeat this
procedure for the B band data replacing this with the instruction #B B
#P 1.13948
This is the refined period time which we found in section 1.2.5
#Z 2444648.3283
#N 0.25
#S 0.0

15

3.0.6

Setting up the configuration file for nightfall

A template file called rcma.cfg will have been copied into the setup. We need
to modify all the entries with the following values (as in type in the numerical
values listed below)
MassRatio r = 6.42 (section 2.1)
Primary Temperature T2 = 5190K(section 2.2.2)
convention of program is to have the opposite convention for Primary and
Secondary as the usual one, thus we mean by this the approximate effective
temperature of the secondary (i.e cooler star)
Secondary Temperature T1 = 7393.076 section (2.2.4)
Absolute Mass M1 + M2 = 1.158173968 section (2.1.4)
Absolute Distance r1 + r2 = 4.812268299 section (2.1.2)
Absolute Period Tf inal = 1.139480 (1.2.5)
3.0.7

Fitting Light and Radial Velocity Curves with Nightfall

We then find that the Chi-Squared is extremely sensitive to changes in s and p.

Where roughly S increases the width of the fit and P the depth. We get a best
fit of
2 = 12.78828
compared to a value of 10 in the manuscript
p = 0.795
s = 0.575
3.0.8
p
0.5
0.6
0.7
0.8
0.7
0.8
3.0.9

Table showing the range of P and S used and the varying ChiSquared
s
0.5
0.5
0.5
0.5
0.6
0.6

2
507.95
205.4
47.37902
176
153
22.33653

Plane of the Orbit Inclined to the line of sight

We try out different values of p and s for different inclination angles. Noting
that the circular motion is not a perfect but only a close fit, and the fact that
if the angle of inclination is much less than 90 degrees we would not observe
any eclipses, we try for 85, 80 and 75 degrees respectively. We also note that
we only have a component radial velocity of sin(i)v where i is the inclination
angle and v is the speed. Thus clearly if i is smaller than 90 we observe only a
portion of the velocity and thus we have underestimated the mass by a factor
of sin3 (i). Thus we divide through by this factor to get the more appropriate
larger mass. The idea is to find a value of i that gives the smallest value of 2
16

and then stick with that value of i and try to find the absolute minimum value
possible.
3.0.10
p
0.83
0.85
0.86
0.87
0.88
0.875
3.0.11
p
0.5
0.6
0.7
0.8
0.7
0.8
3.0.12
p
0.45
0.6
0.7
0.8
0.9
0.9
1
1
1
0.99
0.98
0.98
0.985

85 degrees
s
0.58
0.59
0.6
0.61
0.62
0.65

2
12.82385
11.74469
10.83151
10.53
10.69
10.5113

80 degrees
s
0.5
0.5
0.5
0.5
0.6
0.6

2
507
208
47.39
176.995
153
22.33653

75 degrees
s
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.6
0.7
0.61
0.61
0.605
0.605

2
883.195
638.97
395
189
48.9
43.602
16.67351
30.95197
13.791
12.791
13.39
12.9
12.885

17

3.0.13

Retry of 85 degrees

p
s
2
0.85
0.6
9.854
0.852 0.601
9.7
0.855 0.605
9.32
0.857 0.607
9.185
0.859 0.609
9.07
0.861 0.611
9.09
0.8611 .611 8.41994
We note out of interest that there is not very much difference between the
chi-squared values obtained at 85 degrees and those at 75 degrees with excessive
manipulation of the values s and p, we also note that the experiment expects a
manual adjustment of the parameters which under the time constraints of the
lab time and in general seem rather short sighted, we propose as an improvement
to the experiment that either the student or the demonstrator create a script
that systematically goes through all the values to find an absolute minimum
and is allowed to run overnight to get the job done.
3.0.14

Final results

We increase the chi-squared by 11.5 giving a value of 19.919 and attempt to

find that value by varying each of the parameters (temperature of the secondary
star, the value p and the value s) until it reaches such a Chi-squared. This is
a 99% confidence interval. We give the results as a range of the variable for
the 99% confidence interval, quoting the Chi-Squared at each limit respectively
(which had to extracted manually as well). Gvein that this was a very laborious
process and that the chi squared was incredibly sensitive to initial conditions,
occasionally in the interests of time keeping a Chi-Squared slightly different to
the value of 19.919 was recorded.
3758 T 5611 (2 = 19.91448/19.8977)
0.887 p 0.8214 (2 = 19.91536/19.96457)
0.58624 s 0.65083 (2 = 20.11648/19.95356)

Interpretation

4.0.15

Conclusion

We find that the mass of the primary is M6} with error bars small in comparision.
However given that it has the temperature and radius of the sun we reason that
it can obviously no longer be a main sequence star. Given that the minimum
solar mass multiple required to get onto the main sequences (fusing of hydrogen
to helium) is 0.08 (and we are at about .16). We infer that it used to be a
main sequence star. At some point when the core was no longer able to burn
hydrogen into helium at a point in the core, it expanded into a shell around the
core, pushing the radius of the star outwards until it hit the point of gravitational
18

equilibrium between itself and the next star. At this point there would have
been a transfer of mass (in the form of hydrogen) from one star to another (this
depends on the roche lobe filling factor which should be around 0.99 which says
how much of this volume has been filled) explaining the observed difference in
mass.

References

List
AS33p-2 Analysis of eclipsing binary data - Oxford University Practical Management Scheme
Blackwell DE & Lynas-Gray AE, 1998 AAS 129, 505
Budding E & Butland R, 2011 MN 418, 1764
Dworetsky MM, 1983, MNRAS 203, 917
Hilditch RW, 2001, An Introduction to Close Binary Stars, Cambridge University Press.
Radhakrishnan KR & Sarma MBK, 1982, Contr. Nizamiah & Japal-Rangapur
Obs. No. 16.
Ribas I, Arenou F & Guinan EF, 2002, AJ 123, 2033
Tomkin J, 1985, ApJ 297, 250
Wichmann R, 2002, NIGHTFALL Users Guide

19