A&A 392, 895-907 (2002)
DOI: 10.1051/0004-6361:20020949

Complex Period Variations in the Binary System IM Aurigae

T. Borkovits1,4 - Sz. Csizmadia2 - T. Hegedüs1 - I. B. Bíró1 - Zs. Sándor2,3 - A. Opitz3


1 - Baja Astronomical Observatory of Bács-Kiskun County, 6500 Baja, Szegedi ut, PO Box 766, Hungary
2 - Konkoly Observatory of the Hungarian Academy of Sciences, 1525 Budapest, PO Box 67, Hungary
3 - Department of Astronomy, Eötvös Loránd University, 1117 Budapest, Pázmány P. sétány. 1/A, Hungary
4 - Visiting astronomer at Konkoly Observatory

Received 5 March 2002 / Accepted 29 May 2002

Abstract
Results of a new light-curve analysis and period variation study of the eclipsing binary IM Aur are presented. Four solutions are given according to the different spectral types of the primary, available in the literature. Multicolour photometry and the period analysis both indicate a close ($\approx $ $ 0\hbox{$.\!\!^{\prime\prime}$ }1$) third physical member. If it is considered as a main-sequence star, then its spectral type is A2V-A8V, with a mass of $m_3\approx1.7{-}2.5~M_\odot$. The light-time effect is modelled by the authors' simultaneous secular and periodic terms fitting code. Clear evidence of a secular period change is also revealed. The detailed effects of the perturbations by a third member of a close binary system are extensively studied, using one of the authors' latest numerical integrator code. Although the amplitudes of the expected O-C changes are too small and also the data set is far undersampled for clear confirmation of the theoretical expectation of the third body-perturbed eclipsing binary O-C curve, some typical patterns are shown for their possible identification amongst the O-C curves of other closer systems (like Algol, $\lambda$ Tau) in the near future.

Key words: methods: numerical - celestial mechanics - binaries: close - binaries: eclipsing - stars: individual: IM Aur


1 Introduction

Eclipsing binary systems are widely used as astrophysical laboratories (recently e.g. Guinan et al. 2000). Due to the Vogt-Russell theorem and Kepler's third law all physical properties of their components can be expressed as a function of the orbital period and the mass-ratio in strict cases. Although generally the situation can be more complicated, however, the period is one of the most important parameters of eclipsing binaries.

Nevertheless, the period is influenced by several mechanisms, e.g. by the radial velocity of the mass centre of the system (see Kopal 1959), the apsidal motion (see Claret 1997 for a review), the light-time effect (see e.g. Frieboes-Conde & Herczeg 1973; Borkovits & Hegedüs 1996), mass transfer, and magnetic activity (see Applegate 1992; Kalimeris et al. 1994b; Kaszás et al. 1998, or more recently Zavala et al. 2002, and the criticism of the magnetic activity-period variations connections in van't Veer 1991). In the present paper we shall investigate the properties of the light-time effect (LITE) in the case of IM Aurigae and it will also be shown how the third light affects the light-curve solution in this particular example. Since the fractional number of triple, quadruple or multiple stars is rather high (e.g. the estimated frequency of triple stellar systems is about half of all the multiple stars, see Batten 1973) it is important to subtract the LITE (if it is considerable) from the O-C curve before an analysis of the physical processes in a system. A closer third body can strongly affect the separation and the orbital eccentricity of the inner pair, and this effect changes the form of the components, the shape and volume of the Roche lobes. This change may give rise to intermittent mass transfer (this can be a possible origin of the observed short time-scale period jump[s]).

It should be noted that the real nature of IM Aurigae was misinterpreted several times. Gülmen et al. (1985) (hereafter GSG85) carried out a light-curve solution without taking into account any third light despite the close ($\approx $ $ 0\hbox{$.\!\!^{\prime\prime}$ }1$) companion (which, however, was unknown at that time). Recently, IM Aur has been included in an "apsidal motion catalogue" (Petrova & Orlov 1999) in spite of the fact that both the primary and secondary minima are moving in the same phase (Bartolini & Zoffoli 1986, hereafter BZ86).

The eclipsing nature of IM Aur was discovered by Strohmeier (1959). The correct period of IM Aur was established by Margoni et al. (1966). Kondo (1966) analyzed his own photoelectric light curves using the method of Russel-Merrill. Period variations were first recognized by Gülmen et al. (1984). BZ86 explained this variation as the presence of LITE on the O-C diagram.

IM Aurigae also shows interesting H${\alpha}$-emission (Vesper et al. 2001) and ultraviolet absorption (Bruhweiler et al. 1986) as traces of certain internal activities in the binary.

The aim of this paper is to study the period variations of IM Aurigae and to interpret them mainly in the frame of the interactions among the three components. A new light-curve solution is also presented in Sect. 2. and a detailed period analysis is given in Sect. 3. Furthermore, numerical studies of the perturbations in the orbital elements are presented in Sect. 4.

2 New light-curve solution

GSG85 published a BV light curve obtained with a single-channel photomultiplier-based photometer. They gave a solution of their own measurements using the Wood and Wilson-Devinney methods. However, there are several circumstances requiring a new light-curve solution.

These are as follows:

(i) The primary component of IM Aur was taken to be a B9 star according to Mammano et al. (1967) and that is why its surface temperature was fixed at 10 500 K in GSG85. However, other studies give an earlier spectral type for the primary (e.g. Roman 1978 gives B6V, Bruhweiler et al. 1986 give B5V or B6V). In our study we accepted B7V from the SIMBAD database for the main component, which yields a higher surface temperature ( T1 = 13 500 K).

(ii) The presence of the third companion was not considered in the former analysis. Nevertheless, according to the LITE solution of BZ86 as well as our new one (see below), this star might produce extra radiation increasing the luminosity of the system. Disregarding this effect would necessarily modify the resulting physical parameters of the eclipsing pair.

Accordingly, we decided to reanalyse the original observational data of GSG85. As an initial parameter set, the solution of GSG85 was applied. We used the Wilson-Devinney (WD) Code (Wilson 1998). Instead of the normal points (which were used in GSG85) the original individual points were used.

The re-analysis of the light curve of IM Aur is difficult due to uncertainty of several factors: the spectral type of the primary, the spectroscopic features of the secondary and the mass ratio. GSG85 obtained a mass ratio q=0.32 employing the single-lined radial velocity curve obtained by Mammano et al. (1967). Since this radial velocity curve is very undersampled and exhibits large scatter, we did not use it directly in our investigation. (Note that each spectral exposure covered almost 0.04 phase.) Instead, several different solutions were obtained. First we fixed the mass ratio at q=0.32 as found in GSG85. The following parameters were adjusted: i, $\Omega_{1,2}$, T2, fractional luminosity of the primary (L1), third light (l3), phase shift ( $\phi_{\rm sh}$). When the adjustments became smaller than their errors a solution was gained. We found that the third light is negative in detached mode. Running the WD Code in Mode 5 (semi-detached configuration) yielded a better fit. A relatively high value of l3 was obtained in this mode (approximately 9%). The results of the solution are presented in Col. "Model 1'' of Table 1. As a next step of our analysis we checked the effect of changing the uncertain parameters q and T1 as well. Throughout Model 2 we accepted continuously the B7V spectral type for the primary, and the mass ratio was adjusted. A significantly larger value was found for the mass ratio as well as for the contribution of the third light (see Col. "Model 2'' in Table 1). As a last step, two further solutions were obtained for the system, assuming B5V (Model 3) and B9V (Model 4) spectral types respectively for the primary. In these models the mass ratio was also fitted. The results are shown in the last two columns of Table 1. As one can see these latter three light-curve solutions do not differ from each other significantly.

