© ESO, 2017

1. Introduction

The interest of binary and multiple systems spans various fields of astrophysics, including stellar formation (initial conditions and formation processes), stellar physics and evolution (accurate stellar parameters determinations and comparison with model predictions), galactic and extra-galactic distance determinations (e.g. Southworth 2012), and cosmology (e.g. type Ia supernovae). Until the end of the twentieth century, binary systems were almost exclusively studied on a case by case basis. The advent of large scale multi-epoch photometric surveys almost three decades ago with the “Expérience pour la recherche d’objets sombres” (EROS-1, 1990−1995; Aubourg et al. 1993; Renault et al. 1998), the “Massive compact halo object” experiment (MACHO, 1992−1999; Alcock et al. 1997), and the “Optical gravitational lensing experiment” (OGLE-I, 1992−1995; Udalski et al. 1992) opened the door to studies based on large databases containing thousands to tens of thousands of eclipsing binaries (EBs) in various stellar populations. Catalogues of EB light curves (LCs) have been published, for example, by the OGLE-III project for the Large Magellanic Cloud (LMC; 26 121 sources, Graczyk et al. 2011), for the Small Magellanic Cloud (SMC; 6138 sources, Pawlak et al. 2013), and for the galactic disk fields (11 589 sources, Pietrukowicz et al. 2013). And very recently, the OGLE team updated the list of EBs in the Magellanic Clouds with new results from the OGLE-IV project (40 204 sources in the LMC and 8401 sources in the SMC, Pawlak et al. 2016).

Another new leap will soon be achieved with ongoing and future very large scale multi-epoch surveys that will further increase the number of EBs as well as the level of completeness to an unprecedented degree. One of those surveys is the European Gaia space mission (Perryman et al. 2001; Gaia Collaboration 2016b), launched in December 2013, the primary aim of which is to determine the three-dimensional positions of over one billion stars in the Galaxy. The preliminary data published in Gaia Data Release 1 (Gaia Collaboration 2016a) reveals the great potential of the Gaia mission in terms of astrometry, photometry, and number of sources surveyed. With its combination of all-sky coverage, of multi-epoch white-band photometry (a mean of ~70 photometry transits per source is expected during its five-year mission), of simultaneous multi-epoch spectro-photometry in blue and red bands, of simultaneous multi-epoch radial velocities and basic astrophysical parameter determinations for the brightest stars, all this in addition to the parallax determinations, the Gaia mission is a golden mine for all fields of astrophysics. In particular, the mission is expected to record the light curves of between half and several million EBs (e.g. Dischler & Söderhjelm 2005; Eyer et al. 2013). Another example of a promising multi-epoch large scale survey is the photometric Large Synoptic Survey Telescope (LSST) project planned to enter science operations in 2022 (Ivezic et al. 2008).

The analysis of hundreds of thousands to millions light (and radial velocity when available) curves from those large scale surveys presents new challenges and requires the development of automated techniques. Within the Gaia Data Processing and Analysis Consortium, our Geneva-led team is responsible for the detection, characterization and classification of variable objects in general (Eyer et al. 2017), and of EBs in particular. To achieve these goals for hundreds of thousands of EB LCs, novel processing and analysis techniques are being explored. The results of these investigations are applied to existing surveys of EBs and simulated Gaia data, and make the object of these series of papers. Two classification techniques have already been explored, using both existing surveys and Gaia simulated data (Süveges et al. 2017; Kochoska et al. 2017). Here, we present a method to characterize eclipse and inter-eclipse properties based on the geometry of EB folded LCs (FLCs).

The study has two goals. The first goal is to provide a set of EB parameters that allows to study the ensemble of EBs on a statistical ground without the need to model the binary system. The second goal is to identify within the large data set binary systems with unexpected properties that could reveal the existence of new configurations. The procedure is based on modeling the geometry of EB LCs using Gaussian functions to model the eclipses and a cosine function to model ellipsoidal variability, if present. The models, which we generically refer to in this paper as the two-Gaussian models, whether they actually contain two, one, or no Gaussian, are described in Sect. 2. The procedure is applied in Sect. 3 to the set of EBs from the LMC identified by the OGLE-III survey. The capability of the two-Gaussian models to achieve the two goals is then addressed in Sect. 4. Conclusions are drawn in Sect. 5.

A table summarizing the EB parameters derived in this study for the OGLE-III EBs of the LMC is made available at the CDS. Its content is described in Appendix A.

2. Two-Gaussian models

We present a description of the geometrical models used to characterize the LCs of EBs (Sect. 2.1), the model computation procedure (Sect. 2.2), and the best-model selection criterion (Sect. 2.3).

2.1. Model description

Folded LC geometries are modeled using a Gaussian function for the eclipses, and a cosine function with a period equal to half the orbital period for ellipsoidal-like variability¹, if present.

The eclipses are modeled with Gaussian functions of the form (1)where index i equals 1 and 2 for the primary (deepest) and secondary (least deep) eclipses, respectively, μ_i, d_i and σ_i being the Gaussian parameters and ϕ the observation phase (i.e. observation time modulo orbital period). The ellipsoidal-like variability, on the other hand, is modeled as , where A_ell is the peak-to-peak amplitude of the ellipsoidal-like variability, and ϕ_0,ell indicates whether the cosine is centered on eclipse 1 (ϕ_0,ell = μ₁) or on eclipse 2 (ϕ_0,ell = μ₂). The two-Gaussian model then writes (C is a constant) (2)Equation (2) includes the mirrors of eclipses 1 and 2 at phases from −2 to +2 in order to take into account the contribution of the tails of the Gaussian functions from adjacent phases due to the periodicity of the eclipses. The model parameters are illustrated in Fig. 1 for three types of EBs.

By convention, we shift LC times such as to locate the primary eclipse at phase 0. We therefore always have (3)even though we may continue to explicitly write μ₁ for clarity in some expressions.

