A&A 398, 1091-1102 (2003)
DOI: 10.1051/0004-6361:20021688

On the detectability of long period perturbations in close hierarchical triple stellar systems

T. Borkovits¹ - B. Érdi² - E. Forgács-Dajka² - T. Kovács^2,

1 - Baja Astronomical Observatory of Bács-Kiskun County, PO Box 766, 6500 Baja, Szegedi út, Hungary
2 - Department of Astronomy, Eötvös Loránd University, 1117 Budapest, Pázmány P. sétány 1/A, Hungary

Received 26 September 2002 / Accepted 12 November 2002

Abstract
We study the possibility of the detection of the low amplitude long (P') period perturbative effect of a distant third companion on the motion of a close binary. We give a new, more accurate analytical formula for this kind of perturbation affecting the moments of the times of minima in eclipsing binaries. The accuracy of this formula is tested by numerical integrations carried out for several initial configurations. We also describe a numerical method based on a non-linear Levenberg-Marquardt algorithm which makes it possible to separate this dynamical effect from the pure geometrical light-time effect in the eclipsing O-C diagram. The capabilities of this new method are demonstrated by the analysis of numerically simulated O-Cs for test systems having physical parameters very similar to Algol and IU Aur. The results show that the above mentioned effect would be detectable in these systems nowadays, observing almost each minima events in a 1-2 year-long interval.

Key words: methods: analytical - methods: numerical - celestial mechanics - binaries: close - binaries: eclipsing

1 Introduction

Several close binary stars have third, distant companion. Due to the presence of this further component, the motion of the binary no longer will be purely Keplerian, but different types of periodic and non-periodic (secular) perturbations would occur. According to the classification of Brown (1936) the periodic perturbations can be divided into the following three groups:

• Short period perturbations. The typical period is equal to the orbital period P of the close binary, while the amplitude has the order (P/P')² (where P' denotes the period of the wide orbit);
• Long period perturbations. This group has a typical period of P', and magnitude of the order (P/P');
• Apse-node terms. In this group the typical period is about P'²/P, and the order of the amplitude reaches unity.

