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. Borkovits1 - B. Érdi2 - E. Forgács-Dajka2 - T. Kovács2,

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')2 (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'2/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', , , , where a12 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:

 (1)

where G denotes the gravitational constant, m3 is the mass of the tertiary, while 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:

 (2)

where M12 means the total mass of the binary, and denote the absolute value of the first two Jacobian position-vectors, (e.g. means the separation of the members of the eclipsing pair, while is the distance between the centre of mass of the binary and the distant third companion), while stands for the direction cosine between and . 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 , , , respectively. The direction cosines between the vector and the axes (as it can be seen in Fig. 1) are as follows:

 (3) (4) (5)

where w and w' denote the true longitude of the secondary and the tertiary measured from the intersection of the orbits, and 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 are

 (6) (7) (8) (9)

where

 (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 transversal component, while 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

 (11)

function, where

 (12)

furthermore v' denotes the true anomaly of the outer body, is the argument of the periastron of the binary's orbit around the centre of mass of the triple system, and M123 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:

 (13) (14) (15) (16)

where the phase angles are
 = (17) = (18)

(In the above expressions refers to the true longitude of the intersection of the two orbits measured along the inner orbit, while 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

 (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 , 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:

 (20)

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

 (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 means the special angular momentum of the binary, e.g. the length of vector , and not the velocity of light.)

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

 (22)

where

 (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 t0 epoch can be determined as

 (24)

Consequently, if the (or ) versus (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 as

 (25)

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

 (26)

while

 (27)

As far as the perturbations are small in the mean motion the relation between the instantaneous period and the Keplerian period P0 of the unperturbed motion can be written as

 (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 component on the instantaneous period. It can be easily seen, that

 (29)

The second derivative of the part of the true longitude has the following form (cf. Milani et al. 1987, Chap. 3.2):

 (30)

where

 (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 , and so it can be shown that the two terms on the rhs of (30) may have the same order. Furthermore, since both and have the same order of magnitude as the small variation in the 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 term of Eq. (30), as in the first expression the multiplicator in the amplitude of cancels the denominator.

Let us define the expessions

 (32) (33)

respectively. Then

 (34)

For the evaluation of the first term we have to express the true anomaly of the third component by . 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:

 (35)

where l'0 is the mean anomaly at the epoch t0. Now using the expressions

 (36) (37)

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

 (38)

Our next task is the calculation of (33). The dependence of the integrand on can be written as
 = (39)

where the further integrands can be evaluated with good approximation as

 (40)

and

 (41)

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

 (42)

where 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 function. As it is well known

 (43)

Since has itself the same order as the components of  , the other quantities in (43) can be treated as constants. Consequently, our approximation for is the following:

 = (44)

Let us turn back to the expression (28) of . It can be seen easily that all of the above calculated perturbative terms have the order of (P/P')2, 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)

 (45)

Consequently, the amplitude of these integrated expressions is multiplied by . 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

 (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:

 = (47)

where

 (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 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 evidently has to be measured from the plane of the sky. However, the two meanings of the 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, , consequently the terms multiplied by  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.

 System m1 m2 P e i 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

 System f(m3) P' e' "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.)

 No. m3 i' 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. . 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:

 (49)

where N is the cycle number, and

 (50)

while
 (51) (52)

where

 (53)

furthermore,
 (54) (55)

and in the latter expressions Jk 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:

 (56) (57) (58) (59) (60) (61)

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

 (62)

where
 (63) (64) (65)

We omitted terms which are multiplied by , as in eclipsing systems the observable inclination of the binary is usually close to , consequently these terms give only a minor contribution. (E.g. even for , .) A comparison between (62) and (49) reveals that in this case the coefficients are as follows:

 (66) (67) (68) (69) (70) (71)

where the amplitude A* is

 (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' a12 i' D m3 MHJD 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' a12 i' D m3
MHJD 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 (m3), the observable inclination of the wide orbit (i'), and the mutual inclination (). It is evident that the dependence on m3 appears in the amplitude A*, while the effect of the inclinations manifests through the following equations of spherical triangles:

 (73) (74) (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: c0 (a zero-point correction), P, A (the amplitude of the light-time terms), e', , l'0, i', and instead of the mutual inclination, . (Here we note, that although P is an adjusted parameter, de fundamental frequency 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 m3 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' a12 i' D m3 MHJD 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' a12 i' D m3 MHJD 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', , than it was obtained from the simple light-time (L) fit. In these runs the value also improved with respect to the corresponding L-fits;
• Significant exceptions arise in fits A3, A6, A12, e.g. at those runs, where . In these cases , while , which is also zero for run A3, while and for runs A6 and A12, respectively. Consequently, in these cases only the -term, and the coefficients a*1, b*1 have significant contributions. On the other hand, in run A15, where was zero a good fit was also found. In this case already the term is the dominant. (Here we note, that in this integration the mutual inclination was , which is very close to that "critical'' value where the -term disappears.) It is interesting, that despite the weaker fits in these cases, the , 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 m3, 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 . 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 , while for IU Aur the same value is .

The three presently known closest eclipsing triple systems are  Tau ( ), DM Per ( ), and VV Ori ( ) (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  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: c0, P, , e', , l'0, i', . Here we list the analytical form of the partial derivatives of the Fourier-coefficients with respect to these parameters. The ak, bk 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.1) (A.2) (A.3) (A.4) (A.5) (A.6) (A.7) (A.8) (A.9) (A.10) (A.11) (A.12) (A.13) (A.14) (A.15) (A.16) (A.17) (A.18) (A.19) (A.20) (A.21) (A.22) (A.23) (A.24) (A.25) (A.26) (A.27) (A.28) (A.29) (A.30)

where
 (A.31) (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.