Eclipse durations w_i (durations expressed in phase) are taken equal to the widths of the Gaussian functions at a magnitude depth of 2% relative to Gaussian depth d_i, that is w_i = 5.6σ_i, with an upper limit of 0.4. This somewhat arbitrary limit is set in order to avoid unphysical large eclipse durations for wide Gaussians. We thus have (4)Eclipse depths are taken equal to the difference between the magnitude at the bottom of the eclipse and the brightest magnitude G_min of the model: (5)Finally, we note that the constant C in Eq. (2) equals G_min only for detached EBs that do not show ellipsoidal-like variability (illustrated in the top panel of Fig. 1). For EB LCs displaying ellipsoidal-like variability (middle panel of Fig. 1) or for contact binaries (bottom panel of Fig. 1), C ≠ G_min.

2.2. Model computation

We fix the orbital period of each EB to the value published in the OGLE-III catalog.

A two-Gaussian model G(ϕ) defined by Eq. (2) is fit to the FLC { y_j(ϕ_j) } of each EB, where j is an index running over all measurements from 1 to the number N_obs of observations. The computation of the model parameters follows three steps: time series outlier removal (Sect. 2.2.1), initial values estimation of the two-Gaussian model parameters (Sect. 2.2.2), and non-linear fitting (Sect. 2.2.3).

2.2.1. Light curve outliers removal

Outlier removal is performed in two steps. First, all measurements with uncertainties greater than 1 mag are removed. Second, isolated measurements having magnitudes at the extremes of the magnitude distribution are removed. To do this, measurements with extreme magnitudes are identified from their deviations from the median magnitude when these exceed a certain number of times the inter-quantile range IQR (ten times at the faint side and two times at the bright side). They are considered to be outliers, and removed from the time series, unless they have similar (magnitude within 30%) neighbors in time (preceding or following measurement in time within a quarter of a day) or in the magnitude distribution (nearest points in the histogram of magnitudes).

2.2.2. Initial value determination of model parameters

Fitting a two-Gaussian model to a time series is very sensitive to the adopted initial values of the parameters. The better the initialization is, the better the convergence of the non-linear fitting algorithm is expected. We therefore proceed in three steps, first to catch the global shape of the FLC, then to detect the two eclipse candidates, and finally to initialize the two-Gaussian model.

Folded light curve smoothing.

We start performing a weighted running average on the FLC, replacing each magnitude value y_j at a given phase ϕ_j by a weighted average of the magnitudes y_k within a [ϕ_j−δϕ,ϕ_j + δϕ] phase window. The weights w_k are taken equal to (6)and the average magnitude is given by (7)with the index k running on all measurements available in the phase window. We take δϕ = 0.01. From this FLC {}, an evenly sampled FLC of 200 points is produced by linear interpolation of the averaged FLC. A smoothed FLC is then computed using the Savitzky-Golay (SG) algorithm (Savitzky & Golay 1964; Gorry 1990; Protopapas et al. 2006), which has the main advantage to preserve quite well the minima and widths of the eclipses. We use the Java implementation of the SG algorithm in the Flanagan library², which consists of a least-squares polynomial regression of degree 3 applied on 2M + 1 points centered on each considered point, M being a parameter which we take equal to 15. The resulting smooth FLC is denoted the SG FLC.

Eclipse identification.

Eclipse candidates are searched for in the SG FLC. We define a threshold magnitude equal to the median magnitude plus the median of the observation uncertainties, and determine a baseline magnitude M_b equal to the median magnitude of all observations brighter than this threshold. We then select the two faintest dips having magnitudes above this baseline as the two eclipse candidates.

Initial value estimation of model parameters.

The initial value of the constant C is set to the baseline magnitude M_b computed in the preceding step. The initial value of μ₁ (μ₂) is set equal to the phase of the measurement closest to the maximum magnitude of the deepest (second deepest) eclipse candidate identified in the SG FLC, while d₁ (d₂) is set to the difference between that maximum magnitude and the baseline magnitude. Finally, σ₁ (σ₂) is taken equal to 0.2 times the phase extent covered by all adjacent measurements around μ₁ (μ₂) fainter than the baseline magnitude.

2.2.3. Non-linear fitting procedure

We use the non-linear fitting algorithm nls of the R Project for Statistical Computing to search the solution to Eq. (2). We estimate the initial values of the parameters as explained in the previous section, swapping the order of the eclipses if necessary to have the first eclipse to be the deepest. A weight of is assigned to each measurement y_i, where ε_i is the uncertainty on y_i. We constrain the solutions to positive Gaussian depths and cosine amplitudes (the later constrain to avoid the non-linear method to converge to a sine solution) by transforming Eq. (2) such that it takes the logarithm of d₁, d₂ and A_ell. If, after convergence, the second eclipse turns out to be deeper than the first eclipse, we swap the two eclipses and search again for a solution because the cosine variability, if present, may impact the solution when the two eclipse candidates are not separated by exactly 0.5 in phase. This procedure ensures a consistent solution with the primary eclipse always numbered 1.

2.3. Model selection

Several models are tested on each LC, and their Bayesian information criterion (BIC) compared to identify the model that best matches the data given the measurement uncertainties. The BIC is computed as (Feigelson & Babu 2012, Eq. (3.54)) (8)where p is the number of model parameters, given in Table 1 for the models considered in this paper, and lnL is the log-likelihood given by (9)with ε_j being the uncertainty on measurement y_j.

Including two Gaussian functions and a cosine in Eq. (2) to model FLCs may lead to an overfit of the data if one or more of the components are insignificant (e.g., if the amplitude of the given component is small compared to the mean uncertainty of the measurements) or spurious (e.g., if the locations of the eclipse candidates were wrongly initialized). We therefore fit several models to each FLC, each having a different combination of components, and retain the one that has the highest BIC. The BIC takes into account the number of degrees of freedom of each model. We therefore avoid overfitting the data with non-significant components (either Gaussians or ellipsoidal variability).

The various models are summarized in Table 1. They comprise:

• a pure constant model (model “C”), representing a LC with nodetectable eclipse or ellipsoidal-like variability;