The graphical fits are shown in Figs. 1 and 2, respectively, while the calculated absolute dimensions can be found in Table 2. These were calculated in the following way. The mass of the primary was estimated from its spectral type via Lang's tables (Lang 1992). Then the semi-major axes for the four different solutions were calculated using Kepler's third law. Finally the "lc'' program of Wilson-Devinney Code was used for estimating the absolute dimensions (see also Wilson 1998).

The major difference of the new solutions from that of GSG85 is that we obtained a semi-detached configuration. (Mardirossian et al. 1980 obtained the same result by reanalyzing the light curve of Kondo 1966.) Furthermore, the amount of the third light in our solutions is in accordance with the value expected from the variation of the moments of the eclipses (see the next section).


 

 
Table 1: New light-curve solutions of IM Aurigae.

Parameters		 Model 1 Model 2 Model 3 Model 4

i		 78 $.\!\!^\circ$109(305) 77 $.\!\!^\circ$068(074) 76 $.\!\!^\circ$769(074) 76 $.\!\!^\circ$878(075)
$\Omega_1$ 3.2699(89) 3.5194(114) 3.5738(116) 3.5177(117)
$\Omega_2$ 2.5100 2.7407 2.7419 2.8021
$\left(\frac{L_1}{L_1+L_2+L_3}\right)_B$ 0.8817(290) 0.8157(179) 0.8159(230) 0.8199(298)
$\left(\frac{L_1}{L_1+L_2+L_3}\right)_V$ 0.8323(154) 0.7693(169) 0.7686(218) 0.7737(282)
$\left(\frac{L_3}{L_1+L_2+L_3}\right)_B$ 0.0819(142) 0.1366(30) 0.1398(36) 0.1369(51)
$\left(\frac{L_3}{L_1+L_2+L_3}\right)_V$ 0.1100(134) 0.1568(35) 0.1590(90) 0.1591(57)
T1 ( K) 13 500 13 500 10 500 15 400
T2 ( K) 6 259(28) 6 487(24) 5 558(19) 6 848(30)
F1 1.00 1.00 1.00 1.00
F2 1.00 1.00 1.00 1.00
g1 1.00 1.00 1.00 1.00
g2 0.32 0.32 0.32 0.32
A1 1.00 1.00 1.00 1.00
A2 0.50 0.50 0.50 0.50
q 0.320 0.431(4) 0.431(4) 0.462(5)
x1,B 0.389 0.389 0.463 0.357
x1,V 0.333 0.333 0.395 0.304
x2,B 0.678 0.678 0.782 0.627
x2,V 0.547 0.547 0.641 0.511
$\sum wr^2$ $4.53 \times 10^{-3}$ $4.27 \times 10^{-3}$ $4.50 \times 10^{-3}$ $4.20 \times 10^{-3}$



  \begin{figure} \par\includegraphics[width=13cm,clip]{3529f1.eps}\end{figure} Figure 1: New light-curve solutions of IM Aurigae in V light. Crosses denote the observed points of GSG85. Continuous line stands for Model 1 solutions, while long dashed, short dashed, and dotted ones represent Models 2, 3, and 4, respectively.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=13cm,clip]{3529f2.eps}\end{figure} Figure 2: New light-curve solutions of IM Aurigae in B light. Crosses denote the observed points of GSG85. Continuous line stands for Model 1 solutions, while long dashed, short dashed, and dotted ones represent Models 2, 3, and 4, respectively.