(We have to note that this classification differs from what is used in the classical planetary perturbation theories. There the first two of our groups called together as "short period'' perturbations, while the "apse-node terms'' are called as "long period'' ones. In the stellar three-body problem this latter classification was used by Harrington 1968, 1969.) These effects can be most easily detected in those triple systems, where the close binary happens to be an eclipsing one. This follows from different reasons. First, the usual orbital periods of the eclipsing binaries are several days, so in favourable cases even the apse-node terms appear in a time-scale of some decades or centuries, which for nowadays almost can be covered at least for a few systems. Furthermore, the variation of the orbital elements may produce very spectacular effects in the characteristics of the eclipses. Here we refer to the variable eclipse-depth at some eclipsing binaries in a time-scale of decades. The well-known examples are SS Lac (Torres 2001, and further references therein), V907 Sco (Lacy et al. 1999), SV Gem (Guilbault et al. 2001), where the eclipses disappeared in the last decades, furthermore, the yet-eclipsing binary IU Aur also shows fast inclination variations (Drechsel et al. 1994). Another important effect in eccentric binaries is the precession of the line of the apsides caused by the third star. Nevertheless, this phenomenon is not so easily observable in triple systems, since the main sources of the apsidal motion in the known cases are the tidal forces arising from the close proximity of the members of such binaries. In the cases of some systems with abnormally slow apsidal motion the superposition of the tide-generated and the third body-forced apsidal motion might explain the discrepancies between the theory and the observations (see e.g. Khodykin & Vedeneyev 1997; Kozyreva et al. 1999 for the binary AS Cam).

A further main advantage of dealing with eclipsing binaries is that all of the aforementioned phenomena (as well as further ones) affect the occurrence of eclipse events too. Apart from several other physical mechanisms which can modify the observable eclipsing minima times (e.g. mass flow in/from the system, tidal forces etc.), the effect of the third body on the eclipsing O-C diagram can be divided into a geometrical and a dynamical part. The geometrical contribution is the well-known light-time effect. This reflects the motion of the eclipsing pair around the centre of mass of the triple system. If its quasi-sinusoidal pattern can be separated from the other distortions of the O-C curve, some of the orbital elements of the wide orbit (P', e', $\omega'$, $\tau '$, $a_{12}\sin{i'}$, where a₁₂ denotes the semi-major axis of the orbit of the eclipsing pair around the common centre of mass of the triple system) can be determined. (Perhaps the easiest way of this calculation was introduced by Kopal 1978.)

The dynamical contributions arise from the different perturbations. The typical amplitudes of these terms are listed in Söderhjelm (1975). During an apse-node cycle the magnitude of the O-C variations can reach even the order of days. Nevertheless, on a time-scale which is significantly shorter than the apse-node period this variation can be manifested e.g. as a parabolic pattern in the O-C curve, and its nature very easily can be misinterpreted (for an extended discussion see Borkovits et al. 2002).

Up to this moment we mainly concentrated only for the largest amplitude apse-node terms. Nevertheless, the O-C diagram might give a unique possibility to detect some kinds of long periodic perturbations. Of course, the O-C curve reflects the long, and even the short periodic variations of the orbital elements in the same way as in the case of the apse-node time-scale perturbations, but the amplitude of such variations usually much smaller than the limit of observability. The only exception (at least in some cases) arises from the direct (long period) perturbation in the mean motion of the close binary, which due to its cumulative effect on the O-C diagram may exceed the limit of detectability.

In this paper we concentrate on this long period contribution of the O-C diagram. In Sect. 2 we give an analytical formula for this term which is valid for arbitrary value of the mutual inclination, although only for nearly circular close orbits. We also compare our result with the outputs of direct numerical integrations. In Sect. 3 we present a numerical process to separate the dynamical and geometrical term of the O-C, and we examine whether and how this dynamical term can be used to determine the real mutual spatial orientation of the orbits. We illustrate our results with numerically generated O-C data. Finally in Sect. 4 we consider the chance of the separability in real triple systems.

2 An analytical formula of the long period perturbation of an O-C curve

By the use of the theory of Harrington (1968, 1969), based on the von Zeipel averaging method of the canonical equations, Söderhjelm (1975, 1982) derived analytical formulae for the long period perturbations in the standard Delaunay variables. Although these formulae are exact up to second order in the (a/a') ratio, their practical use is limited, at least in their original forms. Mayer (1990) gave a simple, useful form with the assumptions that the elements of the wide orbit are constant, the close orbit is circular, and the relative orientation of the two planes is invariant. Nevertheless, as it will be shown, the expression of Mayer (1990) maybe somewhat inaccurate. More explanation will come later.

In the following we present a corrected new formula, which is valid for the same assumptions. Such a solution could be get easily from the original formulae of Söderhjelm as it was done by Mayer. Despite this fact we follow a different way. Instead of the perturbing potential we depart from the perturbing force, and we calculate directly the perturbations in the eclipsing period in the function of the true anomaly of the eclipsing pair along the wide orbit. As it will be shown this method is more effective and faster for this particular problem, than the usual methods, furthermore, this automatically helps to avoid that kind of inaccuracies which occured in the former solution.

2.1 The expansion of the perturbing force

	Figure 1: The spatial orientation of the orbital planes. See text for details.
Open with DEXTER

Using the mass-point approximation, the perturbing force acting upon the close binary is:

$$% \vec{f}=Gm_3\left(\frac{\vec{r}_{23}}{r_{23}^3}-\frac{\vec{r}_{13}}{r_{13}^3}\right),$$ (1)

where G denotes the gravitational constant, m₃ is the mass of the tertiary, while $\vec{r}_{i3}$ stands for the position vector between the ith component of the binary and the third star. The above expression, as it is well-known, can be expanded into a series of Lagrangian polynomials of the following form:

 

$$% \vec{f}=\frac{Gm_3}{\rho_2^3}\left\{\left[\sum_{n=0}^\infty\lef... ...n\left(\frac{\rho_1}{\rho_2}\right)^nP_n(\lambda)\right]^3\vec{r}_{13}\right\},$$ (2)

where M₁₂ means the total mass of the binary, $\rho_1$ and $\rho_2$ denote the absolute value of the first two Jacobian position-vectors, (e.g. $\rho_1$ means the separation of the members of the eclipsing pair, while $\rho_2$ is the distance between the centre of mass of the binary and the distant third companion), while $\lambda$ stands for the direction cosine between $\vec{\rho}_1$ and $\vec{\rho}_2$. Let us define a Cartesian coordinate system whose origin is at the centre of mass of the binary, and the three axes are parallel with the vectors $\vec{\rho}_1$, $\vec{c}\times\vec{\rho}_1$, $\vec{c}~=\vec{\rho}_1\times\vec{\dot\rho}_1$, respectively. The direction cosines between the vector $\vec{\rho}_2$ and the axes (as it can be seen in Fig. 1) are as follows:

   

$$\lambda=\cos{w}\cos{w'}+\sin{w}\sin{w'}\cos{i}_{\rm m},$$ (3)

$$\mu=-\sin{w}\cos{w'}+\cos{w}\sin{w'}\cos{i}_{\rm m},$$ (4)

$$\nu=\sin{w'}\sin{i}_{\rm m},$$ (5)

where w and w' denote the true longitude of the secondary and the tertiary measured from the intersection of the orbits, and $i_{\rm m}$ is the mutual inclination. According to these the three (radial, transversal and normal) components of the perturbing force in the first order of the ratio $\rho_1/\rho_2$ are

    

$$f_{{\rm r1}}=\frac{3}{8}\frac{Gm_3}{\rho_2^2}\frac{\rho_1}{\rho_2}\left[\left(1+I^2\right)\cos\left(2w'-2w\right) \right.$$

$$\quad \left. +\left(1-I^2\right)\cos\left(2w'+2w\right)+2\left(1-I^2\right)\cos2w\right],$$ (6)

$$f_{{\rm r2}}=\frac{3}{4}\frac{Gm_3}{\rho_2^2}\frac{\rho_1}{\rho_2}\left[\left(1-I^2\right)\cos2w'+\frac{2}{3}P_2(I)\right],$$ (7)

$$f_{\rm t}=\frac{3}{8}\frac{Gm_3}{\rho_2^2}\frac{\rho_1}{\rho_2}\left[(1+I)^2\sin(2w'-2w) \right.$$

$$\quad \left. -(1-I)^2\sin\left(2w'+2w\right)-2\left(1-I^2\right)\sin2w\right],$$ (8)

$$f_{\rm n}=\frac{3}{4}\frac{Gm_3}{\rho_2^2}\frac{\rho_1}{\rho_2}\left[2\cos{w}\sin2w'\sin{i}_{\rm m} \right.$$

$$\quad \left. +\left(1-\cos2w'\right)\sin{w}\sin2i_{\rm m}\right],$$ (9)

where

$$% I=\cos{i}_{\rm m}.$$ (10)

(As it can be seen, the radial force component is divided into two parts. The first one contains terms depending on w, and very similar to the $f_{\rm t}$ transversal component, while $f_{{\rm r2}}$ does not depend on the revolution of the eclipsing binary.)

In what follows we will refer the orbital elements to a plane perpendicular to the line of sight of the observer, and going accross the centre of mass of the binary. We will call it as the plane of the sky. It is clear that the distance of this plane from the observer varies in time, according to the

$$% \Delta z=r\sin{i'}\sin\left(v'+\omega'\right)$$ (11)

function, where

$$% r=\frac{m_3}{M_{123}}\rho_2,$$ (12)

furthermore v' denotes the true anomaly of the outer body, $\omega'$ is the argument of the periastron of the binary's orbit around the centre of mass of the triple system, and M₁₂₃ stands for the total mass of the triple. As it is well-known this motion is the source of the light-time effect detected in several triple systems.

Using the above mentioned true anomaly, v', of the outer body, and the true longitude, u, of the secondary measured from the plane of the sky, the force-components have the following forms:

    

$$f_{{\rm r1}}=\frac{3}{8}\frac{Gm_3}{\rho_2^2}\frac{\rho_1}{\rho_2}\left[\left(1+I^2\right)\cos\left(2v'-2u-\alpha\right) \right.$$

$$\quad \left. +\left(1-I^2\right)\cos\left(2v'+2u-\beta\right)+2\left(1-I^2\right)\cos2\left(u-u_{\rm m}\right)\right],$$ (13)

$$f_{{\rm r2}}=\frac{3}{4}\frac{Gm_3}{\rho_2^2}\frac{\rho_1}{\rho_2}\left[\left(1-I^2\right)\cos2\left(v'-v'_{\rm m}\right)+\frac{2}{3}P_2(I)\right],$$ (14)

$$f_{\rm t}=\frac{3}{8}\frac{Gm_3}{\rho_2^2}\frac{\rho_1}{\rho_2}\left[(1+I)^2\sin\left(2v'-2u-\alpha\right) \right.$$

$$\quad \left. -(1-I)^2\sin\left(2v'+2u-\beta\right)-2\left(1-I^2\right)\sin2\left(u-u_{\rm m}\right)\right],$$ (15)

$$f_{\rm n}=\frac{3}{4}\frac{Gm_3}{\rho_2^2}\frac{\rho_1}{\rho_2}\left\{2\cos(u-u_{\rm m})\sin2\left(v'-v'_{\rm m}\right)\sin{i}_{\rm m} \right.$$

$$\quad \left. +\left[1-\cos2\left(v'-v'_{\rm m}\right)\right]\sin(u-u_{\rm m})\sin2i_{\rm m}\right\},$$ (16)

where the phase angles are

$$% \alpha 2v'_{\rm m}-2u_{\rm m},$$ = (17)

$$\beta 2v'_{\rm m}+2u_{\rm m}.$$ = (18)

(In the above expressions $u_{\rm m}$ refers to the true longitude of the intersection of the two orbits measured along the inner orbit, while $u'_{\rm m}=v'_{\rm m}+\omega'$ has the same meaning for the outer one.)

2.2 Calculation of the O-C

Using the expansions (13)-(16) the analytical form of the O-C can be calculated very easily. To do that we depart from the well-known fact, that at the moment of a mid-minimum

$$% u=\pm\frac{\pi}{2}+2k\pi,$$ (19)

where k is an integer. (Strictly speaking the above equation is valid exactly only if the eccentricity of the eclipsing binary is zero, or the visible inclination is 90 $\hbox{$^\circ$ }$, nevertheless, in our treatment the first condition is practically fulfilled.) Let us define the so called instantaneous period of the binary in the following way:

$$% {\cal{P}}=\frac{2\pi}{\dot{u}},$$ (20)

where as it is well-known (see e.g. Milani et al. 1987, Chap. 3.2):

$$% \dot{u}=\frac{c}{\rho_1^2}-\dot\Omega\cos{i}.$$ (21)

Here we note an important fact. The first term on the rhs. is independent from the plane of reference, while the second one has different values using different reference planes. (Please keep in your mind that in the above expression, as well as in the following ones $\vec{c}$ means the special angular momentum of the binary, e.g. the length of vector $\vec{c}$, and not the velocity of light.)

Furthermore, let us denote by $\overline{P}_i$ the elapsed time between the ith and the (i+1)th eclipsing (let's say primary) minima. (Hereafter we refer to $\overline{P}$ as eclipsing period.) Then it can be seen easily that

$$% \overline{P}_i=\int_i^{i+1}{\cal{P}}{\rm d}\phi,$$ (22)

where

$$% \phi=\frac{u-u_0}{2\pi}\cdot$$ (23)

According to this the eclipsing period is the average of the instantaneous period during a revolution. Consequently the occurrence of the Nth primary minimum after a t₀ epoch can be determined as

$$% \sum_{k=0}^{N-1}\overline{P}_k=\int_0^N{\cal{P}}{\rm d}\phi.$$ (24)

Consequently, if the ${\cal{P}}$ (or $\dot{u}$) versus $\phi$ (or u) dependence is known, the theoretical form of the O-C curve can be calculated formally by an integration. (Of course, the O-C curve is not a continuous function, only at the integer values of the independent variable has physical meaning.)

To get this relation we rewrite $\dot{u}$ as

$$% \dot{u}=n_0\left(1+\frac{\dot{u}_{\rm p}}{n_0}+\frac{\dot{u}_{\rm n}}{n_0}\right),$$ (25)

where n₀ stands for the mean motion of the unperturbed two-body revolution in a fixed (let's say t =t₀) moment, and

$$% \dot{u}_{\rm p}=\frac{c}{\rho_1^2}-n_0,$$ (26)

while

$$% \dot{u}_{\rm n}=-\dot\Omega\cos{i}.$$ (27)

As far as the perturbations are small in the mean motion the relation between the instantaneous period ${\cal{P}}$ and the Keplerian period P₀ of the unperturbed motion can be written as

$$% {\cal{P}}=P_0\left[1-\frac{P_0}{2\pi}(\dot{u}_{\rm p}+\dot{u}_{\rm n})\right]\cdot$$ (28)

(In what follows we omit the "0'' subscripts from the initial values of the quantities referring to the unperturbed motion.)

First we calculate the effect of the $\dot{u}_{\rm p}$ component on the instantaneous period. It can be easily seen, that

$$% \frac{{\rm d}\dot{u}_{\rm p}}{{\rm d}\phi}=\frac{{\rm d}\do... ...rm d}t}\frac{{\rm d}t}{{\rm d}\phi}=\ddot{u}_{\rm p}{\cal{P}}.$$ (29)

The second derivative of the $u_{\rm p}$ part of the true longitude has the following form (cf. Milani et al. 1987, Chap. 3.2):

$$% \ddot{u}_{\rm p}=\frac{f_{\rm t}}{\rho_1}-2\mu\frac{e\sin{v}}{\rho_1^3},$$ (30)

where

$$% \mu=GM_{12},$$ (31)

and v denotes the true anomaly of the secondary component. Here we have to note an important fact. For the first sight the presence of the second term in Eq. (30) contradicts our previous assumption, that the orbit of the binary is circular. In fact it is not true. Although the eccentricity is close to zero, it cannot be permanently exactly zero in perturbed systems. Even if at some moment the close orbit was circular in the next moment due to the perturbing forces it would not be that. (For the possible astrophysical importance of this small non-zero eccentricity especially for semi-detached systems see Eggleton et al. 1998.) So in nearly circular systems the eccentricity oscillates between zero and a small value (typically some ten-hundredthousandths, see e.g. our numerical integrations for the system IM Aur in Borkovits et al. 2002). In that case e approximately has the same magnitude as $\dot{e}$, and so it can be shown that the two terms on the rhs of (30) may have the same order. Furthermore, since both $f_{\rm t}$ and $e\sin{v}$ have the same order of magnitude as the small variation in the $\rho_1$ radius, the denominators in (30) can be replaced by a constant average distance, which is the Keplerian semi-major axis a. On the other hand, we note that this is important only in the case of the $e\sin{v}$ term of Eq. (30), as in the first expression the multiplicator $\rho_1$ in the amplitude of $f_{\rm t}$ cancels the denominator.

Let us define the expessions

  

$$\Delta\dot{u}_1=\frac{1}{a}\int_0^\phi f_{\rm t}{\rm d}\phi',$$ (32)

$$\Delta\dot{u}_2=-2\frac{\mu}{a^3}\int_0^\phi e\sin{v}{\rm d}\phi',$$ (33)

respectively. Then

$$% \dot{u}_{\rm p}\approx(\dot{u}_{\rm p})_0+P(\Delta\dot{u}_1+\Delta\dot{u}_2).$$ (34)

For the evaluation of the first term we have to express the true anomaly of the third component by $\phi$. This can be done in two steps. First, we can change from the true anomaly to the mean anomaly by the use of the expansions of Cayley (1861), and after that we approximate the mean anomaly l' by the following formula:

$$% l'=2\pi\frac{P}{P'}\phi+l'_0,$$ (35)

where l'₀ is the mean anomaly at the epoch t₀. Now using the expressions

  

$$\int_0^N\cos\left(jl'\pm{ku}\right){\rm d}\phi= \pm\frac{1}{2\pi}... ...\sin\left(jl'\pm{ku}\right)\left(1\mp\frac{j}{k}\frac{P}{P'}\right)\right]_0^N,$$ (36)

$$\int_0^N\sin\left(jl'\pm{ku}\right){\rm d}\phi=\mp\frac{1}{2\pi}\... ...[\cos\left(jl'\pm{ku}\right)\left(1\mp\frac{j}{k}\frac{P}{P'}\right)\right]_0^N$$ (37)

the evaluation of (32) is trivial, and we get that

$$% \Delta\dot{u}_1(\phi)=\frac{1}{4\pi}\left[\frac{f_{{\rm r1}}}{a}\right]^\phi_0\left[1+{\cal{O}}\left(P/P'\right)\right].$$ (38)

Our next task is the calculation of (33). The dependence of the integrand on $\phi$ can be written as

$$% e\sin{v}(\phi) \left[(e\cos\omega)_0+\int_0^\phi\frac{{\rm d}}{{\rm d}\phi}(e\cos\omega){\rm d}\phi'\right]\sin{u}(\phi)$$ =

$$-\left[(e\sin\omega)_0+\int_0^\phi\frac{{\rm d}}{{\rm d}\phi}(e\sin\omega){\rm d}\phi'\right]\cos{u}(\phi),$$ (39)

where the further integrands can be evaluated with good approximation as

$$% \frac{{\rm d}}{{\rm d}\phi}(e\cos\omega)\approx\frac{1}{a}\frac{P^2}{2\pi}\left(f_{\rm r}\sin{u}+2f_{\rm t}\cos{u}\right),$$ (40)

and

$$% \frac{{\rm d}}{{\rm d}\phi}(e\sin\omega)\approx\frac{1}{a}\frac{P^2}{2\pi}\left(2f_{\rm t}\sin{u}-f_{\rm r}\cos{u}\right)$$ (41)

(see e.g. Milani et al. 1987, Chap. 3.2). Performing the integrations we obtain that

$$% \Delta\dot{u}_2(\phi)=\Delta\dot{u}_2^*(\phi)+\frac{2}{3\pi... ...{\rm r1}}}{a}(\phi)-\frac{1}{\pi}\frac{f_{{\rm r2}}}{a}(\phi),$$ (42)

where $\Delta\dot{u}_2^*$ contains the constant terms and those depending upon only u (via some trigonometric functions). These terms will not give any contribution to the values of the O-C, since they have the same values in every minima, which fact is the direct consequence of (19).

Finally we give the $\dot{u}_{\rm n}=\dot{u}_{\rm n}(\phi)$ function. As it is well known

$$% \dot{u}_{\rm n}=-\frac{f_{\rm n}\sin{u}}{\rho_1}\frac{\rho_1^2}{c}\cot{i}.$$ (43)

Since $f_{\rm n}\sin{u}$ has itself the same order as the components of  $\dot{u}_{\rm p}$, the other quantities in (43) can be treated as constants. Consequently, our approximation for $\dot{u}_{\rm n}$ is the following:

 

$$% \dot{u}_{\rm n} -\frac{3P}{16\pi}\cot{i}\frac{Gm_3}{\rho_2^3}\left\{2\sin{i_{\rm m}}\sin{u_{\rm m}}\sin2\left(v'-v'_{\rm m}\right) \right.$$ =

$$+\sin2i_{\rm m}\cos{u_{\rm m}}\left[1-\cos2\left(v'-v'_{\rm m}\right)\right]$$

$$+2\sin{i_{\rm m}}\sin\left(2u-u_{\rm m}\right)\sin2\left(v'-v'_{\rm m}\right)$$

$$-\left. \sin2i_{\rm m}\cos\left(2u-u_{\rm m}\right)\left[1-\cos2\left(v'-v'_{\rm m}\right)\right]\right\}.$$ (44)

Let us turn back to the expression (28) of ${\cal{P}}$. It can be seen easily that all of the above calculated perturbative terms have the order of (P/P')², which is in the order of 10^-2-10^-4 even for the closest hierarchical systems. So, our expansion is verified. A further integration of (28) gives the analytical form of the effect of the long period perturbations on the O-C curve. We keep only the terms which depend also on v'. (The constant terms will give a linear contribution to the O-C, and so they will build up into the observed eclipsing period, while terms which contain pure trigonometric functions of u will disappear.) First let's treat the terms which depend on only the true anomaly v' of the tertiary. For these the integration can be carried out directly with respect to v', using the expression (cf. e.g. Roy 1988, p. 292)

$$% {\rm d}{\phi}=\frac{1}{2\pi}\frac{P'}{P}\frac{\rho_2^2}{a'^2\left(1-e'^2\right)^{\frac{1}{2}}}{\rm d}v'.$$ (45)

Consequently, the amplitude of these integrated expressions is multiplied by $P'/P\gg1$. On the other hand, according to (36) and (37) terms which contain both v' and uafter the integration will have the same order of magnitude than before. Consequently, the terms which depend on purely the orbital motion of the tertiary will be dominant. Keeping only these terms we obtain that

 

$$% {\rm O}{-}{\rm C} \approx \frac{3}{8\pi}\frac{m_3}{M_{123}}\frac{P^2}{P'}\left(1-e'^2\right)^{-\frac{3}{2}}\biggl\{\left(1-I^2\right)\biggl\{\sin2\left(v'-v'_{\rm m}\right)$$

$$+e'\left[\sin\left(v'-2v'_{\rm m}\right)+\frac{1}{3}\sin\left(3v'-2v'_{\rm m}\right)\right]\biggl\}$$

$$+\left(2I^2-\frac{2}{3}\right)\left(v'-l'+e'\sin{v'}\right)$$

$$-\frac{1}{2}\cot{i}\sin{i_{\rm m}}\biggl\{\cos{i_{\rm m}}\cos{u_{\rm m}}\big\{\sin2\left(v'-v'_{\rm m}\right)$$

$$+e'\left[\sin\left(v'-2v'_{\rm m}\right)+\frac{1}{3}\sin\left(3v'-2v'_{\rm m}\right)\right]$$

$$-2\left(v'-l'+e'\sin{v'}\right)\big\}+\sin{u_{\rm m}}\biggl\{\cos\left(2v'-2v'_{\rm m}\right)$$

$$+e'\left[\cos\left(v'-2v'_{\rm m}\right)+\frac{1}{3}\cos\left(3v'-2v'_{\rm m}\right)\right]\biggl\}\biggl\}\biggl\}\cdot$$ (46)

(Kepler's third law has been used for the transformation of the amplitude.)

2.3 Comparison with other analytical and numerical calculations

For a comparison of our result with the formula of Mayer (1990) we enclose here his solution:

 

$$% {\rm O}{-}{\rm C}_{{\rm Mayer}} \frac{3}{8\pi}\frac{m_3}{M_{123}}\frac{P^2}{P'}\left(1-e'^2\right)^{-\frac{3}{2}}\biggl\{(2-Z)\biggl\{\sin2\left(v'+\omega'\right)$$ =

$$+e'\left[\sin\left(v'+2\omega'\right)+\frac{1}{3}\sin\left(3v'+2\omega'\right)\right]\biggl\}$$

$$+\left(Z-\frac{2}{3}\right)\left(v'-l'+e'\sin{v'}\right)\biggl\},$$ (47)

where

$$% Z=\cos{i_{\rm m}}+\cos^2i_{\rm m}.$$ (48)

(We used our notations instead of the original ones, furthermore, some obvious misprints were corrected here.) The fundamental difference between (47) and (46) manifests in the phase of the trigonometric terms. The phasing would be identical if $\omega'$ in Mayer's paper would be measured from the intersection of the two orbital planes. Nevertheless, he used the same notation for the argument of the periastron in the light-time contribution, where $\omega'$ evidently has to be measured from the plane of the sky. However, the two meanings of the $\omega'$ would be identical only if the observational and the dynamical system of references were the same, or if the two orbital planes intersected each other in the plane of the sky. As it is well-known the calculation of the perturbational problems is usually carried out in the dynamical frame of reference, where the fundamental plane is the invariable plane of the system. In the case of the hierarchical triple stellar systems the net angular momentum of the system mainly concentrates in the wide orbit (see e.g. Eq. (26) of Söderhjelm 1975), consequently the plane of the wide orbit is very close to the invariable plane, and in the immovable wide orbit approximation (which was used by Mayer 1990) the two planes become identical. The other discrepancies also arise from the same problem. If the plane of reference is the plane of the wide orbit, $u_{\rm m}\equiv0$, consequently the terms multiplied by  $\sin{u_{\rm m}}$ will disappear.

	Figure 2: The long period dinamical contribution of O-Cs calculated by numerical integration, furthermore, with the analytical formulae presented in this paper, and in Mayer (1990). Upper panels: low mutual inclinations. Middle panels: medium mutual inclinations. Lower panels: high mutual inclinations. (For the exact input parameters see Tables 2, 3.)
Open with DEXTER

In order to illustrate the accuracy of our result, and to compare it with the formula of Mayer (1990) we carried out several numerical integrations with different initial conditions. The description of the integrator can be found in Borkovits et al. (2002). The only alteration applied here is, that the sampling of the Jacobian coordinates and velocities is done after the integration step closest to the center of an eclipse, and not to the vicinity of the periastron. Only mass-point approximation was applied. As initial parameters the physical properties and orbital elements of two well-known close triple systems were chosen (see Tables 1-3). As it can well be seen in Fig. 2, in the exact coplanar case (upper left panel), as well as in the case, when the two orbital planes intersect each other on the plane of the sky (upper right panel) both Mayer's and our results give similarly accurate approximations, while in the other cases the differences are significant.

 

 

Table 1: The initial parameters of the close systems. (The masses are given in solar mass, the period in days, and the angular orbital elements in degrees.) The non-arbitrary parameters are taken from Söderhjelm (1980), Lestrade et al. (1993) for Algol, and from Drechsel et al. (1994) for IU Aur.

System	m1	m2	P	e	i	$\Omega$	u
"Algol AB''	3.7	0.8	2.8673	0.0	82.3	52	60
"IU Aur AB''	21.3	14.4	1.811474	0.0	88.0	60	90

 

 

Table 2: The fixed initial parameters of the wide systems. The mass function f(m₃) is calculated from the amplitude of the O-C curve, and is given in solar mass. The period P' is given in days, while the periastron passage $\tau '$ in HJD-2 400 000.

System	f(m3)	P'	e'	$\tau '$
"Algol AB-C''	0.125	679.9	0.23	50 000.0
"IU Aur AB-C''	1.89	294	0.54a	50 000.0

^a In the (last) run I10 e'=0.24 was chosen. (See text for details.)

 

 

Table 3: Initial parameters which varied in different runs, and three calculated initial quantities. ("A'' runs refer to Algol-like, while "I'' runs to IU Aur-like system.) The angular elements are given in degrees, while m₃ in solar mass. The quantities in the last three columns are calculated. Furthermore, $A_{\rm p}=A^*(1-e'^2)^{-3/2}$ denotes the amplitude of the perturbative term (without its angular dependence), while $A_{\rm L}=A(1-e'^2\cos^2\omega')^{1/2}$ is the same for the light-time contribution.

No.	m3	i'	$\Omega'$	$\omega'$	$A_{\rm p}/A_{\rm L}$	$i_{\rm m}$	$u'_{\rm m}$
A1	1.7	82.3	142	60	0.10	89.0	82.4
A2	1.7	82.3	97	60	0.10	44.6	86.8
A3	1.7	82.3	52	60	0.10	0.0	-
A4	1.7	82.3	7	60	0.10	44.6	266.8
A5	1.7	82.3	187	60	0.10	-47.4	273.2
A6	1.7	82.3	232	60	0.10	-15.4	257.8
A7	1.7	82.3	277	60	0.10	-47.4	93.2
A8	1.7	82.3	322	60	0.10	89.0	97.6
A9	1.7	82.3	322	150	0.10	89.0	97.6
A10	2.0	60.0	142	60	0.11	86.2	83.3
A11	2.0	60.0	97	60	0.11	47.6	108.5
A12	2.0	60.0	52	60	0.11	22.3	180.0
A13	4.2	30.0	142	60	0.18	83.3	86.1
A14	4.2	30.0	97	60	0.18	62.2	127.6
A15	4.2	30.0	52	60	0.18	52.3	180.0
I1	17.5	88.0	150	5	0.14	89.0	88.0
I2	17.5	88.0	105	5	0.14	45.0	89.1
I3	17.5	88.0	60	5	0.14	0.0	-
I4	17.5	88.0	150	60	0.12	89.0	88.0
I5	17.5	88.0	105	60	0.12	45.0	89.1
I6	17.5	88.0	60	60	0.12	0.0	-
I7	27.8	45.0	150	60	0.16	88.6	88.6
I8	27.8	45.0	105	60	0.16	58.4	123.9
I9	27.8	45.0	60	60	0.16	43.0	180.0
I10	17.5	88.0	105	60	0.08	45.0	89.1

3 Separation of the dynamical term from the O-C

In this section first we show how the presence of the dynamical term can influence the usual method of light-time solutions, and then, we give a numerical method to separate the two terms, which can improve the accuracy of the light-time solution, and, furthermore, may give additional information about the spatial orientation of the triple system.

3.1 The effect of the dynamical term on the light-time solution

A usual way of calculation of the light-time solution is based on the fact that there are some very simple relations (at least in the first and second order in e') between the orbital elements of the wide orbit and the first two or three pairs of coefficients of the Fourier-expansion of the light-time curve, where the fundamental frequency is the period ratio, e.g. $2\pi P/P'$. Consequently, if the harmonic coefficients of the O-C were determined by some numerical methods (typically by weighted least-squares fit), then the orbital elements could be calculated in a very simple way.

For the sake of completeness we describe here the most important formulae after Kopal (1978, Chap. V). In the case of the pure light-time effect the mathematical form of the O-C curve is:

$$% {\rm O}{-}{\rm C}=\sum_{k=1}^\infty\left[a_k\sin(k\nu N)+b_k\cos(k\nu N)\right]-\sum_{k=1}^\infty b_k,$$ (49)

where N is the cycle number, and

$$% \nu=2\pi\frac{P}{P'},$$ (50)

while

$$a_k=A\left[g_k\left(e'\right)\cos\omega'\cos{kl'_0}-h_k\left(e'\right)\sin\omega'\sin{kl'_0}\right],$$ (51)

$$b_k=A\left[g_k\left(e'\right)\cos\omega'\sin{kl'_0}+h_k\left(e'\right)\sin\omega'\cos{kl'_0}\right],$$ (52)

where

$$% A=\frac{a_{12}\sin{i'}}{c},$$ (53)

furthermore,

$$g_k\left(e'\right)=2\sqrt{1-e'^2}\frac{J_k\left(ke'\right)}{ke'},$$ (54)

$$h_k\left(e'\right)=\frac{2}{k}\frac{{\rm d}J_k\left(ke'\right)}{{\rm d}ke'},$$ (55)

and in the latter expressions J_k represents the Besselian function of the kth order. (We note, that in (53) c stands for the velocity of light.) Considering a quadratic approximation in the outer eccentricity the non-zero coefficients are as follows:

  

$$a_1=A\left[\left(1-\frac{3e'^2}{8}\right)\cos\left(l'_0+\omega'\right)-\frac{e'^2}{4}\cos\omega'\cos{l'_0}\right],$$ (56)

$$b_1=A\left[\left(1-\frac{3e'^2}{8}\right)\sin\left(l'_0+\omega'\right)-\frac{e'^2}{4}\cos\omega'\sin{l'_0}\right],$$ (57)

$$a_2=A\frac{e'}{2}\cos\left(2l'_0+\omega'\right),$$ (58)

$$b_2=A\frac{e'}{2}\sin\left(2l'_0+\omega'\right),$$ (59)

$$a_3=A\frac{3e'^2}{8}\cos\left(3l'_0+\omega'\right),$$ (60)

$$b_3=A\frac{3e'^2}{8}\sin\left(3l'_0+\omega'\right).$$ (61)

Using the expansions of Cayley (1861) Eq. (46) also can be easily expanded into trigonometric series of the mean anomaly l', as

 

$$% {\rm O}{-}{\rm C} \approx \frac{3}{8\pi}\frac{m_3}{M_{123}}\frac{P^2}{P'}\biggl\{{\cal{C}}\cos2\left(l'+\omega'\right)+{\cal{S}}\sin2\left(l'+\omega'\right)$$

$$+e'\biggl[{\cal{M}}\sin{l'}-{\cal{C}}\cos\left(l'+2\omega'\right)-{\cal{S}}\sin\left(l'+2\omega'\right)$$

$$+\frac{7}{3}{\cal{C}}\cos\left(3l'+2\omega'\right)+\frac{7}{3}{\cal{S}}\sin\left(3l'+2\omega'\right)\biggl]\biggl\}+{\cal{O}}\left(e'^2\right),$$ (62)

where

$${\cal{C}}=-\left(1-I^2\right)\sin2u'_{\rm m},$$ (63)

$${\cal{S}}=\left(1-I^2\right)\cos2u'_{\rm m},$$ (64)

$${\cal{M}}=6I^2-2.$$ (65)

We omitted terms which are multiplied by $0.5\cot{i}$, as in eclipsing systems the observable inclination of the binary is usually close to $90\hbox{$^\circ$ }$, consequently these terms give only a minor contribution. (E.g. even for $i=80\hbox{$^\circ$ }$, $0.5\cot{i}<0.09$.) A comparison between (62) and (49) reveals that in this case the coefficients are as follows:

  

$$a^*_1=A^*e'\left[{\cal{M}}\cos{l'_0}-{\cal{S}}\cos\left(l'_0+2\omega'\right)+{\cal{C}}\sin\left(l'_0+2\omega'\right)\right],$$ (66)

$$b^*_1=A^*e'\left[{\cal{M}}\sin{l'_0}-{\cal{S}}\sin\left(l'_0+2\omega'\right)-{\cal{C}}\cos\left(l'_0+2\omega'\right)\right],$$ (67)

$$a_2^*=A^*\left[{\cal{S}}\cos2\left(l'_0+\omega'\right)-{\cal{C}}\sin2\left(l'_0+\omega'\right)\right],$$ (68)

$$b_2^*=A^*\left[{\cal{S}}\sin2\left(l'_0+\omega'\right)+{\cal{C}}\cos2\left(l'_0+\omega'\right)\right],$$ (69)

$$a_3^*=A^*\frac{7e'}{3}\left[{\cal{S}}\cos\left(3l'_0+2\omega'\right)-{\cal{C}}\sin\left(3l'_0+2\omega'\right)\right],$$ (70)

$$b_3^*=A^*\frac{7e'}{3}\left[{\cal{S}}\sin\left(3l'_0+2\omega'\right)+{\cal{C}}\cos\left(3l'_0+2\omega'\right)\right],$$ (71)

where the amplitude A^* is

$$% A^*=\frac{3}{8\pi}\frac{m_3}{M_{123}}\frac{P^2}{P'}\cdot$$ (72)

(It will be seen in the next section, that in the case of systems interesting for us, for moderate outer eccentricities, the order of A^* is about Ae'. This is the reason why the quadratic terms were held in (56)-(61).) It is evident that a numerical modelling of the O-C curve in the form (49) yields to coefficients which are the sums of the corresponding (56)-(61) and (66)-(71) terms.

Now some qualitative remarks can be easily done about the effect of the dynamical terms on a usual light-time solution. First, if the third star revolves in a circular orbit, which is coplanar with the orbit of the binary, the dynamical terms diminish, e.g. the geometrical terms are unaffected in that case. Furthermore, for a non-coplanar, but circular outer orbit, the amplitude of the light-time solution remains invariant, at least as far as the quadratic term is not counted. On the other hand, in the case of this configuration the usual determination of the outer eccentricity (from the second Fourier-coefficients) may give an error of several 10 percents. Finally, if the outer orbit is significantly eccentric, both the mass-function (via the amplitude), and the outer eccentricity is affected.

 

Table 4: The results of the parameter search for different runs in the "Algol''-like system. The "L''-rows list the results of the simple light-time solutions. Quantities marked with ^* are calculated for the input values of i', which are also listed in the "L''-rows in parenthesis. (For the other input parameters see Tables 2, 3.)

 

Run	No.	e'	$\omega'$	$\tau '$	a12	i'	D	$i_{\rm m}$	m3	$\chi^2$
			$\hbox{$^\circ$ }$	MHJD	$10^6~{\rm km}$	$\hbox{$^\circ$ }$	$\hbox{$^\circ$ }$	$\hbox{$^\circ$ }$	$M_\odot$	10-7
A1	L	0.22	101	50 071	128*	(82)	(90)		2.0*	848
	1	0.23	57	49 993	140	65	63	62	2.2	678
	2	0.23	57	49 993	130	99	247	113	2.0	719
	3	0.21	148	50 157	226	35	14	49	4.2	837
A2	L	0.22	77	50 034	131*	(82)	(45)		2.0*	359
	1	0.22	59	49 999	128	91	219	140	2.0	216
	2	0.22	59	49 999	128	89	320	40	2.0	217
A3	L	0.25	58	50 005	136*	(82)	(0)		2.1*	266
	1	0.16	54	49 995	175	131	158	141	3.0	400
	2	0.21	73	50 027	141	69	185	151	2.2	637
A4	L	0.22	78	50 036	131*	(82)	(315)		2.0*	397
	1	0.22	59	49 998	135	72	41	41	2.1	205
	2	0.22	58	49 998	129	92	222	138	2.0	216
	3	0.22	58	49 998	129	88	318	42	2.0	216
A5	L	0.25	80	50 039	131*	(82)	(135)		2.0*	514
	1	0.23	56	49 995	143	64	44	46	2.3	365
	2	0.23	56	49 994	130	84	227	132	2.0	411
A6	L	0.24	59	50 006	135*	(82)	(0)		2.1*	237
	1	0.21	64	50 011	134	76	194	155	2.1	493
	2	0.21	66	50 014	134	103	347	25	2.1	554
	3	0.21	67	50 015	148	61	12	24	2.4	564
A7	L	0.17	76	50 031	130*	(82)	(225)		2.0*	448
	1	0.23	59	49 999	128	88	42	43	2.0	258
	2	0.23	60	50 000	135	109	221	139	2.1	266
A8	L	0.17	103	50 077	128*	(82)	(270)		2.0*	907
	1	0.23	58	49 995	144	63	298	62	2.3	607
	2	0.23	57	49 995	129	84	68	68	2.0	645
	3	0.17	168	50 197	235	147	162	129	4.5	722
	4	0.17	167	50 197	187	42	199	122	3.2	729

Run	No.	e'	$\omega'$	$\tau '$	a12	i'	D	$i_{\rm m}$	m3	$\chi^2$
			$\hbox{$^\circ$ }$	MHJD	$10^6~{\rm km}$	$\hbox{$^\circ$ }$	$\hbox{$^\circ$ }$	$\hbox{$^\circ$ }$	$M_\odot$	10-7
A9	L	0.36	166	50 030	139*	(82)	(270)		2.2*	679
	1	0.23	144	49 990	131	87	68	67	2.0	796
	2	0.23	144	49 989	151	61	299	61	2.4	799
	3	0.30	186	50 065	153	115	324	49	2.4	1084
	4	0.29	186	50 066	185	131	146	135	3.1	1169
A10	L	0.22	106	50 081	145*	(60)	(90)		2.3*	968
	1	0.23	57	49 994	142	64	68	67	2.2	778
	2	0.23	57	49 994	133	107	253	108	2.1	811
	3	0.23	155	50 170	231	34	15	50	4.4	958
A11	L	0.16	72	50 025	149*	(60)	(45)		2.4*	408
	1	0.05	147	50 165	212	37	20	48	3.9	238
	2	0.22	62	50 004	128	95	39	41	2.0	276
	3	0.22	63	50 004	139	66	323	39	2.2	302
A12	L	0.26	53	49 993	154*	(60)	(0)		2.5*	280
	1	0.17	40	49 968	177	47	25	41	3.0	265
	2	0.17	39	49 967	149	119	333	45	2.4	289
	3	0.23	71	50 023	144	65	180	147	2.3	635
	4	0.23	73	50 025	166	52	359	31	2.7	744
A13	L	0.24	133	50 129	252*	(30)	(90)		5.0*	1820
	1	0.23	58	49 996	199	40	84	81	3.5	1256
	2	0.34	174	50 204	206	37	198	118	3.7	2332
	3	0.25	59	49 997	543	167	222	108	20.1	9377
A14	L	0.07	355	49 876	252*	(30)	(45)		5.0*	1474
	1	0.22	60	50 002	186	43	133	111	3.2	601
	2	0.23	61	50 003	229	213	42	121	4.3	685
	3	0.21	270	49 713	201	322	33	114	3.6	1325
	4	0.21	268	49 712	250	150	151	123	4.9	1330
A15	L	0.35	31	49 950	264*	(30)	(0)		5.3*	953
	1	0.23	64	50 011	202	321	356	122	3.6	724
	2	0.32	6	49 903	151	60	221	126	2.4	1307

3.2 The numerical method of the separation

The separation of the dynamical terms from the light-time curve may give two important advantages. The first is the determination of more accurate orbital elements, mainly the "projected'' third-body mass. The second one is the possibility to determine the relative spatial orientation of the two orbital planes. This latter arises from the fact, that the dynamical terms - through the orbital elements which determine the light-time orbit - depend on the third body mass (m₃), the observable inclination of the wide orbit (i'), and the mutual inclination ($i_{\rm m}$). It is evident that the dependence on m₃ appears in the amplitude A^*, while the effect of the inclinations manifests through the following equations of spherical triangles:

$$\sin{i_{\rm m}}\cos{u'_{\rm m}}=-\cos{i}\sin{i'}+\sin{i}\cos{i'}\cos\left(\Omega'-\Omega\right),$$ (73)

$$\sin{i_{\rm m}}\sin{u'_{\rm m}}=\sin{i}\sin\left(\Omega'-\Omega\right),$$ (74)

$$\cos{i_{\rm m}}=\cos{i}\cos{i'}+\sin{i}\sin{i'}\cos\left(\Omega'-\Omega\right).$$ (75)

In order to make this separation we developed a computer code which is based on a non-linear Levenberg-Marquardt algorithm (see Press et al. 1989, Chap. 14.4). In the present state the code adjusts the following eight parameters: c₀ (a zero-point correction), P, A (the amplitude of the light-time terms), e', $\omega'$, l'₀, i', and instead of the mutual inclination, $D\equiv\Omega'-\Omega$. (Here we note, that although P is an adjusted parameter, de fundamental frequency $\nu$ is constant during the iteration.) The partial derivatives of the coefficients (56)-(71) are listed in the Appendix.

The method of the computation is the following. As a first step we determine the light-time solution from the O-C curve in the usual way. These orbital elements are used as input parameters for the Levenberg-Marquardt method. For the remaining two parameters (i', D) several initial trial-values are applied automatically. (Using the mass-function, the mass m₃ and the amplitude A^* are calculated in each iteration step.) If a solution is convergent, the program saves the final parameters, and takes the following set of the angular (i', D) parameters.

 

Table 5: The results of the parameter search for different runs in the "IU Aur''-like system. The "L''-rows list the results of the simple light-time solutions. Quantities marked with ^* are calculated for the input values of i', which are also listed in the "L''-rows in parenthesis. (For the other input parameters see Tables 2, 3.)

 

Run	No.	e'	$\omega'$	$\tau '$	a12	i'	D	$i_{\rm m}$	m3	$\chi^2$
			$\hbox{$^\circ$ }$	MHJD	$10^6~{\rm km}$	$\hbox{$^\circ$ }$	$\hbox{$^\circ$ }$	$\hbox{$^\circ$ }$	$M_\odot$	10-7
I1	L	0.40	6	50 000	150*	(88)	(90)		16.0*	2806
	1	0.09	18	50 010	411	160	359	72	71.3	2289
	2	0.47	4	49 999	143	89	318	42	15.1	2694
	3	0.47	4	49 999	144	85	42	42	15.2	2705
I1S	L	0.40	6	50 000	150*	(88)	(90)		16.0*	2691
	1	0.09	15	50 007	402	159	358	71	69.3	2187
	2	0.47	5	49 999	143	88	319	41	15.1	2641
	3	0.47	5	49 999	143	86	40	41	15.2	2665
I2	L	0.49	7	50 003	184*	(88)	(45)		21.0*	2418
	1	0.53	7	50 001	159	88	205	154	17.1	3560
	2	0.53	7	50 001	160	83	25	26	17.2	3639
I3	L	0.55	6	50 002	220*	(88)	(0)		26.7*	3506
	1	0.63	10	50 004	172	94	201	159	19.6	6624
	2	0.63	11	50 004	172	90	159	159	19.4	6840
I4	L	0.46	73	50 004	169*	(88)	(90)		18.8*	3157
	1	0.46	71	49 999	165	86	322	38	18.1	2619
	2	0.46	70	49 999	165	89	38	38	18.0	2623
I5	L	0.50	65	50 006	183*	(88)	(45)		20.8*	2823
	1	0.43	76	50 009	183	108	339	29	20.5	3286
	2	0.43	77	50 010	189	67	20	29	21.6	3469

Run	No.	e'	$\omega'$	$\tau '$	a12	i'	D	$i_{\rm m}$	m3	$\chi^2$
			$\hbox{$^\circ$ }$	MHJD	$10^6~{\rm km}$	$\hbox{$^\circ$ }$	$\hbox{$^\circ$ }$	$\hbox{$^\circ$ }$	$M_\odot$	10-7
I6	L	0.54	55	50 005	203*	(88)	(0)		23.9*	3043
	1	0.39	54	50 002	260	132	155	134	34.2	3501
	2	0.39	54	50 002	247	129	35	47	31.7	3536
	3	0.48	82	40 020	235	55	180	143	29.3	5968
	4	0.49	84	50 021	247	51	359	37	31.5	6275
I7	L	0.43	78	50 006	233*	(45)	(90)		28.9*	3544
	1	0.47	67	49 998	167	77	314	46	18.4	3224
	2	0.47	67	49 998	164	83	134	133	18.0	3239
I8	L	0.43	62	50 000	249*	(45)	(45)		31.7*	2540
	1	0.47	80	50 009	187	66	160	148	21.3	3657
	2	0.47	81	50 009	194	62	341	32	22.4	3809
I9	L	0.52	53	50 000	262*	(45)	(0)		34.2*	3092
	1	0.33	42	49 990	303	216	332	123	42.1	3579
	2	0.49	81	50 015	224	128	6	41	27.1	5511
	3	0.49	82	50 016	236	48	354	40	29.2	5815
I10	L	0.23	75	50 009	184*	(88)	(45)		20.9*	1243
	1	0.22	56	49 997	190	73	42	44	22.1	1341
	2	0.22	56	49 997	188	102	337	44	21.6	1343
	3	0.15	106	50 038	271	138	345	51	36.6	1409

The results of our numerical fitting can be found in Tables 4, 5. In the case of the "Algol-like'' system we deduced the following conclusions:

• In most cases more accurate orbital elements were gained for the elements e', $\omega'$, $\tau '$ than it was obtained from the simple light-time (L) fit. In these runs the $\chi^2$ value also improved with respect to the corresponding L-fits;
• Significant exceptions arise in fits A3, A6, A12, e.g. at those runs, where $D_{{\rm initial}}=0\hbox{$^\circ$ },180\hbox{$^\circ$ }$. In these cases ${\cal{C}}\equiv0$, while ${\cal{S}}=\sin^2i_{\rm m}$, which is also zero for run A3, while ${\cal{S}}\approx0.07$ and ${\cal{S}}\approx0.14$ for runs A6 and A12, respectively. Consequently, in these cases only the ${\cal{M}}$-term, and the coefficients a^*₁, b^*₁ have significant contributions. On the other hand, in run A15, where $D_{{\rm initial}}$ was zero a good fit was also found. In this case already the ${\cal{S}}\approx0.62$ term is the dominant. (Here we note, that in this integration the mutual inclination was $52\hbox{$.\!\!^\circ$ }3$, which is very close to that "critical'' value where the ${\cal{M}}$-term disappears.) It is interesting, that despite the weaker fits in these cases, the $D\approx0\hbox{$^\circ$ }$, $180\hbox{$^\circ$ }$ values were reproduced well;
• The determination of the visible inclination of the outer orbit is less accurate. This is of course not so surprising. We note that the dependence of the tertiary mass on the visible inclination (i') is very weak in the high inclination region. A comparison of the upper or the middle panels in Fig. 2 shows this clearly. Nevertheless, for the last three runs (in the low visible inclination region), where a smaller variation in the visible inclination i' results already significant variation in m₃, really a smaller inclination, and consequently larger third body mass was found;
• Finally, we can conclude that the better fits were reached when the mutual inclination had a medium value.

Considering the "I'' runs with the parameters of the IU Aur system our results are less satisfactory (see Table 5). Although the light-time parameters were improved in more than the half of these runs, the accuracy is far from that achieved in the "A'' runs. What could be the reason? There are three significant differences between the configurations of the two systems. The first is the large outer eccentricity in the IU Aur system, while the other two are the larger masses of the stars, and the smaller P'/P ratio, e.g. the closeness of this triple. The effects of these properties for the behaviour of our solution are very complex. A purely mathematical effect is that due to the larger eccentricity our expression gives a weaker approximation. On the other hand, we note, as perhaps the most important physical effect, that because the latter system is more compact, the apse-node time scale of the largest amplitude perturbations is significantly shorter, and these disturbances may manifest on the O-C diagram within years. (E.g. according to Drechsel et al. 1994 the period of the nodal regression is about 300 years in the IU Aur system.)

In order to study these two above mentioned phenomena we carried out the following two tests. In the last run (I10) we changed the outer eccentricity for the value e'=0.24, while the other parameters had the same values as in run I5. A significant improvement was found in most of the parameters (although a false result also occurred). Furthermore, we calculated the O-C solution of run I1 for a shorter (about half of the original) time interval. (These results are listed in the "I1S'' rows of Table 5.) We did not find significant discrepancy with respect to the original "I1'' solutions. We can conclude, that the larger eccentricity plays a more important role in our weaker results in this case.

4 Conclusions

In this paper we examined the possibility of detecting certain kind of perturbations, which are manifesting themselves on a very short time scale (even during a yearly observing window of the eclipsing variable), but usually omitted due to their low amplitude. This question naturally has different aspects. These are as follows:

• There are only a few known systems where the amplitude of the dynamical term in the O-C exceeds significantly the present observational accuracy. We can easily estimate the maximum-value of the P'/P ratio which is necessary to fulfill this condition. Supposing that the three components have equal masses, and the third body revolves in a circular orbit, according to Eq. (72) the amplitude may exceed the 10^-3 day limit, if P'/P<40P, at least when the two orbits are perpendicular to each other. If the outer orbit has significant eccentricity (say e'=0.50), or the system consists of stars with different masses, then this limit in the function of the mutual inclination may even grow up to $P'/P\approx150P$. From that point of view the two systems for which the numerical integrations were carried out in this paper are placed near the limit of the detectability, as for Algol $P'/P\approx83P$, while for IU Aur the same value is $P'/P\approx91P$.

  The three presently known closest eclipsing triple systems are $\lambda$ Tau ( $P'/P\approx2.11P$), DM Per ( $P'/P\approx13.17P$), and VV Ori ( $P'/P\approx54.41P$) (see e.g. the catalogue of Chambliss 1992). In the cases of the latter two binaries unfortunately only a very few times of minima observations can be found in the literature, although we could expect the largest effects at these stars. On the other hand, we have to note that as the amplitude of the light time-effect decreases with P'^2/3, at these stars already the detection of the pure geometrical effect is also a challenge. In the most interesting case of $\lambda$ Tau, the amplitude of the perturbative terms may be larger with an order of magnitude than that of the light-time terms. Furthermore, at this system due to the large amount of proximity our initial assumptions may loose their validity.

• The second aspect is mainly a technical question. It refers to the observing strategy. This was written already in the conclusion of Borkovits et al. (2002), but, for the sake of the completeness, we repeat it. In order to have any chance for the detection of this phenomenon frequent and accurate timings are necessary. It is desirable to cover a few revolutions of the distant object as densely as possible. The shorter the time interval of such coverage, the smaller the apse-node time scale or secular variations in the orbital elements which could modify the results.
• Finally, the mathematical modelling of perturbations are necessary in order to extract all information from the observations.

In this paper we concentrated mainly on this third item. We calculated a new analytical formula which gives the long period perturbations of the times of minima in eclipsing binaries. We found that this formula is very similar to the earlier expression of Mayer (1990), nevertheless, some errors are corrected. Using this expression we developed a numerical method to separate this dynamical effect from the pure geometrical light time effect. We tested the capabilities of our model by the analysis of numerically simulated O-C curves. In the case of the test runs for outer orbits with moderate eccentricity, significantly better solutions were found than in the larger eccentricity cases.

Naturally the next step would be the application of the method for observed O-C diagrams of real systems. Unfortunately, up to now there are not any O-C diagrams with sufficient accuracy for the possible target systems. This is why we plan to observe some of the few such systems in the near future to collect as many new times of minima as possible.

Appendix A: Partial derivatives for the Levenberg-Marquardt algorithm

As it was mentioned in Sect. 3, the following eight parameters can be adjusted in our non-linear Levenberg-Marquardt code: c₀, P, $A\equiv a_{12}\sin{i'}/c$, e', $\omega'$, l'₀, i', $D\equiv\Omega'-\Omega$. Here we list the analytical form of the partial derivatives of the Fourier-coefficients with respect to these parameters. The a_k, b_k symbols refer to the coefficients of the geometrical light-time effect, while a^*_k, b^*_k denote the dynamical terms. (We note that in the following: $A_k\equiv a_k+a^*_k$, $B_k\equiv b_k+b^*_k$.)

$$\frac{\partial{A_k}}{\partial A}=\frac{a_k}{A}+3\frac{M_{12}}{M_{123}+M_{12}}\frac{a^*_k}{A},$$ (A.1)

$$\frac{\partial{B_k}}{\partial A}=\frac{b_k}{A}+3\frac{M_{12}}{M_{123}+M_{12}}\frac{b^*_k}{A},$$ (A.2)

$$\frac{\partial{A_1}}{\partial{e'}}=-A\frac{e'}{2}\left[\frac{3}{2}\cos\left(l'_0+\omega'\right)+\cos\omega'\cos{l'_0}\right]+\frac{a^*_1}{e'},$$ (A.3)

$$\frac{\partial{B_1}}{\partial{e'}}=-A\frac{e'}{2}\left[\frac{3}{2}\sin\left(l'_0+\omega'\right)+\cos\omega'\sin{l'_0}\right]+\frac{b^*_1}{e'},$$ (A.4)

$$\frac{\partial{A_2}}{\partial{e'}}=\frac{a_2}{e'},$$ (A.5)

$$\frac{\partial{B_2}}{\partial{e'}}=\frac{b_2}{e'},$$ (A.6)

$$\frac{\partial{A_3}}{\partial{e'}}=\frac{2a_3+a^*_3}{e'},$$ (A.7)

$$\frac{\partial{B_3}}{\partial{e'}}=\frac{2b_3+b^*_3}{e'},$$ (A.8)

$$\frac{\partial a_1}{\partial\omega'}=A\left[\left(A.1-\frac{3e'^2}{8}\right)\sin\left(l'_0+\omega'\right)-\frac{e'^2}{4}\sin\omega'\cos{l'_0}\right],$$ (A.9)

$$\frac{\partial b_1}{\partial\omega'}=-A\left[\left(A.1-\frac{3e'^2}... ...right)\cos\left(l'_0+\omega'\right)+\frac{e'^2}{4}\sin\omega'\sin{l'_0}\right],$$ (A.10)

$$\frac{\partial a^*_1}{\partial\omega'}=-2A^*e'\left[{\cal{C}}\cos\left(l'_0+2\omega'\right)+{\cal{S}}\sin\left(l'_0+2\omega'\right)\right],$$ (A.11)

$$\frac{\partial b^*_1}{\partial\omega'}=2A^*e'\left[{\cal{S}}\cos\left(l'_0+2\omega'\right)-{\cal{C}}\sin\left(l'_0+2\omega'\right)\right],$$ (A.12)

$$\frac{\partial A_2}{\partial\omega'}=-\left(b_2+2b^*_2\right),$$ (A.13)

$$\frac{\partial B_2}{\partial\omega'}=a_2+2a^*_2,$$ (A.14)

$$\frac{\partial A_3}{\partial\omega'}=-\left(b_3+2b^*_3\right),$$ (A.15)

$$\frac{\partial B_3}{\partial\omega'}=a_3+2a^*_3,$$ (A.16)

$$\frac{\partial A_k}{\partial l'_0}=-k\left(b_k+b^*_k\right),$$ (A.17)

$$\frac{\partial B_k}{\partial l'_0}=k\left(a_k+a^*_k\right),$$ (A.18)

$$\frac{\partial A_1}{\partial i'}=-3\frac{M_{12}}{M_{123}+2M_{12}}a^*_1\cot{i'}$$

$$\qquad\quad +A^*e'\sin2i_{\rm m}\big\{\cos{u'_{\rm m}}\left[6\cos{l'_0}+\cos\left(l'_0+2\omega'\right)\right]$$

$$\qquad\quad +\sin{u'_{\rm m}}\sin\left(l'_0+2\omega'\right)\big\},$$ (A.19)

$$\frac{\partial B_1}{\partial i'}=-3\frac{M_{12}}{M_{123}+2M_{12}}b^*_1\cot{i'}$$

$$\qquad\quad +A^*e'\sin2i_{\rm m}\big\{\cos{u'_{\rm m}}\left[6\sin{l'_0}+\sin\left(l'_0+2\omega'\right)\right]$$

$$\qquad\quad -\sin{u'_{\rm m}}\cos\left(l'_0+2\omega'\right)\big\},$$ (A.20)

$$\frac{\partial A_2}{\partial i'}=-3\frac{M_{12}}{M_{123}+2M_{12}}a^*_2\cot{i'}$$

$$\qquad\quad -A^*\sin2i_{\rm m}\big[\cos{u'_{\rm m}}\cos2\left(l'_0+\omega'\right)$$

$$\qquad\quad +\sin{u'_{\rm m}}\sin2\left(l'_0+\omega'\right)\big],$$ (A.21)

$$\frac{\partial B_2}{\partial i'}=-3\frac{M_{12}}{M_{123}+2M_{12}}b^*_2\cot{i'}$$

$$\qquad\quad -A^*\sin2i_{\rm m}\big[\cos{u'_{\rm m}}\sin2\left(l'_0+\omega'\right)$$

$$\qquad\quad -\sin{u'_{\rm m}}\cos2\left(l'_0+\omega'\right)\big],$$ (A.22)

$$\frac{\partial A_3}{\partial i'}=-3\frac{M_{12}}{M_{123}+2M_{12}}a^*_3\cot{i'}$$

$$\qquad\quad -A^*\frac{7e'}{3}\sin2i_{\rm m}\big[\cos{u'_{\rm m}}\cos\left(A.3l'_0+2\omega'\right)$$

$$\qquad\quad +\sin{u'_{\rm m}}\sin\left(3l'_0+2\omega'\right)\big],$$ (A.23)

$$\frac{\partial B_3}{\partial i'}=-3\frac{M_{12}}{M_{123}+2M_{12}}b^*_3\cot{i'}$$

$$\qquad\quad -A^*\frac{7e'}{3}\sin2i_{\rm m}\big[\cos{u'_{\rm m}}\sin\left(3l'_0+2\omega'\right)$$

$$\qquad\quad -\sin{u'_{\rm m}}\cos\left(3l'_0+2\omega'\right)\big],$$ (A.24)

$$\frac{\partial A_1}{\partial D}=2A^*e'\big\{I_1\left[6\cos{l'_0}-\cos\left(l'_0+2\omega'\right)\right]-I_2\sin\left(l'_0+2\omega'\right)$$

$$\qquad\quad -\cos{i'}\left[{\cal{C}}\cos\left(l'_0+2\omega'\right)+{\cal{S}}\sin\left(l'_0+2\omega'\right)\right]\big\},$$ (A.25)

$$\frac{\partial B_1}{\partial D}=2A^*e'\big\{I_1\left[6\sin{l'_0}-\sin\left(l'_0+2\omega'\right)\right]+I_2\cos\left(l'_0+2\omega'\right)$$

$$\qquad\quad -\cos{i'}\left[{\cal{C}}\sin\left(l'_0+2\omega'\right)-{\cal{S}}\cos\left(l'_0+2\omega'\right)\right]\big\},$$ (A.26)

$$\frac{\partial A_2}{\partial D}=2b^*_2\cos{i'}+2A^*\left[I_1\cos2\left(l'_0+\omega'\right)+I_2\sin2\left(l'_0+\omega'\right)\right],$$ (A.27)

$$\frac{\partial B_2}{\partial D}=-2a^*_2\cos{i'}+2A^*\left[I_1\sin2\left(l'_0+\omega'\right)-I_2\cos2\left(l'_0+\omega'\right)\right],$$ (A.28)

$$\frac{\partial A_3}{\partial D}=2b^*_3\cos{i'}+A^*\frac{14e'}{3}\big[I_1\cos\left(3l'_0+2\omega'\right)$$

$$\qquad\quad +I_2\sin\left(3l'_0+2\omega'\right)\big],$$ (A.29)

$$\frac{\partial B_3}{\partial D}=-2a^*_3\cos{i'}+A^*\frac{14e'}{3}\big[I_1\sin\left(3l'_0+2\omega'\right)$$

$$\qquad\quad -I_2\cos\left(3l'_0+2\omega'\right)\big],$$ (A.30)

where

$$I_1=-I\sin{i'}\sin{i_{\rm m}}\sin{u'_{\rm m}},$$ (A.31)

$$I_2=I\sin{i'}\sin{i_{\rm m}}\cos{u'_{\rm m}}.$$ (A.32)

Acknowledgements

We thank Dr. László Szabados and Mr. Szilárd Csizmadia for their comments on the manuscript, as well as Mrs. Rita Borkovits-József and Mr. Imre Barna Bíró for the linguistic corrections. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

References

 
• Borkovits, T., Csizmadia, Sz., Hegedüs, T., et al. 2002, A&A, 392, 895 In the text NASA ADS  
• Brown, E. W. 1936, MNRAS, 97, 62 In the text NASA ADS  
• Cayley, A. 1861, Mem. Roy. Astron. Soc., 29 In the text  
• Chambliss, C. R. 1992, PASP, 104, 663 In the text NASA ADS  
• Drechsel, H., Haas, S., Lorenz, R., & Mayer, P. 1994, A&A, 284, 853 In the text NASA ADS  
• Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, ApJ, 499, 853 In the text NASA ADS  
• Guilbault, P. R., Lloyd, C., & Paschke, A. 2001, IBVS, No. 5090 In the text  
• Harrington, R. S. 1968, AJ, 73, 190 In the text NASA ADS  
• Harrington, R. S. 1969, Cel. Mech., 1, 200 In the text NASA ADS  
• Khodykin, S. A., & Vedeneyev, V. G. 1997, ApJ, 475, 798 In the text NASA ADS  
• Kopal, Z. 1978, Dynamics of Close Binary Systems (D. Reidel, Dordrecht) In the text  
• Kozyreva, V. S., Zakharov, A. I., & Khaliullin, Kh. 1999, IBVS, No. 4690 In the text  
• Lacy, C. H. S., Helt, B. H., & Vaz, L. P. R. 1999, AJ, 117, 541 In the text NASA ADS  
• Lestrade, J.-F., Phillips, R. B., Hodges, M. W., & Preston, R. A. 1993, ApJ, 410, 808 In the text NASA ADS  
• Mayer, P. 1990, BAIC, 41, 231 In the text  
• Milani, A., Nobili, A. M., & Farinella, P. 1987, Non-gravitational perturbations and satellite geodesy (Adam Hilger, Bristol) In the text  
• Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1989, Numerical Recipies in Pascal (Cambridge University Press, Cambridge) In the text  
• Roy, A. E. 1988, Orbital Motion (Adam Hilger, Bristol) In the text  
• Söderhjelm, S. 1975, A&A, 42, 229 In the text NASA ADS  
• Söderhjelm, S. 1980, A&A, 89, 100 In the text NASA ADS  
• Söderhjelm, S. 1982, A&A, 107, 54 In the text NASA ADS  
• Torres, G. 2001, AJ, 121, 2227 In the text NASA ADS