• a model including only an ellipsoidal component (model “CE”);

• models including only one Gaussian, without (model “CG”) or with (model “CGE”) an ellipsoidal component. For the CGE model, the ϕ_0,ell parameter in Eq. (2) is taken equal to the μ value of the eclipse candidate;

• models including two Gaussians, without (model “CG12”) or with (models “CG12E1” and “CG12E2”) an ellipsoidal component. Models “CG12E1” and “CG12E2” distinguish cases where the cosine of the ellipsoidal component is centered on the first (ϕ_0,ell = μ₁) or second (ϕ_0,ell = μ₂) eclipse candidate, respectively. They differ from one another only for eccentric systems for which μ₂−μ₁ ≠ 0.5.

The initial values of the model parameters for each of these models are taken from the set of initial values computed in Sect. 2.2.2. When testing models with only one Gaussian, if the procedure described in Sect. 2.2.2 identifies two eclipse candidates, two sets {CG, CGE} of models are tested, one set {CG1, CG1E1} with the Gaussian (and cosine when present) centered on the first eclipse candidate and another set {CG2, CG2E2} with the Gaussian (and cosine when present) centered on the second eclipse candidate. These two sets are not distinguished in Table 1, where models CG1 and CG2 (or CG1E1 and CG2E2) are indistinguishably referred to as model CG (or CGE).

Having computed, for a given FLC, all the models mentioned above, their BIC values are compared with one another, and the model with the highest BIC is retained. One exception to this rule concerns models with ellipsoidal-like variability that contain a wide Gaussian. These models are retained for comparison with the other models only if the phase duration of the eclipse(s) is(are) shorter than a given limit w_max,ell after model convergence, that is if σ < w_max,ell/ 5.6 for models with one Gaussian, or if both σ₁ < w_max,ell/ 5.6 and σ₂ < w_max,ell/ 5.6 for models having two Gaussians (we take w_max,ell = 0.4). Otherwise, the model is automatically rejected in favor of the models without ellipsoidal component. This condition is imposed in order to lift the degeneracy of some EB-type binaries for which the FLC can be modeled, for example, by either a CG12 model (correct model) or by a CGE model that includes a wide Gaussian on top of an ellipsoidal-like variability (fake model where the wide Gaussian added to one of the two depths of the ellipsoidal-like variability mimics a two-eclipse EB-type binary with non-equal depths). Tests performed on some EB-type binaries showed that the BIC value of the CGE model could indeed be larger than that of the CG12 value, thereby leading to a wrong model selection. Similarly to this example for models with one Gaussian, models with two Gaussians and an ellipsoidal component (models CG12E1 and CG12E2) are a-posteriori discarded if the phase duration of either Gaussian is larger than w_max,ell.

3. Application to OGLE-III LMC eclipsing binaries

We apply in this section the two-Gaussian model to the I-band LCs of the 26121 EBs in the LMC published by the OGLE-III survey (Graczyk et al. 2011). The survey operated from July 2001 to May 2009. Each LC has, in the mean, 500 measurements, 90% of which are in the photometric I band³.

We fix the orbital periods to the values listed by Graczyk et al. (2011). The LCs are then fit by a two-Gaussian model following the procedure described in Sect. 2, selecting the best model among a combination of Gaussian and cosine functions based on the BIC analysis (see Sect. 2.3). The computation takes less than 1 sec per source on a single 2.7 GHz CPU.

The vast majority (85%) of EB LCs are modeled with two Gaussians (with or without an additional ellipsoidal component), and 14% of LCs are modeled with only one Gaussian (with or without an ellipsoidal component). The remaining 1% EBs have their LCs modeled with only an ellipsoidal component, except for two cases for which the highest BIC model is a pure constant.

The success of our procedure to model the FLCs of the majority of OGLE-III EBs with two Gaussians does not guarantee reliable model components. The efficiency to identify true components in the LCs can be hampered by several effects, including time sampling, measurement uncertainties, wrong initial guess of eclipse locations, or additional intrinsic variability in one or both stars of the binary system. Some Gaussian or ellipsoidal-like components in the final models may therefore be missing, spurious, or wrong. We therefore devote this section to analyze model results.

We first present in Sect. 3.1 the phase coverage properties of the OGLE-III LMC EBs. Phase coverage depends only on observation time sampling and orbital period, but can impact model results. The significance and reliability of model components are then studied in Sect. 3.2. The question of model degeneracy is addressed in Sect. 3.3, and the quality of models is analyzed in Sect. 3.4. Finally, Sect. 3.5 mentions the table published at the CDS that provides the results of this paper for the OGLE-III LMC EBs.

3.1. Phase coverage

A good phase coverage of the eclipses is essential to correctly model the LCs. This, however, depends on the observation time sampling, on eclipse duration and on orbital period.

For a ground-based survey like OGLE, a significant phase gap is expected for orbital periods close to a multiple of the day (since the star is observable only during nights) or of the year (since the star is visible only at specific times of the year). The phases are then clumped in groups. Figure 2 plots the largest phase gap recorded in each OGLE-III LMC EBs FLCs versus their orbital period, highlighting in black crosses the FLCs that have phase-clumped data. The degree of clumpiness is evaluated based on the distribution of phase intervals between two successive measurements in the FLC. We do this after having shifted the phases of the FLCs by a constant value such as to move the largest phase interval to the end of the [0, 1] phase interval (this is done in order to have a quantity that is independent of the reference time used to fold the time series). We then define the phase clumpiness as the fraction of phase intervals that have durations less than 1/(N_obs−1). The phase clumpiness computed in this way is expected to be around 0.5 for sources regularly sampled in phase, and close to 1 for highly-clumped distributions in phase. Sources with a phase clumpiness above 0.75 are highlighted with black crosses in Fig. 2. Their two-Gaussian model parameters may be at fault due to missing data.