Open with DEXTER

From the light-curve solution in different colours, the colour index of the tertiary component can be determined from the colour index of e.g. the primary by the following expression:

\begin{displaymath}(B-V)_3 = -2.5\left[\log\left(\frac{L_{3,B}}{L_{3,V}}\right)-\log \left(\frac{L_{1,B}}{L_{1,V}} \right) \right] + (B-V)_1, \end{displaymath} (1)

where (B-V)i refers to the ith component's colour index, $L_3\approx4\pi l_3$ where l3 is the third light obtained from the light-curve solution, and the luminosities of the components are given in arbitrary units as used by the Wilson-Devinney Code. Getting the (B-V)3 colour index of the tertiary from the above formula for the four models the spectral type and the mass of the third star can be determined e.g. from the tables of Lang (1992). These parameters of the tertiary are listed in the last three rows of Table 2. (We have arbitrarily assumed that the third companion is a main-sequence star.)


  
Table 2: Absolute dimensions of IM Aur.
\begin{table} \begin{displaymath} \begin{array}{lllll} \hline\hline \noalign{\sm... ...dot) & 1.7 & 2.2 & 2.0 & 2.54 \\ \hline \end{array}\end{displaymath}\end{table}


  
Table 3: Times of photoelectric minima of IM Aurigae (HJD -2 400 000). References: (1) Kondo (1966), (2) BZ86, (3) Margoni et al. (1966), (4) Rossati (1967), (5) Dworak (1974), (6) Dworak (1976), (7) Güdür et al. (1982), (8) Gülmen et al. (1984), (9) Diethelm (1982), (10) Isles (1985), (11) Bruhweiler et al. (1986), (12) Diethelm (1985), (13) Isles (1988), (14) Isles (1990), (15) Mayer et al. (1991), (16) Hübscher & Lichtenknecker (1988), (17) Isles (1991), (18) MitVS (1989), (19) Hübscher et al. (1990), (20) Hübscher et al. (1991), (21) Diethelm (1992), (22) Hübscher et al. (1993), (23) Diethelm (1995), (24) Agerer & Hübscher (1996), (25) Agerer & Hübscher (1999), (26) unpublished, (27) Csizmadia et al. (2001), (28) Bíró & Borkovits (2000), (29) Borkovits et al. (2001), (30) present paper.
\begin{table} \begin{displaymath} \begin{array}{llll} \begin{array}{lll} \hline\... ...81 &{\rm s}&30\\ \hline \end{array} \end{array} \end{displaymath} \end{table}
$^{a,\ b}$ These minima times were omitted from the calculations. (See text for details.)

3 Period variations

All photoelectric minima found in the literature are listed in Table 3. Some of them were omitted from our calculations for various reasons. If we had strongly different times for the same or very close eclipse events, we decided to exclude those with extreme residuals. Other minima were omitted as well despite the lack of other data in their vicinity, due to their very deviating values. We do not think that all omitted minima are due to observational errors. For example, one of us (Sz.Cs.) observed an unexpected sudden fading of the system (this event is marked with superscript b in Table 3). But even if some of the cancelled points may represent real values, they probably refer to other processes which are out of the scope of our present study.

The time that has elapsed from the discovery of the light-time effect of IM Aur by BZ86 to the present is almost the same as the time span between the first published photoelectric minima by Kondo (1966) and their work. Now we are in the position to realize that the O-C curve shows rather complex behaviour (see Fig. 3). A quasi-sinusoidal curve is superimposed on another kind of period variation. At first glance, this latter feature could be either a continuous secular period decrease or an abrupt period change. Due to the short time interval over which the times of the minima are distributed and the relatively large amplitude of the periodic variation we are not yet able to choose among these interpretations. Despite this fact, we applied the usual quadratic description for modelling the long-term behaviour of our curve. (As is well known, if the period of a system varies by a constant rate, this yields a quadratic O-C curve.) In the last decade other mathematical methods have been developed for the analysis of O-C curves as well (see e.g. Kalimeris et al. 1994a; Jetsu et al. 1997). The main reason why we did not apply any of these methods is that we used "a priori'' knowledge on the trace of the LITE on the O-C diagram. The analysis of the O-C curve was done in the same manner as in Borkovits & Hegedüs (1996) with one fundamental difference, namely that the parabola was not subtracted previously, instead its coefficients were fitted simultaneously with the Fourier-coefficients of the periodic terms. Thus, the shape of our diagram was modelled by the following form:

 \begin{displaymath}f=c_0+c_1E+c_2E^2+\sum_{j=1}^3(a_j\sin{j\omega E}+b_j\cos{j\omega E}), \end{displaymath} (2)

