© 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 P₂< 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 P₂< 1000 days. The lull of systems at these intermediate P₂ 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 ${\mathit{P}}_{\mathrm{2}}\gtrsim \mathrm{4.7}(\frac{{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathrm{1}}}{)}^{\mathrm{1}\mathit{/}\mathrm{10}}\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{e}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{5}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{2}}{\mathrm{)}}^{\mathit{S}}}{\mathit{P}}_{\mathrm{1}}\mathit{,}$(1)where m₀ and m₁ are masses of inner binary, m₂ is a mass of the third body, e₂ 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 e₂ 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 e₂ from the period ratio P₂/P₁, 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:

• eclipse time variations (ETVs; Irwin 1952; Mayer 1990; Borkovitset al. 2016);

• presence of a third light in the LC solution or a third spectrum in the overall spectrum;

• visual resolution of the third body and its orbit;

• variation of γ-velocity of the binary;

• additional eclipses in a few rare cases (Carter et al. 2011; Borkovits et al. 2015; Alonso et al. 2015);

• variation of LC amplitude as a result of orbital plane precession (Söderhjelm 1975b; Borkovits et al. 2016, 2011; Breiter & Vokrouhlický 2015).

In the most favorable cases (e.g., high L₃/ (L₁ + L₂) 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 ${\mathit{P}}_{\mathrm{nodal}}\mathrm{~}{\mathit{P}}_{\mathrm{2}}^{\mathrm{2}}\mathit{/}{\mathit{P}}_{\mathrm{1}}$ 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 P₁ ≲ 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 P₂/P₁ and short period P₂ 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.

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 i₀ is given by the orbital geometry (Söderhjelm 1975b) $\begin{array}{ccc}& & \cos {\mathit{i}}_{\mathrm{0}}\mathrm{=}\cos \mathit{I}\cos {\mathit{i}}_{\mathrm{1}}\mathrm{\pm }\sin \mathit{I}\sin {\mathit{i}}_{\mathrm{1}}\cos \mathrm{\Omega ̇}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\mathit{,}\\ & & \mathrm{\Omega ̇}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{{\mathit{P}}_{\mathrm{nodal}}}\mathit{,}\end{array}$(2)where I is an angle between system invariant plane and a plane tangent to the celestial sphere, i₁ is an angle between the invariant plane and the orbit of the EB, t₀ 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 P_nodal is a period of the nodal line precession¹. 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 i₀ ∈ (I−i₁,I + i₁). Measuring time dependence of i₀(t), we can use Eq. (2) to determine remaining parameters, which are deemed constant. These include i₁, I, and t₀. In what follows, we are mainly interested in solution of i₁ 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 i₁ → π−i₁ and I → π−I and the orbital solution is ambiguous.

In case that eccentricity of the orbit of the inner pair is e₁ ≃ 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 $\begin{array}{ccc}& & \mathrm{\Omega ̇}\mathrm{\simeq }\frac{\mathrm{2}\mathit{\pi }}{{\mathit{P}}_{\mathrm{1}}}\frac{\mathrm{3}}{\mathrm{4}{\mathit{\eta }}_{\mathrm{2}}^{\mathrm{3}}}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{M}}_{\mathrm{2}}}(\frac{{\mathit{P}}_{\mathrm{1}}}{{\mathit{P}}_{\mathrm{2}}}{)}^{\mathrm{2}}\cos \mathit{j}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\gamma }}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{\gamma }\cos \mathit{j}}\mathit{,}\\ & & \mathit{\gamma }\mathrm{=}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{2}}}\frac{{\mathit{l}}_{\mathrm{1}}}{{\mathit{l}}_{\mathrm{2}}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{2}}}\frac{{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{2}}{\mathit{P}}_{\mathrm{2}}{\mathit{a}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{2}}{\mathit{M}}_{\mathrm{1}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{1}}{\mathit{a}}_{\mathrm{2}}^{\mathrm{2}}}\mathit{,}{\mathit{\eta }}_{\mathrm{2}}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{2}}^{\mathrm{2}}}\mathit{,}\\ & & {\mathit{M}}_{\mathrm{1}}\mathrm{=}{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{M}}_{\mathrm{2}}\mathrm{=}{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}\mathit{,j}\mathrm{=}{\mathit{i}}_{\mathrm{1}}\mathrm{+}{\mathit{i}}_{\mathrm{2}}\mathit{,}\end{array}$(3)where i₂ is an angle between the invariant plane and the orbital plane of the third body, a_k is the semimajor axis and l_k 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 e₂ = 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 sini₂ = γsini₁, the mutual inclination of both orbits j = i₁ + i₂ can be rewritten as $\mathit{j}\mathrm{=}{\mathit{i}}_{\mathrm{1}}\mathrm{+}\arcsin \mathrm{\left(}\mathit{\gamma }\sin {\mathit{i}}_{\mathrm{1}}\mathrm{\right)}\mathit{.}$(4)Eliminating a₁ and a₂ from Eq. (3), we obtain for γ the useful relation $\mathit{\gamma }\mathrm{=}\frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{2}}}\frac{{\mathit{m}}_{\mathrm{0}}{\mathit{m}}_{\mathrm{1}}}{{\mathit{M}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}(\frac{{\mathit{P}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{2}}}{{\mathit{P}}_{\mathrm{2}}{\mathit{M}}_{\mathrm{1}}}{)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(5)The light curve solution together with spectra of components of the binary gives us the masses of the inner binary components m₀ and m₁. Therefore, γ is only a function of m₂, P₂, and e₂, and thus . The Eq. (2), applied to the observed dependence i₀(t), gives a correlated least-squares fit of parameters i₁ 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 P_nodal. In accordance with the functional dependence of , each solution can be transformed to the (m₂,P₂) parameter space, which leads to restrictions on the third body mass and its orbital period. Influence of unknown eccentricity e₂ on an estimation of m₂ and P₂ is rather small, because $\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{2}}^{\mathrm{2}}}\mathrm{\simeq }\mathrm{1}$ even for e₂ values as large as 0.2−0.3. In addition, maximal e₂ 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 m₂ can be estimated.