Another useful phase coverage related quantity is eclipse phase coverage by observations. We estimate the eclipse coverage by considering 11 equal phase intervals within the eclipse width [μ_i−w_i,μ_i + w_i] and by computing the fraction of these intervals that have at least one measurement. A value of 0 would mean that no observation is available within the considered interval – this can happen if observations are only available at the very borders of the eclipse candidate –, while a value of 1 means that measurements are available in all eleven phase bins. Eclipse candidates with insufficient eclipse coverage may have wrong model parameters. Such eclipse candidates are usually, but not always, spurious. About 92% of the sources in our sample have their eclipse candidates covered by observations over more than 70% of their durations. Eclipses with less than a few percent coverage are usually narrow in phase, irrespective of their period. This is shown in the eclipse width versus orbital period diagram displayed in Fig. 3, where eclipses with a phase coverage of less than 50% are highlighted in color.

3.2. Model components significance

The reliability of the Gaussian and ellipsoidal components found by the two-Gaussian model procedure is analyzed in this section. The analysis is done in Sect. 3.2.1 based on the BIC values obtained for different models. The significance of the eclipses and of the ellipsoidal component are then considered in Sects. 3.2.2 and 3.2.3, respectively.

3.2.1. Reliability of two-Gaussian model components

The reliability of a model component detected by the two-Gaussian procedure can be estimated by comparing the BIC value of the given model to the BIC value of the alternative model without the given component. The difference between these two BIC values, BIC(with component) − BIC(without component), is denoted by Δ_componentBIC. The distribution of Δ_Ecl1BIC for the primary eclipse candidate, of Δ_Ecl2BIC for the secondary eclipse candidate, and of Δ_EllBIC for the ellipsoidal component are shown in Fig. 4 for the various models (from models containing two Gaussians and an ellipsoidal component in the top panel to models containing only an ellipsoidal component in the bottom panel). The larger the Δ_componentBIC difference is, the larger the probability is for the given component to be significant and non-spurious. In the great majority of cases, the reliability of the eclipse candidates is very good, with 91% (66%) of models with two Gaussians satisfying Δ_ecl1BIC > 100 (Δ_ecl2BIC > 100) for the primary (secondary) eclipse irrespective of the presence of an ellipsoidal component, and 81% of models with one Gaussian satisfying Δ_eclBIC > 100 irrespective of the presence of an ellipsoidal component.

The significances of the various components in models containing two Gaussians are shown in Fig. 5, where Δ_ecl2BIC is plotted versus Δ_ecl1BIC with Δ_EllBIC shown in color for models with an ellipsoidal component. Two examples with highly significant eclipses are shown in the top panels of Fig. 6. Small values of Δ_ecl1BIC or Δ_ecl2BIC, on the other hand, point to unreliable primary or secondary eclipse candidates, respectively. Fortunately, this concerns only a small fraction of the eclipse candidates, as seen in Fig. 5. This feature must however be kept in mind when studying the ensemble of EBs with the two-Gaussian model results.

A similar analysis can be done on models containing one Gaussian. The distribution of Δ_eclBIC and Δ_ellBIC for the ones containing an ellipsoidal component is displayed in Fig. 7, and two examples with highly significant components are shown in the two bottom panels of Fig. 6.

3.2.2. Eclipse significance

It is instructive to further analyze the characteristics of the least reliable (according to Δ_eclBIC) eclipse candidates. One expects these candidates to have Gaussian depth d_i comparable to, or smaller than, the measurement uncertainties. This is confirmed in Figs. 8 and 9, which show the ratio of the Gaussian depth over mean measurement uncertainty inside the eclipse versus Gaussian width σ_i for primary (i = 1) and secondary (i = 2) eclipse candidates, respectively. Sources that have Δ_ecliBIC < 100 are shown in color in Figs. 8 and 9; they are seen to have the lowest ratios.

Figures 8 and 9 further show a dependency of eclipse reliability on eclipse width. Narrow eclipses require larger ratios than wide eclipses do in order to be significant, because narrow eclipses contain, on the mean, less measurements than wide eclipses. Therefore, the narrower the eclipse is, the deeper it must be to be reliably detected.

A (small) fraction of eclipse candidates have Gaussian depth to mean measurement uncertainty ratios that are off the bulk distribution in the versus σ_i diagrams shown in Figs. 8 and 9. They concern very narrow eclipse candidates (σ ≲ 10^-3) with ratios that can reach above 100. They are mainly eclipse candidates that lack sufficient observations inside the eclipse. As a result, the Gaussian depth cannot be well constrained, and an unrealistically deep Gaussian is adopted by the model fitting algorithm with a concomitant large uncertainty. This is verified in Fig. 10, which shows the relative uncertainty d_1,err/d₁ of the Gaussian depth versus eclipse significance Δ_ecl1BIC of the primary eclipse candidates of all models containing two Gaussians. The eclipse coverage factor, shown in color for eclipses that have a coverage less than 50%, is seen to be small for the models with small eclipse significance and/or with large relative uncertainty of the Gaussian depth. They usually correspond to spurious eclipse candidates.

3.2.3. Ellipsoidal component significance

The amplitude A_ell of the ellipsoidal component can be relatively large compared to the Gaussian depth of the primary eclipse candidate. The ratio A_ell/d₁ is greater than 0.1 for 69% (82%) of models containing two (one) Gaussians, and greater than 0.5 for still 15% (29%) of cases. The histograms of the distributions of this ratio are shown in Fig. 11.

The cosine function used in our models describes any type of ellipsoidal-like variability that would be present in the LCs. For circular systems, the two eclipses are separated from each other by 0.5 in phase, and the ellipsoidal component, if present, is centered on both Gaussians (i.e. cos4πμ₁ = 1 and cos4πμ₂ = 1). The case of elliptical systems containing two eclipses needs additional investigation. For these systems, an ellipsoidal-like variability added to the two-Gaussian model would have the cosine centered on one of the two Gaussians and displaced relative to the other Gaussian. Figure 12 plots the amplitude of the ellipsoidal component versus phase separation | μ₂−μ₁ | between the two Gaussians of all CG12E1 and CG12E2 models. It shows that the majority of models containing an ellipsoidal component have either a near-circular orbit (| μ₂−μ₁ | close to 0.5) or a small ellipsoidal component (A_ell < 0.05 mag). Investigation of the few non-circular model candidates with a significant ellipsoidal component show that one of their Gaussians may be a spurious candidate (Δ_ecl1BIC or Δ_ecl2BIC < 50, shown with crosses in Fig. 12).

