EDP Sciences
Free Access
Issue
A&A
Volume 609, January 2018
Article Number A46
Number of page(s) 30
Section Catalogs and data
DOI https://doi.org/10.1051/0004-6361/201730655
Published online 05 January 2018

© ESO, 2018

1. Introduction

It is well known that substantial part of all types of star in the solar neighborhood form binary or multiple stellar systems (e.g., Abt 1983; Guinan & Engle 2006). These systems include stars of all spectral types and all stages of their evolution. Proper description of their dynamical characteristics may result in determination of their physical parameters, and thus provide clues to their formation paths (e.g., Goodwin & Kroupa 2005).

Despite the huge effort of astrophysicists in the field of stellar multiplicity in recent years, some mysteries about multiple systems still remain. For instance, recent results of thorough analysis of data from the Kepler mission have shown that there is a significant drop in population of triple systems with the period of a third component P2< 200 days (Borkovits et al. 2016). That is in accord with earlier results of Tokovinin (2014) who pointed out a similar drop of systems with the P2< 1000 days. The lull of systems at these intermediate P2 values appears to be real and subject to severe selection effects. The theoretical explanation, however, still remains unknown.

Another uncertainty persists in the problem of dynamical stability. Mardling & Aarseth (2001) derived a theoretical limit on the stability of coplanar hierarchical triple systems as (1)where m0 and m1 are masses of inner binary, m2 is a mass of the third body, e2 is an eccentricity of the outer orbit, and the parameter s = 1.8. We follow the same notation hereafter. Sterzik & Tokovinin (2002) later improved this limit according to numerical simulations and showed that exponent s in relation (1) has a different value, s = 1.35, to the best conformity between model and simulated data. However, it turned out later that there are only a few systems close to this limit and that the vast majority of observed triple systems fulfill an even stricter empirical criterion with the value of the exponent s = 3.0 (Tokovinin 2004, 2007). The lack of observed systems with high e2 might be caused by unmodeled dynamical effects (Tokovinin 2004) but generally the empirical stability criterion is still not satisfactorily explained. Nevertheless, recent analysis of 222 compact triples located in the original Kepler field of view (FOV) performed by Borkovits et al. (2016) showed that the previously derived empirical criterion may rather be a result of some observational effects and all observed triples fulfils the theoretical criterion by Sterzik & Tokovinin (2002, see Fig. 1. We point out that relation (1) is very useful for an estimation of maximal e2 from the period ratio P2/P1, even if one has no information about individual body mass, because the period ratio dependence on masses is rather weak.

In spite of the intense research, there are still only a few multiple systems with precisely determined orbits and physical parameters of all components, especially out of the Milky Way galaxy. Therefore, studying the known multiple systems, as well as a search for new ones (especially those systems, which are close to the stability limit) is crucial for obtaining sufficient statistics and comparison of theoretical simulations with observational data.

Special cases, when the inner pair of a multiple system is an eclipsing binary, are excellent laboratories for stellar evolution. Light curve analysis together with thorough study of radial velocities are able to reveal physical parameters (such as luminosities, masses, and radii) of components of the eclipsing binary (EB), as well as its orbit (Southworth 2012; Torres et al. 2010). One can also estimate some parameters of a third body which can manifest itself via various physical effects. These effects are listed below:

In the most favorable cases (e.g., high L3/ (L1 + L2) ratio or almost coplanar orbits of EB and the third component) a combination effect can occur, which makes the estimation of the third body parameters more accurate. All these effects, except for the last one, are widely used for detection of new multiple systems. Effects of orbital precession led to discovery of a new triple only in relatively few cases mentioned in Sect. 3. Timescale of the orbital precession is usually very long compared to the time span of the most extended observations that are available (about 100 yr). Consequently, only those triples with sufficiently short periods of the outer orbits can be detected in available data sets. Because of this, there are still relatively few known systems showing orbital precession, despite the fact that all binaries with P1 ≲ 1 day probably have a third companion which caused shrinkage of an initially wide orbit via combination of Kozai cycles with tidal friction (Kozai 1962; Kiseleva et al. 1998; Eggleton & Kiseleva-Eggleton 2001; Fabrycky & Tremaine 2007) and orbital precession occurs in all cases when orbits are non-coplanar. An important step forward has been made recently by Rappaport et al. (2013) and Borkovits et al. (2016), who found 42 new compact triples showing orbital precession in the original FOV of the Kepler satellite. Unfortunately there are far fewer such systems outside the Milky Way galaxy (see Sect. 3).

Only the compact triple systems with low ratio P2/P1 and short period P2 can be discovered by analysis of amplitude variation of light curve (LC) with typical timespan of several dozens of years. Interestingly, those systems fall exactly into the area of unclear lack of triples in Fig. 1, which makes this method very suitable for detection these systems and extending of their poor statistics. In this work, we aim to develop suitable methods for detecting those systems and their thorough analysis.

thumbnail Fig. 1

Left: orbital period of the third companion P2 versus orbital period of the inner binary P1 for 222 triple systems within the original Kepler field (Borkovits et al. 2016). Area with the lack of triples is highlighted with yellow. The stability limit is calculated for the parameter s = 1.8 in Eq. (1), with the assumption that m0 = m1 = m2 = 1, for the mean value of orbital eccentricity of the outer orbit e2 ⟩ ≃ 0.35. Right: eccentricity of the outer orbit as a function of the orbital period ratio. Dash-dot line is stability criterion according to Mardling & Aarseth (2001), dashed line (Sterzik & Tokovinin 2002) and red solid line is empirical criterion according to observations of Tokovinin (2007). Eighty-eight triples from Multiple Star Catalog are plotted with black crosses (Tokovinin 1997) and 222 triples in the original Kepler field are plotted with black dots (Borkovits et al. 2016).

Open with DEXTER

2. Orbital-plane precession

While the orbital-plane precession results in variation of inclination of both orbits, only changes of inclination of the eclipsing binary are observable as changes of LC amplitude. Time dependence of observed inclination of the inner binary i0 is given by the orbital geometry (Söderhjelm 1975b) (2)where I is an angle between system invariant plane and a plane tangent to the celestial sphere, i1 is an angle between the invariant plane and the orbit of the EB, t0 is a moment when the nodal line passes the plane tangent to the celestial sphere, is an angular velocity of the nodal line precession and Pnodal is a period of the nodal line precession1. We recall that the invariant plane of the system is normal to its total and conserved orbital angular momentum. It is clearly seen that the observed inclination oscillates in the interval i0 ∈ (Ii1,I + i1). Measuring time dependence of i0(t), we can use Eq. (2) to determine remaining parameters, which are deemed constant. These include i1, I, and t0. In what follows, we are mainly interested in solution of i1 and . Obviously, because of their non-linear appearance in (2), there are potentially severe correlations in their solution. We note also that Eq. (2) is invariant to the transformation i1πi1 and IπI and the orbital solution is ambiguous.

In case that eccentricity of the orbit of the inner pair is e1 ≃ 0, which is fulfilled for a vast majority of short period binaries due to circularization, the precession rate of the nodal line in the invariant system may be derived analytically with a sufficient accuracy for most of the triple systems. We have (3)where i2 is an angle between the invariant plane and the orbital plane of the third body, ak is the semimajor axis and lk is the orbital angular momentum of the respective orbits (Breiter & Vokrouhlický 2015). In these quantities, the indices k = 1 and k = 2 represent the inner and outer orbit, respectively. Angular velocity of the nodal line precession remains constant so long as the conditions e2 = const. and j = const. are fulfilled.

Modeling the variation of an EB inclination allows us to estimate the mass of the third body in the following way. Because sini2 = γsini1, the mutual inclination of both orbits j = i1 + i2 can be rewritten as (4)Eliminating a1 and a2 from Eq. (3), we obtain for γ the useful relation (5)The light curve solution together with spectra of components of the binary gives us the masses of the inner binary components m0 and m1. Therefore, γ is only a function of m2, P2, and e2, and thus . The Eq. (2), applied to the observed dependence i0(t), gives a correlated least-squares fit of parameters i1 and . The size and shape of the area with possible solutions in the parameter space depend on the data time span and also on the period Pnodal. In accordance with the functional dependence of , each solution can be transformed to the (m2,P2) parameter space, which leads to restrictions on the third body mass and its orbital period. Influence of unknown eccentricity e2 on an estimation of m2 and P2 is rather small, because even for e2 values as large as 0.2−0.3. In addition, maximal e2 can be estimated according to the stability restrictions (see the right panel of Fig. 1).

A LC solution leads to another restriction on the third body mass. A third body contributes with another light to the LC and it can be found during LC analysis in some cases. Then the third body mass can be estimated due to mass-luminosity relation. Even in cases when the third light is not detected, at least an upper limit on m2 can be estimated.

Additional restriction on the third body mass can be obtained by analysis of ETVs. Orbital period P2 is usually under the time resolution of long-term photometric surveys, because . Therefore, ETVs with period P2 cannot usually be detected. Nevertheless, dispersion in eclipse timing residual diagram at least limits the maximal amplitude of ETVs. ETVs include classical Røemer’s delay (so called light-time effect (LTE)), caused by finite travel time of light (Irwin 1952; Mayer 1990), and so called dynamical delay, which means a physical variation of P1 as a result of the gravitational interaction of components of the inner binary with the third body (Rappaport et al. 2013; Borkovits et al. 2003, 2011, 2015). Rappaport et al. (2013) showed that while Røemer’s delay is dominant for systems with a long period of the third component and a short period of the inner binary, dynamical delay is dominant for compact systems with a short period of the third component.

In this work, we deal with systems with P1 of the order of days. Therefore, for period of a third body P2< 200 days (see Fig. 7 in Rappaport et al. 2013), dynamical delay dominates and the amplitude of ETVs is given by relation (Borkovits et al. 2003) (6)which allowed us to consider a restriction on relation m2 = f(P2) because Aphys is bound from the dispersion of eclipse timing residual diagram and masses of the inner binary components can be estimated from LC solution. Conversely, for P2> 200 days classical Røemer’s delay dominates and ETVs amplitude is given as (7)where the third body inclination i3 oscilates within (Ii2,I + i2) and its maximal values are known from the time dependence of i0 (Irwin 1952; Mayer 1990)2. We note that factor 173.15 holds when ALTE and P2 are given in days and masses are in units of solar mass. Outer orbit eccentricity e2 is supposed to be unknown. However, Tokovinin (2007) showed that compact triple systems with high eccentricity of outer orbit e2 becomes unstable and there is a natural limit on e2 for each ratio of period P2/P1 (right panel of Fig. 1). According to selection effects of our methods, we can expect period ratios of analyzed systems within the interval P2/P1 ∈ (5,100) and therefore e2 ≤ 0.6 in case of all systems analyzed below (see Fig. 3 in Tokovinin 2006).

Table 1

Known inclination changing EBs in the Milky Way galaxy (42 EBs within the Kepler FOV are not included).

3. Summary of known systems showing orbital precession

Only 53 systems showing orbital precession have been discovered in the Milky Way galaxy so far. While 11 of them have been found by different observers within various fields of view in the sky (see Table 1), 42 systems have been identified within the original Kepler field (Borkovits et al. 2016). Away from the Milky Way galaxy, 17 more systems are located in the Large Magellanic Cloud (LMC; Graczyk et al. 2011; Zasche & Wolf 2013) and there is only one known system in the Small Magellanic Cloud (SMC; Pawlak et al. 2013).

Despite thorough analysis of the Kepler systems in our Galaxy, only one system outside the Milky Way galaxy has been carefully studied – MACHO 82.8043.171 (Zasche & Wolf 2013). Studying multiple systems in external galaxies and improving statistics of known systems could help to reveal inter-galactic differences in star formation, which can be affected by different metallicity or other parameters and processes (Davies et al. 2015).

Large databases of medium quality lightcurves of EBs in the LMC and SMC, from photometric surveys like OGLE (Udalski 2003; Udalski et al. 1997, 2008, 2015; Szymanski 2005) or MACHO (Bennett et al. 1993; Cook et al. 1995), allow us to develop new methods to find candidate triple systems in the Magellanic Clouds. The goal of our study is to identify new LMC and SMC compact multiple systems via detection of amplitude variations in archival LCs. These amplitude variations could be caused by an orbital precession due to the presence of a third body. Because of rather short timescale of photometric surveys (in the order of ten years) and typically long timescales of orbital precession it is possible to find only systems with short P2 and small ratio P2/P1. However, those systems should be exactly in the area with lack of triples in the P1P2 distribution (yellow area in the Fig. 1) near stability limit, which makes this simple method suitable for identifying new compact triples.

4. Methods

LC of eclipsing binaries observed by the OGLE III photometric survey have sufficient precision, low scattering and good homogeneity over the whole term of observation. This makes this database suitable for finding new compact multiple systems in the Magellanic Clouds. The OGLE III database contains LCs of 26 121 EBs in the LMC (Graczyk et al. 2011) and 6138 EBs in the SMC (Pawlak et al. 2013) (mostly in Johnson I photometric band) which are available online3 (Szymanski 2005)4.

thumbnail Fig. 2

Demonstration of Method A. Left: example LC of the system OGLE-LMC-ECL-01350 divided to three intervals. Green line represents linear fit, black lines show the limits for removal of the outliers (−3σ, +7σ). Right: phase curve for the whole LC fitted with a Fourier series of the 5th order (red line). Green line marks detected primary minimum, gray lines represent detected region of the drop in brightness.

Open with DEXTER

EB light curve amplitude variation can be caused by several effects, for example, physical variations of luminosity of the EB components, stellar spots, changes of overlapping surfaces during eclipses as a result of apsidal motion, etc. Therefore, when the maximum luminosity of both components remains constant and the orbit is circular, the only possible explanation of amplitude changes is that the inclination angle varies. Both conditions can be easily checked from phased light curve (PLC).

In principle, there may be two basic methods for detection of LC amplitude variations:

  • Method A – measurement of the eclipse depth of a PLC in differenttime intervals;

  • Method B – measurement of the amplitude of a non-phased LC in different epochs.

A typical number of data points N in an OGLE III light curve is within the interval (400,600), from which only a small fraction is located near eclipses (in the case of the most frequent EB type – detached binaries). Because of that, the number of time intervals has to be small to get a reasonable fit of minima shapes. Therefore, we had no ambition to detect other but linear amplitude changes on the timescale of OGLE III observations (8 yr), in the case of Method A.

4.1. Data preparation

Ephemerides and classification of each EB were taken from OGLE III catalog (Graczyk et al. 2011). Only those LCs that contain sufficient amount of data points (>250) were taken into account. Those LCs, where the primary minima depth were smaller than the standard deviation of magnitudes, were discarded. It should be mentioned that it obviously add an artificial selection effect to the analysis, because shallow minima may be caused by massive and luminous third companion. But this restriction is necessary to avoid false detection due to high scattering of some LCs.

For the purpose of removing outliers, each LC was fitted with a linear function f(t) and all data points outside the intervals (f(t)−3σ,f(t) + 7σ) in the case of detached binaries, (f(t)−3σ,f(t) + 6σ) in the case of semidetached binaries and (f(t)−3σ,f(t) + 3σ) in the case of overcontact binaries, were removed. We note that σ, in this case, means a standard deviation computed from the whole LC. These asymmetrical intervals were chosen with respect to different characters of LCs of individual EB types so that the whole LC (except outliers) lies inside the interval.

4.2. Method A

As the first step of this method, the whole LC of a given EB is divided into several intervals which contain approximately 150 data points. A number of intervals depends on the number of data points in the whole LC. As shown in the left panel of Fig. 2, a typical LC has been divided into three or four intervals.

In the next step, the LC in each interval is phased according to the catalog ephemerides, in order to obtain a dependence of magnitude in the Johnson I band on the phase I(Φ). The core of this method is based on fitting of primary minimum on each phased interval of the LC with a proper phenomenological function, which is able to describe well the depth of an eclipse. There is a class of mathematical functions which generate similar dependencies such as shapes of real LCs around eclipse (Andronov 2012a,b; Mikulášek et al. 2012; Mikulášek 2015). For our purposes we used the following formula (taken from Andronov 2012a) (8)where I0 is a vertical shift in magnitude, A is the depth of the eclipse, Φ0 is a phase of the time of minimum, d represents minima width and C is a parameter determines a shape of the eclipse. An example of fitted minima on each interval of the whole LC is shown in the Fig. 3.

However, Eq. (8) describes well only a small part of the PLC around minimum and each PLC had to be properly reduced before fitting. As a reasonable compromise between robustness of the algorithm and its computing speed, preliminary fitting of the PLCs with Fourier series of the fifth, fourth, and second orders were performed in the case of detached, semidetached and overcontact binaries, respectively (see the right panel of Fig. 2). The minimum of the Fourier model corresponds to the minimum of given PLC and the first maxima of this model on both sides of an eclipse define the whole part of drop in brightness5 (see the right panel of Fig. 2).

Minima depths (the parameter A from the Eq. (8)) were registered on each interval of the whole LC and finally, the linear model of its time dependence was computed. For each EB, the slope of the linear fit and R2 parameter6 were tested and if some specific values were exceeded, the EB was marked as suspicious of the inclination change. Further details about setting these values are described in Sect. 4.4.

thumbnail Fig. 3

Demonstration of Method A. Fitting of primary eclipses of OGLE-LMC-ECL-01350 with function (8) on three intervals of LC, corresponding to the left panel of Fig. 2. The increase of the minima depth is apparent.

Open with DEXTER

4.3. Method B

In the second method, the whole LC was split into eight fixed intervals with respect to the seasons of observation (see Fig. 4). In contrast with Method A, Method B is less demanding as far as the number of data points in a given interval is concerned, because it is not necessary to have a lot of points around the eclipse and that allows division of the LC to more intervals. After LC splitting, two data points7 corresponding to the lowest brightness were selected from each of the intervals of the LC and than all chosen data points were fitted with the linear model according to the σ-clipping method with a rejection of each point above the 2σ limit. As in the case of Method A, the detection criteria are the certain values of the slope of linear fit and the R2 parameter. The setting of those values is described in Sect. 4.4.

The advantage of this simple method is that there is no need to know a precise orbital period of given EB and also there is no need to delimit a part of the LC with an eclipse. That makes this method faster than the Method A. On the other hand, what is actually fitted in this case, is not exactly the minimum depth, but only local LC amplitude, which makes this method more sensitive to outliers than slower fitting minima in each part of LC.

thumbnail Fig. 4

Demonstration of Method B. Example LC of the system OGLE-LMC-ECL-01350 splitted into 8 intervals according to observational seasons. The lowest two points inside each interval are plotted with green circles, the green line shows the linear fit.

Open with DEXTER

4.4. Parameter thresholds

Specific values of slope and R2 in the case of both methods served as thresholds for preliminary distinction whether changes of inclination occur in case of given EB or not. Setting of certain thresholds on these quantities is difficult because of the high variability of LC shape and their scattering that strongly depends on the brightness of a particular EB. However, the estimation of the threshold can be made in the parameter space defined by the 17 known systems in the OGLE III LMC database (Graczyk et al. 2011).

These systems, together with two others, which have been found independently by the second author of this paper, are shown in the parameter space of both methods in the Fig. 5. However, careful examination of the previously known systems from Graczyk et al. (2011) showed that in some cases the change of the amplitude of LC could rather be an artefact than a real feature. Furthermore, in some cases, amplitude variation is so fast that a significant part of LC is completely without eclipses. As described above, our “linear” methods are usually not able to detect such systems (e.g., OGLE-LMC-ECL-17212, 17972, see Fig. 5) and parameters thresholds were set irrespective of these systems8.

thumbnail Fig. 5

All eclipsing binaries from the OGLE III LMC database (gray dots) plotted in the (slope, R2) parameter space of both methods. Nineteen previously known systems showing variations of the LC amplitude are marked with blue circles (artefacts or systems with rapid changes of amplitude) and green triangles (systems with approximately linear amplitude change on the timescale of observation). Our setting of thresholds is shown as a black line.

Open with DEXTER

Systems with approximately linear time dependence of eclipse depth on the OGLE III observational time scale (see Fig. 5) served for setting the parameters thresholds. The area of interesting systems in the parameter space was chosen to exclude as many EBs as possible without any amplitude change. Finally, our empirical thresholds for the Method A were | Δmag/Δt | ≥ 5 × 10-5mag day-1 and R2 ≥ 0.75. However, some well covered LCs with very slow changing of amplitude could exist and for these systems the slope threshold is too strong. For this reason the restrictions were alleviated, and systems with | Δmag/Δt | ≥ 1 × 10-5mag day-1 and R2 ≥ 0.90 were marked as suspicious (see the Fig. 5). For Method B the following thresholds were set: | Δmag/Δt | ≥ 2 × 10-5mag day-1 and R2 ≥ 0.50.

Table 2

Detected systems with change of LCs amplitudes in the LMC.

Table 3

Detected systems with change of LCs amplitudes in the SMC.

In order to identify systems with shorter timescales of orbital precession (Pnodal ~ 10 yr) it was neccessary to modify our methods so that they could fit time dependence of LC amplitude by polynomials of higher orders. Method A cannot be modified in this way because it demands a relatively high number of data points around eclipses. But in Method B, LC is divided to eight intervals and the number of fitted data points is mostly 16, which allowed us to apply this modification.

Considering a scattering and coverage of LCs, polynomials up to the fourth order were used. Systems fulfilling the following criteria were marked as suspected of being of higher order amplitude variations:

  • at least one higher-order polynomial fit (second, third, or fourth)was significantly better9 than a linearfit;

  • at least one higher-order polynomial fit had to have AIC (Akaike 1974), AICc (Sugiura 1978) and BIC (Schwarz 1978) information criteria value lower than a linear fit;

  • χ2of the higher-order polynomial fit was smaller than 50 to avoid systems with highly scattered data points around eclipses.

5. Results of detection of new systems with amplitude variation

Our methods described in Sect. 4 were applied to all light curves of eclipsing binaries from the OGLE III LMC+SMC. Thresholds set on the Method A and the linear regime of the Method B were exceeded in 8.1% and 3.1% of the LMC EBs, respectively. In the SMC it was in 3.9% and 2.8%. In the higher-order polynomial regime of Method B, 1.5% of LMC EBs and 1.4% of SMC EBs were marked as suspected of LC amplitude variations. All of the positive detections were subjects of thorough visual inspections and most of them were rejected as false positives.

False detections were mostly caused by too few data points around eclipses, which affects results of both methods. This is, for example, the case of OGLE-LMC-ECL-10172, which is a nice example of the system with rapid variation of LC amplitude during the term of OGLE III observations. However, the orbital period of the eclipsing binary is P1 = 1.993644 d, which makes impossible to observe eclipses in some seasons from one observing site and an artificial change of the LC amplitude occurs.

Results of Method A strongly depend on a precision of the catalog value of a binary orbital period, and also on time stability of the light curve. In some cases, the mid-eclipse times of an EB vary in time due to apsidal motion or ETVs. That leads to “blurriness” of a PLC, and results obtained by Method A may have been irrelevant and had to be rejected.

In addition, some systems had to be rejected because of variations of maximal brightness (magnitude at quadratures), which may be caused by a motion of surface spots, physical variability of one or both components of a binary, or by a drift of zero points in the measuring equipment.

After visual inspection and combining results of both methods 51 and 21 candidates remained in the LMC and SMC samples, respectively. These systems are listed in Tables 2 and 3 together with cross-identifications with the OGLE II and MACHO surveys and their lightcurves are shown in Figs. A.1 and A.2. Apart from the 13 previously known systems10, there are 38 new systems in LMC, clearly showing amplitude variations while maximal brightness remains constant. From the 21 systems detected within the SMC, only one was previously known (Pawlak et al. 2013).

6. Analysis of individual systems

Several detected systems showing LC amplitude variation were selected and analyzed thoroughly. The analysis was performed especially for systems with additional archival data (MACHO or OGLE II/IV) available.

For four selected systems, new charge-coupled device (CCD) photometry was obtained at the La Silla Observatory in Chile, with 1.54-m Danish telescope in the Johnson I photometric band to secure consistency with the OGLE data. For one system – OGLE-SMC-ECL-1532 – archival CCD frames, taken in the Johnson R photometric band with the Danish telescope, were used. For aperture photometry we developed and used Python 2.7 scripts with usage of the photutils package and differential magnitudes were obtained. The individual data (HJD vs. Δm) are listed in Table 4.

Analysis of LCs of individual systems was performed in PHOEBE (Prša & Zwitter 2005) program. Obtaining of spectra of the stars outside of the Milky Way galaxy is complicated and requires long exposure time even when using the largest telescopes. For this reason, a radial velocity curve is not available for any studied system and the mass ratio was fixed as q = 1 in the first iteration of the LC modeling. When the photometric model was inconsistent with this assumption, it was recalculated with new q value estimated from bolometric magnitudes of components. Synchronous rotation F1 = F2 = 1 of both components was assumed and limb darkening coefficients were obtained from the square-root model which is more suitable for hot stars than the logarithmic model (Diaz-Cordoves & Gimenez 1992). Bolometric albedos and gravity darkening were fixed as A1,2 = g1,2 = 1 which is fulfilled for stars with T> 7200 K whose subsurface layers are in radiative equilibrium. Without any specroscopic data for given stars, solar metallicity was assumed and fixed.

The basic model was computed on a subset with the lowest data scattering and the best coverage and on the other subsets inclinations i0 and luminosities L1,2 were fitted only. With the assumption that physical parameters of components remain constant, that led us to obtain a time dependence of binary inclination.

To improve the orbital period of the binary, primary and secondary minima were computed with the use of slightly modified AFP method (for the description of the original AFP method see Zasche et al. 2014) and eclipse timing residual diagram for each system was constructed. Modification of the original AFP method was neccessary due to the fact that LC amplitudes of our systems vary with time. Therefore, we had to fit not only a phase and magnitude shift of the model curve, but also its “contraction”, which represents amplitude variation of the LC.

The final fixed parameter was primary temperature T1 which had to be estimated from photometric indices because of a lack of another spectral information about given star in the LMC and SMC. However, correct estimation of the temperature is usually quite tricky. There are many photometric catalogs of the Magellanic Clouds with various color indices for a given star, which sometimes lead to different temperatures. In the case of the hot stars in our sample, relative differences in temperatures are up to 20%. Therefore, temperatures and masses of an individual component of an EB can be computed wrongly, which can also affect an estimated mass of the third body. But precision of estimation of the nodal period remains unaffected. For each system, all available color indices were collected and for the T1 estimate the most probable value was selected. Each color index was also corrected from an effect of interstellar reddening, according to the relation (BV)0 = f((BV),(UB)) (Johnson & Morgan 1953), a map of interstellar reddening in the LMC (Haschke et al. 2011) and mean reddening in the direction toward the SMC (Massey et al. 1995).

Table 4

Photometric observations of individual systems.

Photometric solutions of LCs of individual systems in the LMC and SMC are in Tables B.1 and B.2, respectively. In these tables, the computed and fixed parameters are marked. Luminosities in V and R Johnson photometric bands were obtained from the MACHO data. MACHO photometry was not originally obtained in standard Johnson passbands but with BMACHO and RMACHO filters instead, and it had to be transformed before the PHOEBE modeling according to the calibration relations in (Bessell & Germany 1999). Examples of LCs of each system, from which the time variations of inclination is apparent, are shown in Fig. G.1.

Precise linear ephemerides of the inner eclipsing binaries in the LMC and SMC are listed in Tables C.1 and C.2, respectively. In some cases, eclipse timing residual diagram of the EB is not linear and ETVs become apparent. For precise modeling of LCs of such systems several sets of linear ephemeris had to be calculated for each part of a given LC. In Tables C.1 and C.2, there are both the best ephemerides on whole time interval of observation and the set of ephemerides for each interval in cases when it was needed. Eclipse timing residual diagrams for every studied system is shown in Figs. 6 and F.1.

Photometric solution of each LC subset of each system enables to derive the time dependence of inclination i0(t) of every EB. In order to obtain Pnodal, physical parameters of a third body and mutual orientation of orbits, each dependence cos(i0) = f(t) was modeled with Eq. (2). The results are listed in Tables 5 and 6 together with 68% confidence intervals. Confidence interval of each fitted parameter was estimated from projection of χ2 of the fit of cos(i0) = f(t) dependence. In Fig. D.1, there is χ2 of the fit shown in cos(i1) parameter space, which provides an insight into how the parameters are limited. In most cases the model is not well limited and many possible solutions have very similar χ2. That is the reason why the results in the Tables 5 and 6 have extreme uncertainties in some cases. However, the most likely parameters of the third body can be estimated. For a given solution of cos(i0) = f(t), one value of mass of a third body m2 and its orbital period P2 can be calculated according to the Eqs. (3)–(5). In Fig. E.1, there are all m2 and P2 computed from our model with 68.3% probability in both possible orientations of orbits according to the “180°” degeneracy of the solution. The area of possible solutions is not covered homogeneously and the number of solutions in each bin is indicated. Each figure is also shown for two extremal third body orbital eccentricities e2 = 0 and e2 = 0.5 according to the argumentation in Sect. 2.

From the distribution of possible m2 and P2, angles i2 and i3 can also be estimated. For each solution of m2 and P2, one value of i2 can be computed from Eqs. (4) and (5). The most probable value and its uncertainty based on distribution, are listed in the Tables 5 and 6 for each system. With the knowledge of i2, the also “observable” inclination of the third body i3 can be computed from the Eq. (2). In Tables 5 and 6, we list their values, but usually with very large uncertainties which make it impossible to estimate real amplitude of radial velocities, which would be extremely useful for planning spectroscopic observations for confirmation of the presence of the third body.

From the LC solution, masses of EB’s components can be estimated with an assumption that both components lie on the main sequence, where the mass-luminosity ratio is relatively well defined. Therefore, analysis of the third light can lead to an estimate of the third body mass. In cases of the most of our analyzed EB the third light did not contribute to the LC with more than 1% of the total light, which is not significant with respect to a precission of photometry. In these cases, however, at least a limit on maximal possible m2 can be estimated and is also marked in Fig. E.1. As mentioned in Sect. 2, additional observational bound is given from the amplitude of ETVs which can be estimated from scattering of eclipse timing residual diagram. All presented systems are relatively compact and thus ETVs are dominated by dynamical term and its amplitude Aphys is only function of m0, m1, m2, period ratio and e2. An upper limit on the third body mass is also shown in Fig. E.1, but in cases of small m2 and short P2 the stability limit is even stricter. All limits on masses and period of third bodies are summarized in Tables 7 and 8.

Table 5

Orbital orientations of LMC systems and fitted parameters of i0 time dependence.

Table 6

Orbital orientations of SMC systems and fitted parameters of i0 time dependence.

Table 7

Masses and periods of studied triples in the LMC.

Table 8

Masses and periods of studied triples in the SMC.

6.1. OGLE-LMC-ECL-16023

OGLE-LMC-ECL-16023 (05:27:04.8669:29:01.6, I = 16.9 mag) is an overcontact binary with early spectral type components. The best estimate based on photometric indices leads to B7V spectral type and temperature of the primary component T1 = 14 000 K. Quite a large amount of photometric data from several databases including our own observations (MACHO, OGLE II/III/IV, DK154) is available, which allowed us to calculate inclination of the binary in 28 time intervals over 23 yr (see Fig. D.1). Therefore, together with relatively short period Pnodal, parametric space is rather well limited. Estimated orbital period of the third body is very short P2< 60 days, which makes this system to be very compact.

However, in eclipse timing residual diagram in the Fig. 6, ETVs with periods of approximately 6500 days are apparent. That is about two orders of magnitude longer than the expected orbital period of a third body which therefore cannot be responsible for this phenomena. Observed variation of minima timings may be due to LTE caused by another body in the system. With respect to this hypothesis, we have fitted the fourth body orbit and results are listed in Table 9, where T0 is Julian date of periastron passage of the hypothetical fourth body. Because of a large uncertainty of the third body mass, the mass function of the fourth component could not be calculated without additional spectroscopic observation.

We note that the third body mass limit, based on the limit of detectable third light in the LC, appears to be rather an approximate limit on both, third and fourth body, masses. This would mean that the maximum orbital period of the third body is even shorter than 60 days.

thumbnail Fig. 6

Eclipse timing residual (observed computed) diagrams for OGLE-LMC-ECL-16023 and 18240 with respect to linear ephemerides listed in Table C.1. Black and white points represent primary and secondary minima, respectively. The red line is the best fit of LTE and apsidal motion in the case of OGLE-LMC-ECL-16023 and 18240, respectively.

Open with DEXTER

Table 9

Fitted parameters of the LTE for OGLE-LMC-ECL-16023.

6.2. OGLE-LMC-ECL-18240

OGLE-LMC-ECL-18240 (05:31:33.6671:14:25.1, I = 17.0 mag) is a detached eclipsing binary, which was found to be also an eccentric one. For this reason, its analysis was a little different. We also included the hypothesis of apsidal motion for the detailed description of its eclipse timing residual diagram analysis (e.g., Gimenez & Garcia-Pelayo 1983). The effect of apsidal motion also affects the depths of both minima, however, the nodal precession is the most dominant contribution. Moreover, the effect of changing depth in eccentric binaries was properly modelled in our solution using the PHOEBE code. The time coverage is still rather poor and the apsidal motion slow, but the change of the periastron is apparent in the data covering about 15 yr. The eccentricity of the inner orbit of the system was found to be of about 0.2 and the apsidal motion period of 149 yr. The complete solution of apsidal motion is listed in Table 10, where Ps stands for sidereal period. Detailed photometric monitoring in the upcoming years would help us to derive its apsidal parameters with higher confidence. Once this result is achieved, the constrained apsidal motion of the inner binary and the orbital precession may provide more severe limits on the third star mass and its period.

Table 10

Fitted parameters of the apsidal motion for OGLE-LMC-ECL-18240.

7. Discussion and conclusions

Multiple stellar systems are important astrophysical laboratories which could help us to understand general mechanisms of mutual N-body dynamic interactions between components, as well as the process of their formation. However, the number of relatively well studied multiple systems still remains low, especially the compact ones manifesting orbital precession. Besides that, results from the Kepler mission show an interesting distribution of orbital periods of triple systems, which is not explained theoretically and more investigation is needed.

In this work, we focused on changing of inclination of eclipsing binaries, and developed new methods, which appear to be suitable for detection of new triples with a small P2/P1 ratio. The presented methods led to an identification of 58 new compact triple candidates within the LMC and SMC which is, together with 14 previously known systems, the largest published sample of inclination changing compact triple candidates out of the Milky Way galaxy.

thumbnail Fig. 7

Area with the lack of triples in P1P2 distribution with the systems studied in this work. As only upper limit on P2 has been estimated, each system is depicted like an arrow with upper limit. Triples found within the Kepler field are marked with a gray dots.

Open with DEXTER

Eight of detected systems were studied thoroughly to determine the basic physical parameters of the eclipsing pair, the third component, and mutual orientation of the orbits. Unfortunately, we found that for precise determination of mutual inclinations a time base of given observations (more than 20 yr in some cases) is still too short and obtained angles were computed with very large uncertainties. In some cases also Pnodal were computed with large uncertainty, but despite that the systems in our sample still belongs to those with rather small Pnodal (compare Tables 7 and 8 in this work with Table 10 in Borkovits et al. 2016, with a large sample of compact triples discovered by Kepler mission). Our analysis also led to restrictions on the upper limit of possible orbital period P2 of the third body. The distribution of the P1 and P2 with the results of our study are shown in the Fig. 7. One can clearly see that our sample is almost completely located within the area of the lack of triples, as reported by Borkovits et al. (2016). Moreover, periods of some systems in our sample might be close to the limit of stability, which is not determined unambiguously (see Mardling & Aarseth 2001; Sterzik & Tokovinin 2002; Tokovinin 2004, 2007). But new observations with longer time bases, in the ideal case both photometric and spectroscopic, are needed to obtain parameters for all identified triples and to improve the precision of the determined parameters for the eight analyzed systems.

It should be noted that the upper limits on m2 and consequently P2 are estimated with relatively rough assumptions based on absence of third light in the light curve solutions and the real upper limits might be slightly different. However, it cannot change the fact that systems presented in this paper are most probably within the region with the lack of triples in the P1P2 distribution, and that makes listed systems as perfect targets for a campaign of photometric observations targeting on minima timing with cadence of the order of weeks, which should lead to precise determination of the third body orbital period. Absence of a third light in the light curve solution for all studied systems could seem quite interesting, because it also means that m2<m0 + m1 but in the Milky Way galaxy the third body mass m2 tends to be similar to the mass of eclipsing pair m0 + m1 and for 81% of triples m2/ (m0 + m1) ratio is greater than 0.2 (Tokovinin 2008; Correia et al. 2006). Comparison of distributions of triple component masses between the Milky Way galaxy and the Magellanic Clouds with different metalicities could be very interesting and useful for theory of multiple stellar systems formation. But in this case the third light absence is most probably result of selection effects of our methods.


1

Description of all symbols used is summarized in Table H.1.

2

We note that the original Irwin’s formula for the Røemer’s delay contains an extra term due to unusual convention of coordinate system. This term is constant as long as the elements of the orbit remains constant, which is clearly not the case of the systems analyzed in this work (see footnoote 2 in Borkovits et al. 2016). However, the amplitude ALTE remains the same for the Irwin’s as well as for the modern convention.

3

ftp.astrouw.edu.pl

4

We note that at the time of the first part of our analysis, the OGLE IV data had not yet been released.

5

An optimal solution was found by repetition of the fitting procedure with several multiples of width obtained from the Fourier model.

6

Defined as , where yi is the observed data value, fi is the predicted value from the fit, is the mean of the observed data. is the weight for each data point, where σi is uncertainty of each minima depth determined by the Levenberg-Marquardt curve-fitting algorithm.

7

Selection of more then one data points with extremal value is for better robustness of the method.

8

However, several such systems have still been detected with described linear methods or with the modified method described below (see Figs. A.1 and A.2 showing all detected systems).

9

Meaning that the reduced χ2 of polynomial fit was at least two times smaller than χ2 of the linear fit.

10

Only 13 from a whole sample of 17 known systems in LMC are detectable with our methods.

Acknowledgments

This work was supported by the Czech Science Foundation grant No. GA15-02112S, and also by the grant MSMT INGO II LG15010. We are also grateful to the ESO team at the La Silla Observatory for their help in maintaining and operating the Danish telescope. We do thank the macho and ogle teams for making all of the observations easily public available. We would like to thank also the EROS-1 team – Jean-Baptiste Marquette, Philippe Schwemling and Marc Moniez, who kindly provided us archival photometric data. Marek Skarka acknowledges the support of the postdoctoral fellowship program of the Hungarian Academy of Sciences at the Konkoly Observatory as a host institution and the financial support of the Hungarian NKFIH Grant K-115709. This research was carried out under the project CEITEC 2020 (LQ1601) with financial support from the Ministry of Education, Youth and Sports of the Czech Republic under the National Sustainability Programme II. This research has made use of the SIMBAD and VIZIER databases, operated at CDS, Strasbourg, France and of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Light curves of EBs with amplitude variation located in the LMC and SMC

thumbnail Fig. A.1

Light curves of EBs with amplitude variation from the OGLE III LMC database.

Open with DEXTER

thumbnail Fig. A.2

Light curves of EBs with amplitude variation from the OGLE III SMC database.

Open with DEXTER

Appendix B: Photometric solutions for selected eclipsing binaries in the LMC and SMC

Table B.1

Eclipsing binaries in the LMC.

Table B.2

Eclipsing binaries in the SMC.

Appendix C: Ephemerides for analyzed eclipsing binaries

Table C.1

Ephemerides for the LMC eclipsing binaries.

Table C.2

Ephemerides for the SMC eclipsing binaries.

Appendix D: Solutions of the time dependencies of inclination for selected systems

thumbnail Fig. D.1

Left: projections of χ2 of fitted dependence cos(i0) = f(t) to cos(i1) plane. Confidence intervals 68.3%, 90.0%, and 99.0% are indicated with gray lines. White points show the best solutions whose parameters are listed in Tables 5 and 6. Right: the best solution of dependence cos(i0) = f(t) according to relation (2).

Open with DEXTER

Appendix E: Possible masses and periods of third bodies

thumbnail Fig. E.1

Possible masses and periods of third bodies. All solutions from 68.3 % area in the Fig. D.1 for possible orientations i1πi1 and IπI are shown. Logarithmic scale of a number of solutions (NS) is shown in blue. The orange line indicates the maximal third body mass according to the limit of a third light from the LC solution. Red solid line indicates the stability limit according to Mardling & Aarseth (2001), dashed line indicates the stability limit according to Sterzik & Tokovinin (2002) and dot-dashed line according to Tokovinin (2004). The blue line shows restrictions on ETVs. All figures are plotted for the extremal third body eccentricities e2 = 0 and e2 = 0.5.

Open with DEXTER

Appendix F: Eclipse timing residual diagrams

thumbnail Fig. F.1

Eclipse timing residual (observed computed) diagrams for selected systems with respect to mean linear ephemerides listed in Tables C.1 and C.2. Black and white points represent primary and secondary minima, respectively.

Open with DEXTER

Appendix G: Light curves of selected systems

thumbnail Fig. G.1

Light curves of selected eclipsing binaries. Only several modeled LCs for each system are shown. For each LC, survey designation, photometric band, and central time of given interval are listed.

Open with DEXTER

Appendix H: Definition of symbols

Table H.1

Definition of symbols.

All Tables

Table 1

Known inclination changing EBs in the Milky Way galaxy (42 EBs within the Kepler FOV are not included).

Table 2

Detected systems with change of LCs amplitudes in the LMC.

Table 3

Detected systems with change of LCs amplitudes in the SMC.

Table 4

Photometric observations of individual systems.

Table 5

Orbital orientations of LMC systems and fitted parameters of i0 time dependence.

Table 6

Orbital orientations of SMC systems and fitted parameters of i0 time dependence.

Table 7

Masses and periods of studied triples in the LMC.

Table 8

Masses and periods of studied triples in the SMC.

Table 9

Fitted parameters of the LTE for OGLE-LMC-ECL-16023.

Table 10

Fitted parameters of the apsidal motion for OGLE-LMC-ECL-18240.

Table B.1

Eclipsing binaries in the LMC.

Table B.2

Eclipsing binaries in the SMC.

Table C.1

Ephemerides for the LMC eclipsing binaries.

Table C.2

Ephemerides for the SMC eclipsing binaries.

Table H.1

Definition of symbols.

All Figures

thumbnail Fig. 1

Left: orbital period of the third companion P2 versus orbital period of the inner binary P1 for 222 triple systems within the original Kepler field (Borkovits et al. 2016). Area with the lack of triples is highlighted with yellow. The stability limit is calculated for the parameter s = 1.8 in Eq. (1), with the assumption that m0 = m1 = m2 = 1, for the mean value of orbital eccentricity of the outer orbit e2 ⟩ ≃ 0.35. Right: eccentricity of the outer orbit as a function of the orbital period ratio. Dash-dot line is stability criterion according to Mardling & Aarseth (2001), dashed line (Sterzik & Tokovinin 2002) and red solid line is empirical criterion according to observations of Tokovinin (2007). Eighty-eight triples from Multiple Star Catalog are plotted with black crosses (Tokovinin 1997) and 222 triples in the original Kepler field are plotted with black dots (Borkovits et al. 2016).

Open with DEXTER
In the text
thumbnail Fig. 2

Demonstration of Method A. Left: example LC of the system OGLE-LMC-ECL-01350 divided to three intervals. Green line represents linear fit, black lines show the limits for removal of the outliers (−3σ, +7σ). Right: phase curve for the whole LC fitted with a Fourier series of the 5th order (red line). Green line marks detected primary minimum, gray lines represent detected region of the drop in brightness.

Open with DEXTER
In the text
thumbnail Fig. 3

Demonstration of Method A. Fitting of primary eclipses of OGLE-LMC-ECL-01350 with function (8) on three intervals of LC, corresponding to the left panel of Fig. 2. The increase of the minima depth is apparent.

Open with DEXTER
In the text
thumbnail Fig. 4

Demonstration of Method B. Example LC of the system OGLE-LMC-ECL-01350 splitted into 8 intervals according to observational seasons. The lowest two points inside each interval are plotted with green circles, the green line shows the linear fit.

Open with DEXTER
In the text
thumbnail Fig. 5

All eclipsing binaries from the OGLE III LMC database (gray dots) plotted in the (slope, R2) parameter space of both methods. Nineteen previously known systems showing variations of the LC amplitude are marked with blue circles (artefacts or systems with rapid changes of amplitude) and green triangles (systems with approximately linear amplitude change on the timescale of observation). Our setting of thresholds is shown as a black line.

Open with DEXTER
In the text
thumbnail Fig. 6

Eclipse timing residual (observed computed) diagrams for OGLE-LMC-ECL-16023 and 18240 with respect to linear ephemerides listed in Table C.1. Black and white points represent primary and secondary minima, respectively. The red line is the best fit of LTE and apsidal motion in the case of OGLE-LMC-ECL-16023 and 18240, respectively.

Open with DEXTER
In the text
thumbnail Fig. 7

Area with the lack of triples in P1P2 distribution with the systems studied in this work. As only upper limit on P2 has been estimated, each system is depicted like an arrow with upper limit. Triples found within the Kepler field are marked with a gray dots.

Open with DEXTER
In the text
thumbnail Fig. A.1

Light curves of EBs with amplitude variation from the OGLE III LMC database.

Open with DEXTER
In the text
thumbnail Fig. A.2

Light curves of EBs with amplitude variation from the OGLE III SMC database.

Open with DEXTER
In the text
thumbnail Fig. D.1

Left: projections of χ2 of fitted dependence cos(i0) = f(t) to cos(i1) plane. Confidence intervals 68.3%, 90.0%, and 99.0% are indicated with gray lines. White points show the best solutions whose parameters are listed in Tables 5 and 6. Right: the best solution of dependence cos(i0) = f(t) according to relation (2).

Open with DEXTER
In the text
thumbnail Fig. E.1

Possible masses and periods of third bodies. All solutions from 68.3 % area in the Fig. D.1 for possible orientations i1πi1 and IπI are shown. Logarithmic scale of a number of solutions (NS) is shown in blue. The orange line indicates the maximal third body mass according to the limit of a third light from the LC solution. Red solid line indicates the stability limit according to Mardling & Aarseth (2001), dashed line indicates the stability limit according to Sterzik & Tokovinin (2002) and dot-dashed line according to Tokovinin (2004). The blue line shows restrictions on ETVs. All figures are plotted for the extremal third body eccentricities e2 = 0 and e2 = 0.5.

Open with DEXTER
In the text
thumbnail Fig. F.1

Eclipse timing residual (observed computed) diagrams for selected systems with respect to mean linear ephemerides listed in Tables C.1 and C.2. Black and white points represent primary and secondary minima, respectively.

Open with DEXTER
In the text
thumbnail Fig. G.1

Light curves of selected eclipsing binaries. Only several modeled LCs for each system are shown. For each LC, survey designation, photometric band, and central time of given interval are listed.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.