where the least squares fitting method has been carried out with different (fixed) frequencies ($\omega$). The coefficients of the best fit are listed in Table 4. The upper panel of Fig. 4 represents the above described best fit, while the lower one shows the corresponding "instantaneous'' period (the epoch-derivative of the fitted curve). The extra humps on the descending branches are due to the overtones arising from two causes. The more important one is due to the large gaps in the data set, whose shortness results in minor deviations from the theoretical amplitude ratio following from the Fourier-expansion of the Keplerian motion (see e.g. Kopal 1978). The other reason is the use of only the first few terms of an infinite expansion.
  
Table 4: The Fourier coefficients, and their errors from best LSQ fit. In the last two rows the $\chi ^2$ can be seen for the LSQ curve, and the calculated LITE curve respectively.
\begin{table} \begin{displaymath} \begin{array}{lrlr} \hline\hline \noalign{\sma... ...^2_{\rm LITE}&0.0011270&&\\ [1mm] \hline \end{array}\end{displaymath}\end{table}

The interpretation of the periodic terms in the O-C variation as LITE seems to be less problematic. We examine this case first.

3.1 Light-time effect

The orbital elements of the close binary in the triple system, as well as the mass funtion of the tertiary component can be read in Table 5.

  
Table 5: The LITE solution of this paper compared with the earlier results. The orbital elements refer to the wide orbit of the eclipsing binary in the triple system.
\begin{table} \begin{displaymath} \begin{array}{lllll} \hline\hline &&{\rm {pres... ...&({M}_\odot)&0.084&0.015&0.124\\ \hline \end{array}\end{displaymath}\end{table}

These values are very close to the results of BZ86. This is due to the fact that the available data were spread over a shorter time interval, thus the secular trend was almost unobservable. (Compare the left half of our Fig. 5 with their Fig. 1!) Unfortunately, the last extremum of the LITE curve (which almost corresponds to the pericenter passage of the wide system) occurred last summer when IM Aurigae was not observable from our geographical longitude. As the sharper extrema were very poorly observed also in the past (see Fig. 6), this yields a source of uncertainty, mainly in the projected semi-major axis of the system (through the amplitude of the LITE), and as a consequence, in the mass-function value of the third companion as well.

As is well known, the unknown mass m3 (as a function of the inclination) can be derived from the following third-order equation:

 \begin{displaymath}m_3^3\sin^3{i'} - m_3^2f(m_3) -2m_3m_{12}f(m_3) -m_{12}^2f(m_3) = 0, \end{displaymath} (3)

where m12 and m3 are the masses of the eclipsing pair and the third body, respectively, while the mass function f(m3) can be expressed by the elements of the LITE solution as

 \begin{displaymath}f(m_3)=\frac{4\pi^2a^{'3}\sin^3{i'}}{GP'^2}, \end{displaymath} (4)

where G is the gravitational constant. The calculated masses of the tertiary for different inclinations are listed in Table 6.


  
Table 6: The mass of the third body for different visible inclinations.
\begin{table} \begin{displaymath} \begin{array}{lllll} \hline\hline \noalign{\sm... ...35)&4.25(36)&3.82(32)&4.94(41)\\ \hline \end{array}\end{displaymath}\end{table}

Although BZ86 assumed that the wide and close orbits are coplanar (the same assumption was used by Mayer 1990) there is no physical (or statistical) reason to expect this. Instead, we can use an indirect determination of the observable inclination of the wide orbit. To do so, we assume that the source of the third light, found in our light curve solution, is identical with the source of the LITE. Comparing the third-body masses determined from the mass function with those calculated from the amount of the third light (see the last row in Table 2), we can conclude that according to Model 1 the third body's orbit is seen almost edge-on, but the other models with adjusted mass ratio give an $i'\approx55\hbox{$^\circ$ }{-}58\hbox{$^\circ$ }$ observable inclination for this orbital plane. (It must be stated again that these conclusions are valid only if the tertiary is a main-sequence star.)

Two important points must be emphasized here. First, due to the large uncertainties in the determination of the third mass using an imperfectly covered O-C curve, as well as other inaccurate parameters, and also the hypothesis in the determination of the third mass from the light curve, we do not think that our result would disclose the exact coplanarity of the two orbital planes. Second, since we have no information about the direction of the nodes of the orbits, from the knowledge of the visible inclinations of the orbits we cannot tell anything about the mutual inclination of the two orbital planes.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{3529F3.eps}\end{figure} Figure 3: The O-C plot of IM Aurigae. Ephemeris from BZ86. (Circles denote the primary minima, while rectangles stand for the secondary ones.)
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{3529F4.eps}\end{figure} Figure 4: The simultaneous LSQ parabola and Fourier fit (above), and its derivative, the instantanous "period'' (below).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{3529F5.eps}\end{figure} Figure 5: The LITE solution superimposed to the parabola (above), and the residuals (below).
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{3529F6.eps}\end{figure} Figure 6: The parabola-subtracted O-C and the LITE solution plotted vs. the orbital phase of the tertiary. (The ephemeris is HJD $2~452~085+1372^{\rm d}E$.)
Open with DEXTER

3.2 Secular variation

The following quadratic ephemeris was yielded by the simultaneous least-squares (LSQ) fit of a second order polynomial and Fourier-terms containing a (fixed) fundamental frequency and its first two harmonics:

 
$\displaystyle {\rm {MIN}}_{\rm {I}}=2438327.7970(3)+1.2472899(8)E -1.73(12)\times10^{-10}E^2.$     (5)

Interpreting this period change we examine first the possibility of an eventual mass-flow, and then the rarely studied question: whether the perturbations of the third companion could cause such kind of variation?

3.2.1 Mass flow in (from) the system

Several Algol-systems with lower mass evolved secondaries show a period decrease (see e.g. Qian 2001, and references therein). This phenomenon suggests some angular-momentum loss event, e.g. mass flow from the system, magnetic, or tidal breaking. As it was mentioned in the Introduction, some mass transfer in the system evidently exists.

If the mass-loss were the only source of the observed period variation we would obtain the following estimation:

\begin{displaymath}\dot{m}_{12}=-m_{12}\frac{\dot{P}}{2P}\approx-m_{12}\frac{c_2}{c_1^2}=2.43\times10^{-7}~{M}_\odot ~{\rm y}^{-1} \end{displaymath} (6)

(see e.g. Tout & Hall 1991). (Here the total mass of the binary (m12) was calculated from Model 1, but the results obtained using m12 values coming from the other light-curve models do not differ significantly.) However, this value is 3-4 orders larger than the usual annual rate of wind-driven mass loss of a typical B5-7V star. Consequently the main source of the mass flow must be the Roche-lobe filling secondary. Due to the convective envelope of the secondary some magnetic effects might be accounted for. Nevertheless, according to Tout & Hall (1991), the angular momentum loss via a magnetic stellar wind is effective only in the case of red giant components. So we have to look for some other interpretation.

3.2.2 Perturbations by the third body

The variations of the orbital elements of a double stellar system due to a distant third body have been studied by several authors. (Recently Ford et al. 2000 gave an octuple-level secular perturbation theory.) The typical periods of the different classes of periodic perturbations are P, P', P'2/P (see e.g. Söderhjelm 1975). As is well-known from the textbooks on celestial mechanics, the longer the period of a perturbation, the larger its amplitude. It is clear that only the third type of these perturbations (the so-called "apse-node'' terms) is interesting for us. As long as the orbit of the close binary is (almost) circular, the effect of the apsidal motion can be neglected, so we can concentrate on the nodal motion. Its typical amplitude (in the O-C curve) is

\begin{displaymath}A_{\rm n}=\frac{P}{2\pi}i_1\cot{I}, \end{displaymath} (7)

where i1 is the inclination of the eclipsing binary to the invariable plane, and I is the observable inclination of the invariable plane (Luyten 1938). The ratio of the angular momenta of the two orbits is

\begin{displaymath}\frac{G_2}{G_1}=\frac{m_{12}m_3}{m_1m_2}\left(\frac{m_{12}}{m... ...ac{P'}{P}\right)^{1/3}\left(\frac{1-e'^2}{1-e^2}\right)^{1/2}, \end{displaymath} (8)

which can be estimated in the present case (both for Model 1 and Model 2) as

\begin{displaymath}\frac{G_2}{G_1}\approx12. \end{displaymath} (9)

From this fact it follows that $I\approx i'$. (In the above expression m123 denotes the total mass of the triple system.) According to the third light, if the tertiary is not a degenerate object the observable inclination of the wide orbit (i') must be larger than  $55\hbox{$^\circ$ }$, so

 \begin{displaymath}A_{\rm n}\leq0\hbox{$.\!\!^{\rm d}$ }14i_1. \end{displaymath} (10)

The approximate period of this nodal regression (see again Söderhjelm 1975) is

\begin{displaymath}P_{\rm n}=\frac{4}{3}\left(1+\frac{m_{12}}{m_3}\right)\frac{P... ...}(1-e'^2)^{3/2}\left(\frac{C}{G_2}\cos{i_{\rm m}}\right)^{-1}, \end{displaymath} (11)

where $C=\sqrt(G_1^2+G_2^2)$ is the total angular momentum. In the case of IM Aurigae

 \begin{displaymath}P_{\rm n}\approx 1.^{\!\!{{\rm a}}}3\times10^4(\cos{i_{\rm m}})^{-1}, \end{displaymath} (12)

using Model 1, and almost the same in the other cases. This (minimal) period is longer by almost $50\%$ compared to the value given by Mayer (1990) due to the smaller mass function and larger masses for the primaries. (It should be noted that according to the numerical tests of Söderhjelm 1982 this formula cannot be used in the high inclination regime, e.g. when $i_{\rm m}>i_{{\rm crit}}\approx40\hbox{$^\circ$ }$.) According to the quadratic term, the amplitude of the O-C in the last four decades is about $0\hbox{$.\!\!^{\rm d}$ }02$ which evidently discloses the long period nodal regression as a possible source of this feature of our curve. Furthermore, as is well known from the celestial mechanical textbooks, the rate of the nodal motion can be written as

 \begin{displaymath}\dot\Omega=\frac{f_{\rm n}}{c\sin{i}}\rho_1\sin{u}, \end{displaymath} (13)

where c is the specific orbital angular momentum, $f_{\rm n}$ is the normal component of the force (e.g. $f_{\rm n}=\frac{\vec{c}\cdot\vec{\ddot\rho}_1}{c}$), $\rho_1$ is the separation, u is the true orbital longitude measured from the node, and i is the observable inclination of the binary. It is evident from (13) that the smaller the inclination the faster the nodal motion, and consequently, the largest variation of the node occurs when the system no longer exhibits eclipses.

In the next section we study the perturbed motion of the binary in details.

4 Numerical integration of the motion

In the previous subsection we treated the perturbations arising from the third companion in the mass-point model (stellar three body problem). Nevertheless, this approximation loses its validity when we take into account the nonspherical shapes of the members of the close subsystem due to the tidal and rotational distortion, arising from the strong gravitational interaction between them. To study these additional effects in the orbital evolution of such a triple system we integrated numerically the equations of motion including the tidal and centrifugal terms into the two-body force.

4.1 The equations of motion

The equations of the orbital motion of the three bodies written in the usual Jacobian coordinates have the following form:

  
$\displaystyle \vec{\ddot\rho}_1$ = $\displaystyle -\frac{Gm_{12}}{\rho_1^3}\vec{\rho}_1+Gm_3\left(\frac{\vec{r}_{23... ...frac{m_{3-j}}{m_j}k_j^{(i)}\left(\frac{R_i}{\rho_1}\right)^{2j+1}\right.\right.$  
    $\displaystyle \left.\left.+\frac{k_2^{(i)}R_i^5}{Gm_i}\left(\frac{\omega_{z'_i}... ...rac{\vec{\omega}_{z'_i}\cdot\vec{\rho}_1}{\rho_1^2}\vec{\omega}_{z'_i}\right\},$ (14)
$\displaystyle \vec{\ddot\rho}_2$ = $\displaystyle -\frac{Gm_{123}}{M_{12}}\left(\frac{m_1}{r_{13}^3}\vec{r}_{13}+\frac{m_2}{r_{23}^3}\vec{r}_{23}\right),$ (15)

where kj(i) is the jth apsidal motion constant for the ith star, Ri is its characteristic radius. Furthermore, $\vec{\omega}_{z'_i}$ stands for the uni-axial rotational angular velocity vector of the current star, while the other symbols have their usual meanings. Equation (14) is coupled with the Eulerian equations of rotation of the two members of the close binary via the centrifugal terms. To examine the possible precession and nutation caused by the third body we integrated (simultaneously) these equations as well. Using the formulae of Tokis (1974) (see also Kopal 1978) the final forms of our differential equations are as follows:
 
$\displaystyle \left[m_i{\cal {R}}_i^2+\frac{\alpha_i}{3}(7K^{(i)}_2+8G^{(i)}_2)... ...ec{\dot\omega}_i}{\omega_{z'_i}^2}\vec{\omega}_{z'_i}\right)=\alpha_i\vec{R}_i,$     (16)

where
 
$\displaystyle \vec{R}_i$ = $\displaystyle K^{(i)}_2 \left[\frac{\vec{\rho}_1\cdot(2\vec{\omega}^*+7\vec{\om... ...\cdot\vec{\omega}_{z'_i}}{\rho_1^2}\vec{\rho}_1\times\vec{\omega}_{z'_i}\right]$  
    $\displaystyle - \left\{K^{(i)}_2 \left\{ \frac{\vec{\rho}_1\cdot[(\vec{\omega}_... ...)}_2\vec{\omega}_i\times\vec{\omega}_{z'_i}-\frac{2}{3}H^{(i)}_2\vec{\omega^*},$ (17)

and ${\cal{R}}_i$ is the gyroscopic radius approximated by the formula of Motz (1952), as

\begin{displaymath}\log{\cal{R}}_i=0.24\log k^{(i)}_2-0.13+\log R_i, \end{displaymath} (18)

and

\begin{displaymath}\alpha_i=\frac{1}{4}k^{(i)}_2m_iR^5_i. \end{displaymath} (19)

Furthermore, K(i)2, and G(i)2 come from the amplitude of the tidal and rotational distortions, namely
K(i)2 = $\displaystyle (1+2k_2^{(i)})\frac{m_{3-i}}{m_i}\frac{1}{\rho_1^3},$ (20)
G(i)2 = $\displaystyle -\frac{\omega_{z'_i}^2}{3Gm_i},$ (21)

while

\begin{displaymath}H_2^{(i)}=\dot{K}_2^{(i)}=-3K^{(i)}_2\frac{\vec{\rho}_1\cdot\vec{\dot\rho}_1}{\rho_1^2}\cdot \end{displaymath} (22)

The angular velocity vectors $\vec{\omega}_i$ represent the rotation of the coordinate system co-rotating with the ith star with respect to the inertial system, while $\vec{\omega}^*$ stands for the angular velocity vector of the coordinate system co-orbiting with the binary. (A preliminary description of the integrator can be found in Borkovits 2001, while a more detailed version will be presented in a forthcoming paper.)

4.2 The method of integration

Our 24 time-dependent variables are as follows:

The values of these variables were calculated by direct integration of the equations of motion by a seventh order Runge-Kutta-Nyström integrator (Fehlberg 1974). Energy and angular momentum were typically conserved to 1 part in 105 and 108 respectively.

The determination of the osculating orbital elements from the cartesian coordinates and velocities imposes no difficulty. Two sets of orbital elements were calculated. Besides the commonly used osculating orbital elements we calculated a second set of "orbital elements'', which give a better description of the real orbit. For the determination of this set of orbital parameters we formally replaced the mass parameter $\mu_1=Gm_{12}$ by the $\mu_1^*=-\vec{\rho}_1\cdot\vec{\ddot\rho}_1\rho_1$ expression. (Cf. Eqs. (21) and (22) in Kiseleva et al. 1998.) The main advantage of this choice is that if the orbit is exactly circular, then these "elements'' give the real orbit of the binary at any moment, which in turn is not true for the osculating elements (e.g. the eccentricity of the osculating orbit will never be zero). But even when the eccentricity is not zero and non-radial forces are present (due to e.g. a third body, inclined rotational axes, or the tidal lag caused by the dissipative forces, which is not treated here), the effect of the evidently dominating radial force component is eliminated from the difference between the real and this "super-osculating'' orbit in the neighbourhood of the calculation of the orbital elements. The above-mentioned difference will be

 
$\displaystyle \Delta\vec{r}(t_0+\Delta t)=\frac{1}{2}(\Delta t)^2[\vec{f}_{\rm ... ...t_0)-\left(\frac{\dot{\mu}_1^*}{\rho_1^3}\vec{\rho}_1\!\right)(t_0)\right]+...,$     (23)

where $\vec{f}_{\rm t}$ and $\vec{f}_{\rm n}$ stand for the transversal and the normal force-components, respectively. If $\mu_1^*$ depends (almost) only upon $\rho_1$, and the contributions of the non-radial force terms are small, then it can be seen easily that if the "super-osculating'' elements were calculated near the periastron (or apastron, i.e. when $\vec{\rho}_1\cdot\vec{\dot\rho}_1\approx0$), the third order terms also diminish in the (23) difference. This was our main reason to calculate the orbital elements at the integration step which is the closest to the periastron of the binary.

In this paper we mainly concentrate on the possible observational consequences of the perturbed motion of the eclipsing binary. For this reason, the angular orbital elements were calculated with respect to the plane of the sky throughout this work.

4.3 Results


  
Table 7: The initial orbital elements (and astrophysical parameters of the binary stars) for the different integration runs. (The angular orbital elements refer to the plane of the sky, with an arbitrary initialdirection.)
\begin{table} \begin{displaymath} \begin{array}{lllll} \hline\hline \noalign{\sm... ...87&0.0087\\ \noalign{\smallskip }\hline \end{array}\end{displaymath}\end{table}

In Figs. 7a,b, 9a,b one can see the variations of the osculating orbital elements in the mass-point approximation both in the low- and the high inclination regime calculated by us,

  \begin{figure} \par\includegraphics[width=5.7cm,clip]{3529F7a.eps}\hspace*{2mm} \includegraphics[width=5.7cm,clip]{3529F7b.eps}\end{figure} Figure 7: The variations of the orbital elements of our system during 10, 10 000revolutions of the outer body ($\approx $37.6, $37.6\times 10^3$ y) according to the point mass model. (High initial mutual inclination.)
Open with DEXTER

and hereafter referred to as the first and second runs, respectively (see the corresponding columns in Table 7). This behaviour is the same (as expected) as that described in Ford et al. (2000) (see their Sect. 3.2, and their Figs. 7 and 8).
  \begin{figure} \par\includegraphics[width=8cm,clip]{3529F8.eps}\end{figure} Figure 8: The effect of the perturbations in the orbital motion of the binary for the times of the eclipsing minima. (The corresponding variations of the orbital elements can be found in Fig. 7a.)
Open with DEXTER

We also calculated the variations of the theoretical O-C diagram of the eclipsing system (see Figs. 8, 10). (The geometrical light-time effect was removed from these curves.) We think these figures do not need any further comment. They are shown here only for illustration of the previous subsection and to serve as comparison to the following results.

For the next runs the members of the binary were no longer treated as mass points. The kj constants were taken from the tables of Claret & Gimènez (1992). These values are approximately one order smaller than those calculated from the polytropic models (see e.g. Finlay-Freundlich 1958, or most recently Eggleton et al. 1998), which shows less significant deviation from the centrally condensed model (which is of course dynamically equivalent with the mass-point approximation). The completely different dynamical behaviour of hierarchical triples with distorted components was studied first by Söderhjelm (1984), and recently by e.g. Eggleton et al. (1998); Eggleton & Kiseleva-Eggleton (2001). In our work we concentrate only on those variations of the orbital elements which have an effect on the times of minima on a relatively short time scale (decades, or a few centuries).


  \begin{figure} \par\includegraphics[width=5.6cm,clip]{3529F9a.eps}\hspace*{2mm} \includegraphics[width=5.6cm,clip]{3529F9b.eps}\end{figure} Figure 9: The variations of the orbital elements of our system during 10 and 10 000revolutions of the outer body according to the point mass model. (Low initial mutual inclination.)
Open with DEXTER

In the third run (Figs. 11a,b, 12) the integration was started with the same initial parameters than in the first one (see Table 7). However, the initial spin angular velocity vectors were completely synchronized to the instantaneous orbital angular velocity vector. Of course, due to the precession of the orbital plane of the binary caused by the inclined tertiary, the direction of the rotational axes of the stars in the close system will no longer be perpendicular to the orbital plane. As a result, a small precession of the spin axes will arise, which results in a further small amplitude oscillation in the orbital inclination as well as in the node. Although the amplitude of this oscillation is some hundredth of degrees in the present configuration, this yields an oscillation in the times of minima of the order of 10-5 days. For the present this variation is beyond the observability limit, but in the case of closer systems (as e.g. Algol, $\lambda$ Tau) we believe that this precession might be detectable.


  \begin{figure} \par\includegraphics[width=8cm,clip]{3529F10.eps}\end{figure} Figure 10: The effect of the perturbations in the orbital motion of the binary for the times of the eclipsing minima. (The corresponding variations of the orbital elements can be found in Fig. 9a.)
Open with DEXTER


  \begin{figure} \par\includegraphics[width=5.8cm,clip]{3529F11a.eps}\hspace*{2mm} \includegraphics[width=5.8cm,clip]{3529F11b.eps}\end{figure} Figure 11: The variations of the orbital elements of our system during 10 and 10 000revolutions of the outer body according to the distorted model. (See the initial elements in column "3rd run'' of Table 7.)
Open with DEXTER

Comparing the results of the extended distorted body- and the mass-point integrations in the low mutual inclination area (Figs. 13a,b, 14 vs. Figs. 9a,b, 10) there are no significant differences in the secular evolution of the orbits. Despite this fact, the short-term behaviour of the O-C curves shows important alterations. Unfortunately, the larger amplitude O-C variations (which might be observable) arise in the physically unrealistic mass-point case, due to the approximately 2.5 times larger amplitude variation in the inner orbital eccentricity as well as a cyclic variation in the binary's semi-major axis (manifesting itself as the parabolic shape of the corresponding O-C curve), which are absent in the distorted case.

5 Conclusions

As a case study, our investigations concerning the complicated period variations of IM Aur give a good example of how to obtain an acceptable O-C curve solution physically consistent with other observational (photometric, spectroscopic, etc.) facts. It is also well known that some previous results (on other systems) are still inconsistent with the dynamical evolutionary deductions. These LITE solutions cannot be acceptable.


  \begin{figure} \par\includegraphics[width=8cm,clip]{3529F12.eps}\end{figure} Figure 12: The effect of the perturbations in the orbital motion of the binary for the times of the eclipsing minima. (The corresponding variations of the orbital elements can be found in Fig. 11.)
Open with DEXTER


  \begin{figure} \par\includegraphics[width=5.3cm,clip]{3529F13a.eps}\hspace*{2mm} \includegraphics[width=5.48cm,clip]{3529F13b.eps}\end{figure} Figure 13: The variations of the orbital elements of our system during 10 and 10 000revolutions of the outer body according to the distorted model. (See the initial elements in column 4th run of Table 7.)
Open with DEXTER

We confirmed the earlier LITE solution for IM Aur, while a further kind of period variation is reported for the first time. The results of analysis of the light curves and of the minima timings are not contradictory, and combining all available information, we can predict a 1.7-2.5 $M_\odot$ mass A8V-A2V type star with an 1372 days long period, eccentric orbit (with observable inclination $i'>55\hbox{$^\circ$ }$).

In spite of the obvious fact that a third (not so distant) physical companion of an eclipsing binary should always cause perturbations on the close eclipsing orbit, this effect is not yet generally studied simultaneously with an O-C analysis. The main reason is probably that the O-C investigators usually are the same as the photoelectric minima time observers, while dynamical examinations of the triple (and other type) stellar systems are usually carried out in separate groups of mainly mathematicians. However, as we may conclude, the one-by-one numerical study of each investigated candidate system, together with the standard O-C analysis procedure, can reveal spectacular variations even on a time scale comparable with the third-body orbit. Although neither the amplitude nor the period of the resulting O-C changes make it possible to show the expected perturbations on the O-C curve (as well as on the light curve) in the specific case of IM Aur, in many other sample systems having much closer or more massive third members, the expected changes should be observable. We would like to emphasize the necessity of taking into account these dynamical effects during any similar calculations (see also at Söderhjelm 1975; Mayer 1990).


  \begin{figure} \par\includegraphics[width=8.4cm,clip]{3529F14.eps}\end{figure} Figure 14: The effect of the perturbations in the orbital motion of the binary for the times of the eclipsing minima. (The corresponding variations of the orbital elements can be found in Fig. 13a.)
Open with DEXTER

Moreover, it would seem that occasional minima-time hunting is not the right technique for the verification of the actual perturbative effects. If the dynamical behaviour is in the focus of a systematic study, frequent accurate timings of a few selected targets with the same instrumentation and using the same reduction methods (making the material as homogenous as possible) are more effective. The variations in some other important orbital elements (as e.g. the apparent inclination of the eclipsing orbit) can also be appreciable as secular changes of certain features of the light curves (e.g. depth and asymmetry of minima, etc.). These effects also necessitate very thorough sampling of the light curves at the best available accuracy.

Using our methods and the developed numerical code, a subsequent analysis of some exciting multiple systems having exotic O-C patterns will be carried out in the near future. It should be very important to take into consideration during the apsidal motion analysis of such eccentric binary systems, which have clear evidence of third physical members, so that e.g. the theoretical k2 parameters could be corrected for the perturbations.

Acknowledgements
This work was supported by the OTKA F030147, and T030743 grants of the National Scientific Research Foundation (Hungary). This research has made use of NASA's Astrophysics Data System Bibliographic Services, as well as the SIMBAD database operated at CDS, Strasbourg, France. We also thank Dr. Dirk Terrell, and Dr. László Szabados for their comments on the manuscript, as well as Dr. Jet Katgert for the final stylistic corrections.

References

 

Copyright ESO 2002