3.3. Model degeneracies

The automated selection of a best model among several ones inevitably raises the question of model degeneracies given the data uncertainties. In the case of two-Gaussian models, model degeneracy can arise because, for example, the cosine function describing the ellipsoidal-like variability can be mimicked by two wide Gaussian functions mirrored over adjacent phases. Likewise, EB- and EW-type EBs can mathematically be modeled, within the measurement uncertainties, by either two wide Gaussians or an ellipsoidal component complemented by a wide Gaussian to account for the different eclipse depths. Degenerate models are expected to have their BIC values close to one another, all of them describing almost equally well the data within the given measurement uncertainties. We therefore estimate the degree of degeneracy between two models A and B by the absolute difference | BIC_A−BIC_B | between their BIC values.

Figure 13 illustrates model degeneracy for EBs that have the CGE model (one Gaussian + ellipsoidal component) selected by the automated procedure. A CGE model can be degenerate with either a CG12 model (two Gaussians) or a CE model (only an ellipsoidal component). The figure plots the degree of degeneracy of the CGE model with a CG12 model on the Y-axis versus the degree of degeneracy with a CE model on the X-axis. Models located on the left part of the diagram (small BIC_CGE−BIC_CE) may be confused with a purely ellipsoidal model, while those in the lower part (small BIC_CGE−BIC_CG12) may be confused with a model containing only two Gaussians. These degeneracies should be taken into consideration when performing statistical studies on an ensemble of EBs. This exercise may be more or less straightforward depending on the type of degeneracy. A degeneracy between CGE and CG models, for example, is not very harmful because CG models are equivalent to CGE models with zero amplitude of the ellipsoidal component. A degeneracy between CGE and CG12 models, on the other hand, is more problematic. In this case the alternative CG12 model would be composed of a narrow deep eclipse and a wide shallow eclipse, the astrophysical origin of which would be more challenging to find.

3.4. Quality of the models

We analyzed in the previous sections the significance of model components. We now want to check how suitably the two-Gaussian models describe the variability patterns present in the FLCs of EBs. To achieve this, we use the Abbe value that quantifies the degree of smooth variability present in a curve (see Mowlavi 2014, and references therein). Given a series of n values y_{j = 1 → n}, the Abbe value is defined by (10)where is the mean of { y_j }.

The Abbe values of the original and residual FLCs are noted and , respectively. A FLC with no visible variability pattern will have an value around 1 (there is no correlation between successive y_{j + 1}−y_j differences), while a very clear and smooth variability pattern will result in a value decreasing to 0 (y_{j + 1}−y_j differences are small compared to the standard deviation of the series). If a model successfully describes a FLC, no variability pattern should subsist in the residual FLC and the Abbe value of the residual LC should be close to 1.

The histograms of and are shown in Fig. 14. The Abbe values of the original FLCs are seen to have an almost flat distribution between 0 and 0.6, and to start to decrease above 0.6. This is expected, since an EB with is more difficult to be identified, and hence has a smaller probability to be in the OGLE-III catalog of EBs in the first place. The Abbe values of the residual FLCs after model subtraction, on the other hand, peaks at 0.95 (thick red histogram in Fig. 14). This reflects the efficiency of the two-Gaussian models to adequately fit the geometry of EB FLCs, thereby increasing the Abbe value from values below 0.8 in the original FLC to values above 0.8 in the residual FLC. It does not guarantee, though, an adequate identification of eclipse and/or inter-eclipse components for the EB, which must be studied using component significances as done in the previous sections. But it reveals that no significant variability pattern remains in the FLC after model subtraction.

A small fraction of EBs have a residual Abbe value below 0.8, as shown by the tail distribution of in Fig. 14. In these cases, a variability pattern that can be significant remains in the FLC after model subtraction. To further analyze these cases, the Abbe value is plotted versus in Fig. 15 for all models containing two Gaussians (irrespective of the presence of an ellipsoidal component). Three regions are identified in the diagram, named region A, B and C.

Region A is defined as . It contains sources with FLC geometries that are well described by the two-Gaussian model within the measurement uncertainties. This represents by far the majority of cases, with 94% of all OGLE-III EB LCs falling in this region. The FLCs of several examples labeled in Fig. 15 are shown in Fig. 16 with, from top to bottom panel, increasing (i.e. decreasing LC signal-to-noise ratio).

Region B is defined as , . In this region of the diagram, the signal-to-noise of the LCs is very high (with a resulting ), and the well defined FLC geometry challenges the two-Gaussian model (as seen from the variability pattern still present in the residual FLC, with ). Only 2.4% of sources modeled with two Gaussians fall in this region of the versus diagram. The two-Gaussian model can fail to adequately describe the FLC geometry for two reasons. First, the LC geometries during the eclipse and inter eclipse phases are more complex than what can be described by simple Gaussian and cosine functions, respectively. The three top FLCs in Fig. 17 illustrate such cases, with the residual FLCs displayed in a panel below each FLC. Nevertheless, the examples show that the two-Gaussian models still successful grasp the main properties of the eclipses despite the simplicity of the models.

The two-Gaussian model can also fail to adequately describe the geometry of an EB FLC if a physical effect other than an eclipse or ellipsoidal-like variability is present in the LC. This is the case, for example, if the system has a total eclipse or a reflection component. An example of each of these two cases occurring in region B is shown in the bottom panels of Fig. 17.

Finally, region C gathers the remaining part of the versus diagram (i.e. , ). Sources in this region should be successfully modeled by a two-Gaussian model, because the relatively low S/N is less demanding of the model than in region B. Failure of the two-Gaussian model to do so would imply either additional physics not accounted for in the two-Gaussian model, or more fundamental issues to be investigated. This region contains thus potentially interesting cases of outliers to be investigated. They will be addressed in Sect. 4.1.1.