Additional restriction on the third body mass can be obtained by analysis of ETVs. Orbital period P₂ is usually under the time resolution of long-term photometric surveys, because ${\mathit{P}}_{\mathrm{2}}\mathrm{~}\sqrt{{\mathit{P}}_{\mathrm{1}}{\mathit{P}}_{\mathrm{nodal}}}$. Therefore, ETVs with period P₂ 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 P₁ 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 P₁ of the order of days. Therefore, for period of a third body P₂< 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) ${\mathit{A}}_{\mathrm{phys}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{8}\mathit{\pi }}\frac{{\mathit{m}}_{\mathrm{2}}}{{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}}\frac{{\mathit{P}}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{P}}_{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathit{,}$(6)which allowed us to consider a restriction on relation m₂ = f(P₂) because A_phys 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 P₂> 200 days classical Røemer’s delay dominates and ETVs amplitude is given as ${\mathit{A}}_{\mathrm{LTE}}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{2}}\sin {\mathit{i}}_{\mathrm{3}}}{\mathrm{(}{\mathit{m}}_{\mathrm{0}}\mathrm{+}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}\mathit{/}\mathrm{3}}}\frac{{\mathit{P}}_{\mathrm{2}}^{\mathrm{2}\mathit{/}\mathrm{3}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{2}}^{\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\omega }}_{\mathrm{2}}}}{\mathrm{173.15}}\mathit{,}$(7)where the third body inclination i₃ oscilates within (I−i₂,I + i₂) and its maximal values are known from the time dependence of i₀ (Irwin 1952; Mayer 1990)². We note that factor 173.15 holds when A_LTE and P₂ are given in days and masses are in units of solar mass. Outer orbit eccentricity e₂ is supposed to be unknown. However, Tokovinin (2007) showed that compact triple systems with high eccentricity of outer orbit e₂ becomes unstable and there is a natural limit on e₂ for each ratio of period P₂/P₁ (right panel of Fig. 1). According to selection effects of our methods, we can expect period ratios of analyzed systems within the interval P₂/P₁ ∈ (5,100) and therefore e₂ ≤ 0.6 in case of all systems analyzed below (see Fig. 3 in Tokovinin 2006).

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 P₂ and small ratio P₂/P₁. However, those systems should be exactly in the area with lack of triples in the P₁−P₂ 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 online³ (Szymanski 2005)⁴.

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) $\mathit{I}\mathrm{\left(}\mathrm{\Phi }\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{0}}\mathrm{+}\mathit{A}{\mathrm{\{}\mathrm{1}\mathrm{-}{\left[\mathrm{(}\mathrm{\Phi }\mathrm{-}{\mathrm{\Phi }}_{\mathrm{0}}\mathrm{)}\mathit{/}\mathit{d}\right]}^{\mathrm{2}}\mathrm{\}}}^{\mathit{C}}\mathit{,}$(8)where I₀ is a vertical shift in magnitude, A is the depth of the eclipse, Φ₀ 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 brightness⁵ (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 R² parameter⁶ 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.

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 points⁷ 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 R² 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.

4.4. Parameter thresholds

Specific values of slope and R² 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 systems⁸.

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 R² ≥ 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 R² ≥ 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 R² ≥ 0.50.

In order to identify systems with shorter timescales of orbital precession (P_nodal ~ 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 better⁹ 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;

• χ²of 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 P₁ = 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 systems¹⁰, 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 F₁ = F₂ = 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 A_1,2 = g_1,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 i₀ and luminosities L_1,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 T₁ 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 T₁ estimate the most probable value was selected. Each color index was also corrected from an effect of interstellar reddening, according to the relation (B−V)₀ = f((B−V),(U−B)) (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).

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 B_MACHO and R_MACHO 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 i₀(t) of every EB. In order to obtain P_nodal, physical parameters of a third body and mutual orientation of orbits, each dependence cos(i₀) = 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 χ² of the fit of cos(i₀) = f(t) dependence. In Fig. D.1, there is χ² of the fit shown in cos(i₁) – 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 χ². 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(i₀) = f(t), one value of mass of a third body m₂ and its orbital period P₂ can be calculated according to the Eqs. (3)–(5). In Fig. E.1, there are all m₂ and P₂ 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 e₂ = 0 and e₂ = 0.5 according to the argumentation in Sect. 2.

From the distribution of possible m₂ and P₂, angles i₂ and i₃ can also be estimated. For each solution of m₂ and P₂, one value of i₂ 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 i₂, the also “observable” inclination of the third body i₃ 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 m₂ 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 A_phys is only function of m₀, m₁, m₂, period ratio and e₂. An upper limit on the third body mass is also shown in Fig. E.1, but in cases of small m₂ and short P₂ the stability limit is even stricter. All limits on masses and period of third bodies are summarized in Tables 7 and 8.

6.1. OGLE-LMC-ECL-16023

OGLE-LMC-ECL-16023 (05:27:04.86−69: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 T₁ = 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 P_nodal, parametric space is rather well limited. Estimated orbital period of the third body is very short P₂< 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 T₀ 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.

6.2. OGLE-LMC-ECL-18240

OGLE-LMC-ECL-18240 (05:31:33.66−71: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 P_s 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.

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 P₂/P₁ 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.

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 P_nodal were computed with large uncertainty, but despite that the systems in our sample still belongs to those with rather small P_nodal (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 P₂ of the third body. The distribution of the P₁ and P₂ 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 m₂ and consequently P₂ 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 P₁−P₂ 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 m₂<m₀ + m₁ but in the Milky Way galaxy the third body mass m₂ tends to be similar to the mass of eclipsing pair m₀ + m₁ and for 81% of triples m₂/ (m₀ + m₁) 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.

---

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

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

Appendix C: Ephemerides for analyzed eclipsing binaries

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

Appendix E: Possible masses and periods of third bodies

Appendix F: Eclipse timing residual diagrams

Appendix G: Light curves of selected systems

Appendix H: Definition of symbols

All Tables

All Figures