Sources falling in region C of the versus diagram are expected to have a large dispersion in their residual LCs. Therefore, their reduced χ², defined by (11)where p is the number of parameters in the model and ε_i is the uncertainty on the magnitude of measurement i, should be large. Figure 18 plots versus . The majority of sources are seen in the figure to have and . A stream of sources is also seen toward lower values with a correlated increase of . This is expected, since values of smaller than 0.7 indicate the presence of residual variability patterns that result in larger values. And indeed, the sources labeled in region C of Fig. 15 lie on or close to this stream of data points in Fig. 18.

Figure 18, however, shows the existence of a subset of sources that have larger than what is expected from the bulk or stream distributions of points in the figure. They indicate the presence of an additional variability component of a different nature, that breaks the smoothly varying pattern of a FLC derived from a strictly periodic variability. These cases will further be studied in Sect. 4.1.2.

3.5. Table summary

All quantities derived in this study for the OGLE-III LMC EBs are published in a table at the CDS. A description of the table content is given in Appendix A.

4. Discussion

We present two application examples of the two-Gaussian models. The first one (Sect. 4.1) aims to identify binary systems in physical configurations incompatible with two-Gaussian models. We refer to these systems as outliers. The second example (Sect. 4.2) shows how the two-Gaussian model results can be used to study statistical properties of the ensemble of EBs. They are given here for illustrative purposes only, a full study of each of these two applications is beyond the scope of this paper.

4.1. Identification of outlying cases

The choice of Gaussian and cosine functions to model the FLC geometry of eclipse and ellipsoidal variability, respectively, defines the set of EB configurations than can be described by the two-Gaussian models. Any deviation from this set of configurations will be detectable through poor model fit quality. We use here the two diagnostic tools presented in Sect. 3.4 to evaluate model fit quality: the versus diagram (Fig. 15) and versus diagram (Fig. 18). We discuss these two diagrams in Sects. 4.1.1 and 4.1.2, respectively.

4.1.1. Outliers from the versus diagram

We select all systems from region C of the versus diagram (Fig. 15). A zoomed version of the figure is shown in Fig. 19. The FLCs of all sources that lie within the lower-right area delimited by the dashed green line in Fig. 19 have been visually inspected and classified in one of the categories described below. Examples are provided in Figs. 20 and 21.

Total eclipse

(blue filled squares in Fig. 19). The presence of a total eclipse manifests itself by a flat bottom in the LC during the eclipse. This is poorly approximated by a Gaussian function, and will result in an value below 0.8. An example is given with source 20061 in Fig. 20. Systems observed with close-to-total eclipses will also have LC geometry that deviates from Gaussian, because of the steep ingress and egress curves. A limiting case is given by a system containing two stars of equal radii in circular orbit, for which the LC during the eclipse will have a V-shaped geometry if the system is seen edge-on.

Semi-detached morphology

(gray filled diamonds in Fig. 19). Systems that have one of the stars filling or close to filling its Roche lobe will display a non-cosine LC shape between the eclipses. An example is given with source 13836 shown in Fig. 20.

Reflection effect

(black filled diamonds in Fig. 19). The LCs of some systems show an out-of-eclipse brightening around the secondary eclipse. An example is given with sources 10156 and 9002 in Fig. 20. This can be due to reflection, where the hotter star heats the surface of the cooler star that faces the hot star. For source 9002, a lag is visible between the phase of the secondary eclipse and the phase at maximum luminosity, which could be caused by stellar rotation. We classify these systems as having a reflection signature in their LC. Their LCs cannot be modeled with a cosine function with half the orbital period used to model ellipsoidal-like variability, but could successfully be described with a cosine function with a period equal to the orbital period (Moe & Di Stefano 2015).

Tidal distortions in eccentric binaries and heartbeat stars

(red filled triangles in Fig. 19). The effect on the LC of tidal distortions in eccentric binaries also appears in region C of the versus diagram. The FLCs of four cases showing LC deformations due to such effects are given in Fig. 20 with sources 5347, 3512, 12035 and 23999. Various LC geometries due to tidal distortions have been reported by Thompson et al. (2012) in the Kepler data.

Large intrinsic scatter

(magenta open circles in Fig. 19). Some systems display a scatter in their residual LC larger than what is expected from the measurement uncertainties. Two such cases are shown in Fig. 21 with sources 16745 and 17782. They are further discussed in Sect. 4.1.2.

Failed convergence

(gray cross and plus signs in Fig. 19). The mismatch results from either a failure to correctly identify the initial locations of the eclipses or to converge on the two-Gaussian model (gray crosses in Fig. 19), or from a wrong initial orbital period (gray plus sign in Fig. 19). Only one clear case of the last category is found in the OGLE-III catalog of LMC EBs, for which the double of the true period is reported in the OGLE-III catalog.

4.1.2. Outliers in the versus diagram

The versus diagram introduced in Sect. 3.4 offers a second interesting tool to identify outlying cases of EB LCs. A value larger than 0.8 for indicates a reasonable fit of the geometry of the initial FLC by the two-Gaussian model. The resulting reduced χ² should then be small. For a fraction of the LCs, however, is still large, as has been seen in Fig. 18. A residual scatter is thus present with an amplitude larger than expected from the measurement uncertainties. Here, we check the status of these stars through a visual check of their LCs.

We select all stars that have and , and identify example cases illustrating various potential origins for the higher-than-expected scatter in the residual FLC. The sources chosen as examples are highlighted in the versus diagram shown in Fig. 22, and their LCs are shown in Figs. 23 and 25. The following cases are identified.

Intrinsic periodic variability.

If one or both stars are intrinsically variable, a residual scatter is naturally expected in the FLC. This may be the case for source 16745 shown in Fig. 21, which has a large (see Fig. 22). Two other cases, sources 22288 and 2547, are shown in the top panels of Fig. 23, with intrinsic variability time scales long enough to be visible in their LCs.

Intrinsic non-periodic variability.

An aperiodic variability can originate, for example, from flares, outbursts, or irregular variability. We visually identified several EBs presenting flares in the selected region of the versus diagram. An example is shown in Fig. 23 with source 7551. It is characterized by bright outlying measurements apparently randomly distributed in the FLC. The LC reveals three flares with time scales of the order of 100 days. The source is blue (V−I = −0.206 mag) and has .

Two examples showing an outburst are given in Fig. 23 with sources 18577 and 2594. The LC shapes of these particular cases resemble that of microlensing events. However, being blue (V−I = −0.153 and −0.047 mag, respectively), they may be blue bumpers (Cook et al. 1995; Wyrzykowski et al. 2011), as also suggested by Graczyk et al. (2011) for source 2594.

Finally, an example of a source with irregular intrinsic variability is shown in Fig. 23 with source 9699. It shows irregular brightening and fading, on time scales of tens of days for the brightenings and hundred of days for the fadings. The EB most probably hosts a Be star with a moderately blue color of V−I = 0.15. It has

Apsidal motions.

Apsidal motion systems result from the rotation of the line of apsides, which is a consequence of non-axial distribution of component mass, leading to torque exerted on the Runge-Lenz vector. It can be effectively modeled by a linear rate of change of the argument periastron, which manifests itself as an eclipse timing variation that causes both eclipses to excurse with respect to one another from their initial position. Thus, both eclipses witness phase shifts, leading to a measurement of an anomalous orbital period when the phase of the system is defined with respect to superior conjunction. We found almost twenty sources showing in-eclipse scatter of the measurements with an amplitude larger than expected from the out-of-eclipse scatter. Three examples are shown in Fig. 24 with sources 10157, 19624 and 19879. For source 10157, the appearance and disappearance of the eclipses is probably due to precession (Graczyk, priv. comm.). Such systems have previously been identified by Graczyk et al. (2011) as transient eclipsing binaries.

Multiple systems.

Multiple systems can reveal themselves through the presence of several periods in the LC for specific orbit configurations with respect to the line of sight. Source 16549 shown in Fig. 25 is an example of a hierarchical, gravitationally bound system that imprints its signature in the LC. The four-body system is composed of two EB components, one with a long period reported in the OGLE-III catalog to be of 164.79 d, and a second one with a short period of 0.818033 d. The LC of the system folded on the long period, shown in Fig. 24, clearly shows a narrow eclipse caused by the long-period binary component. The period could actually be double this value, which would then reveal the presence of two eclipses in the FLC. An analysis of the residual LC performed by (Graczyk et al. 2011, see in particular their Fig. 11) reveals the presence of the additional short-period, EB-type, contact system. The contact binary introduces a scatter of ~0.15 mag in the residual LC of the long-period system (which has a primary depth of 0.68 mag), that translates to using our two-Gaussian model for the long-period system.

Disks.

The presence of disks around one or both stars in a binary system can reveal itself in the LC geometry in and/or around the eclipses. A nice example is given by source 17782 displayed in Fig. 21. The source has been discussed by Graczyk et al. (2011) who conclude on the presence of a disk that contributes to a wide plateau in the primary eclipse superposed on a narrower stellar eclipse, with a morphology of the disk-induced eclipse that changes with time (see their Fig. 13). The source is easily identified as an outlier in the versus diagram of Fig. 22, with .

Misclassification.

A large value can also result from a misclassification of the source. Sources 957 and 7435, for example, shown in Fig. 23, have variable light variation amplitudes with time and may be long period variables (LPVs) instead of EBs. This would be consistent with the red color of the two sources) (V−I = 1.87 and 2.21 mag, respectively) and their long periods of variability. We note that the period of source 957 would then be ~207 d if it was a LPV, that is half the quoted value of ~414 d in the OGLE-III catalog. In addition, a variability on time scales of several dozens of days is visible in the FLC and LC of this source.

Potential data reduction issues.

Finally, the versus diagram can also serve as a diagnostic tool for data reduction quality. Problems in data reduction will lead to artificially increased scatter in the residual LC. The source will then appear as an outlier in the diagram, like the other sources analyzed above. An example is shown in Fig. 23 with source 11591, which displays a doubling of the LC toward the end of the OGLE-III survey. Few such cases have been identified from our visual inspection of the selected region of outliers in the versus diagram.

4.2. Statistical analysis

We illustrate in this section the usage of the two-Gaussian models by analyzing the projected orbital eccentricities and eclipse widths of all models containing two significant eclipse candidates. We filter the initial data set of all OGLE-III EBs of the LMC in several steps. We first select all EBs for which the two-Gaussian model successfully describes the geometry of the FLC (i.e. , see Sect. 4.1.1) and which have a scatter in the residual LC (see Sect. 4.1.2) smaller than . This set contains 92% of the initial OGLE-III catalog of LMC EBs. We then take all EBs that are modeled with two Gaussians (i.e. for which the CG12, CG12E1 or CG12E2 has the largest BIC). This represents 85% of the previous set. Finally, we restrict to models having significant eclipse candidates. We use the significance criterion based on the Δ_Ecl1BIC and Δ_Ecl2BIC quantities introduced in Sect. 3.2. For the illustrative purposes of this section, we retain only models which have Δ_Ecl1BIC > 50 and Δ_Ecl2BIC > 50 (see histograms of these quantities in Fig. 4). This represents 77% of the previous set of models containing two Gaussians. In total, our final set of EBs containing two significant eclipse candidates contains 15681 sources.

We take the deviation | μ₂−μ₁−0.5 | of the eclipse separation in phase with respect to 0.5 as a proxy for the projected eccentricity. This quantity is plotted versus the orbital period in the top panel of Fig. 26. The circularization of the orbit as the period shrinks is well visible in the figure. A zoom at short periods is shown in the bottom panel of that figure. The number of eccentric binaries decreases drastically at periods below 2 days, and all binary systems are found to be circular for periods shorter than ~1.2 d.

Figure 27 shows the widths w₁ and w₂ of the eclipses (see Eq. (4)). As expected, the majority of EBs have primary and secondary eclipse widths that are within a factor of two of each other (i.e. within the area delimited by the dashed lines in Fig. 27) . Some examples that deviate from this rule are shown in Fig. 28. The two cases in the top two panels have very narrow eclipses, and are on an elliptic orbit. The third case from top has an out-of-primary-eclipse geometry that is modeled by a wide second Gaussian, possibly representing an eccentric system with reflection. The last case in the bottom panel of Fig. 28 shows an ellipsoidal variability with large amplitude, possibly representing an eccentric system with ellipsoidal variability. The LC of source 10282 may also result from a system with a disk and a strong reflection effect (Garczyk, priv. comm.).

The difference of eclipse width distributions between eccentric and circular systems is illustrated in Fig. 29. The top panel of the figure shows systems with | μ₂−μ₁−0.5 | < 0.002, which favors⁴ systems in circular orbits. The bottom panel shows all systems with | μ₂−μ₁−0.5 | > 0.05, which selects eccentric systems.

5. Conclusions

This work has shown the potential of Gaussian and cosine functions to model the geometry of EB LCs resulting from eclipse and ellipsoidal-like variability, respectively. Using the two-Gaussian model parameters, we were successful in achieving our two goals, that is to identify outliers in a large set of EBs, and to provide a data base for the study of EB parameters on a statistical ground.

Key to these achievements are two diagrams introduced in Sect. 3.4. The first is the versus diagram (Fig. 15) that enables to identify outliers in terms of deviation of the FLC geometry from what can be modeled with a combination of Gaussian and cosine functions. The second diagram is the versus diagram (Fig. 18) to identify EBs that contain additional intrinsic variability other than that resulting from the binary nature of the source. These two diagrams have been exploited in Sect. 4.1 to identify potentially interesting binary systems.

Section 4.2 has then briefly illustrated how the two-Gaussian model results can be used to study the properties of EBs on a statistical ground. An inevitable challenge of automated procedures is to minimize as much as possible the contamination of statistical conclusions by the presence of non-physical components in the models. We presented in Sect. 3.2 a method based on BIC analysis to estimate the significance of each component in the two-Gaussian models. In particular, the significances of primary (Δ_Ecl1BIC) and secondary (Δ_Ecl2BIC) eclipse candidates can be used to filter out models that have a high probability to contain spurious eclipses.

The results of our two-Gaussian models for the OGLE-III EBs of the LMC are available at the CDS. A description of its content is given in Table A.1.

This work constitutes a basis for the establishment of an automated pipeline to process Gaia LCs. Gaia LCs will have, on the mean, about 70 measurements on a five-year mission. The efficiency of the two-Gaussian model to characterize the LCs of the EBs seen by Gaia has been addressed by Kochoska et al. (2017). In that study, the two-Gaussian method has been applied to both the original Kepler LCs and to the Kepler set of LCs resampled with Gaia cadence using the Gaia scanning law at the sky position of the Kepler EBs, and considering a five-year time span for the Gaia mission. The study reveals that 2/3 of the Kepler EBs are detectable by Gaia due to the presence of a sufficient number of observations in the eclipses. The study further shows that, when this is the case, the two-Gaussian method is successful in characterizing the LC geometry of the EBs. Kochoska et al. (2017) further propose a classification scheme of the detectable sources based on the morphological type indicative of the light curve.

Several improvements to the two-Gaussian model are foreseen for the next steps. They comprise the inclusion of an additional component in the models to describe reflection. We also pursue our exploratory works of automated classification techniques initiated with the works of Kochoska et al. (2017) and Süveges et al. (2017). On the path to these realizations, the various procedures will be applied and tested on existing data from surveys such as Hipparcos and the recently-released OGLE-IV, as well as on Gaia-simulated data.

---

References

Appendix A: Table description of the two-Gaussian model results

Table A.1 summarizes the content of the electronic table giving the two-Gaussian models for all the OGLE-III LMC EBs. In the electronic version, a “NA” is published when a quantity is not applicable for a given source, for example for the parameters of a secondary eclipse when the model contains only one Gaussian. The table contains

• the source ID number and the orbital period given in the OGLE-IIIcatalog for the LMC EBs;

• the two-Gaussian model chosen by our automated procedure based on the BIC (see Sect. 2.3);

• the epoch of primary eclipse minimum and the two-Gaussian parameters defined in Eqs. (1) and (2). The uncertainties associated to the parameters are taken from the covariance matrix returned by the non-linear fitting procedure;

• the widths and depths of the eclipses, defined by Eqs. (4) and (5), respectively;

• the maximum phase gap, the phase clumpiness, and the eclipse phase coverages defined in Sect. 3.1;

• the significances Δ_componentBIC of the two-Gaussian model components described in Sect. 3.2;

• the depths of the eclipse candidates relative to the mean measurement uncertainties inside the eclipses, discussed in Sect. 3.2.2;

• the Abbe and reduced χ² values of the FLCs, introduced in Sect. 3.4 to evaluate the overall quality of the fits;

• and the BIC values of all the two-Gaussian models evaluated for each EB. Models containing one Gaussian, named CG (CGE) in Table 1 when they do not (they do) contain an ellipsoidal-like variability, are divided in Table A.1 into CG1 (CG1E1) and CG2 (CG2E2) depending on whether the unique Gaussian is centered on the first or second eclipse candidate identified in the initial value determination step of model parameters (see Sect. 2.2.2). These two initial parameter sets are tested in succession when evaluating for the best model. The distinction between between CG1 (CG1E1) and CG2 (CG2E2) models is also necessary for the computation of the eclipse significances Δ_ecl1BIC and Δ_ecl2BIC of models containing two Gaussians. If the initial value determination procedure of model parameters finds only one eclipse candidate, models CG2 and CG2E2 are non-existent. The BIC values of some models may also be inexistent if the non-linear procedure fails to converge.

All Tables

All Figures