# Early-type Eclipsing Binaries with Intermediate Orbital Periods

###### Abstract

We analyze 221 eclipsing binaries (EBs) in the Large Magellanic Cloud with B-type main-sequence (MS) primaries ( 4 - 14 M) and orbital periods = 20 - 50 days that were photometrically monitored by the Optical Gravitational Lensing Experiment. We utilize our three-stage automated pipeline to (1) classify all 221 EBs, (2) fit physical models to the light curves of 130 detached well-defined EBs from which unique parameters can be determined, and (3) recover the intrinsic binary statistics by correcting for selection effects. We uncover two statistically significant trends with age. First, younger EBs tend to reside in dustier environments with larger photometric extinctions, an empirical relation that can be implemented when modeling stellar populations. Second, younger EBs generally have large eccentricities. This demonstrates that massive binaries at moderate orbital periods are born with a Maxwellian “thermal” orbital velocity distribution, which indicates they formed via dynamical interactions. In addition, the age-eccentricity anticorrelation provides a direct constraint for tidal evolution in highly eccentric binaries containing hot MS stars with radiative envelopes. The intrinsic fraction of B-type MS stars with stellar companions / 0.2 and orbital periods = 20 - 50 days is (7 2)%. We find early-type binaries at = 20 - 50 days are weighted significantly toward small mass ratios 0.2 - 0.3, which is different than the results from previous observations of closer binaries with 20 days. This indicates that early-type binaries at slightly wider orbital separations have experienced substantially less competitive accretion and coevolution during their formation in the circumbinary disk.

###### Subject headings:

binaries: eclipsing, close; stars: massive, formation, evolution, statistics## 1. Introduction

It has long been understood that the main-sequence (MS) binary star fraction increases with primary mass (Abt, 1983; Duquennoy & Mayor, 1991; Fischer & Marcy, 1992; Raghavan et al., 2010; Duchêne & Kraus, 2013, etc.). Indeed, most massive stars with 10 M will interact with a stellar companion before they explode as core-collapse supernovae (Sana et al., 2012). Throughout the decades, there have been significant advances in the detection of close and wide companions to massive stars (Wolff, 1978; Garmany et al., 1980; Levato et al., 1987; Abt et al., 1990; Shatsky & Tokovinin, 2002; Kouwenhoven et al., 2007; Sana et al., 2012; Rizzuto et al., 2013; Kobulnicky et al., 2014). However, the intrinsic properties of binary companions to early-type primaries, e.g. their eccentricity and mass-ratio distributions, remain elusive at intermediate orbital periods. The major goal of this work is to help fill this particular portion of the parameter space.

Eclipsing binaries (EBs) offer a key to the accurate measurement of the binary properties of early-type stars. Large photometric surveys, such as the third phase of the Optical Gravitational Lensing Experiment (OGLE-III), have discovered tens of thousands of EBs (Graczyk et al., 2011; Pietrukowicz et al., 2013). These populations of EBs are orders of magnitude larger than previous binary samples. Despite the geometrical selection effects, we can still achieve large sample statistics to reliably infer the intrinsic binary fraction and properties at intermediate orbital periods. We emphasize that EBs can probe a unique portion of the binary parameter space unavailable to other observational techniques.

In Moe & Di Stefano (2013, hereafter Paper I), we incorporated OGLE catalogs of EBs in the Large and Small Magellanic Clouds (LMC and SMC, respectively) as well as Hipparcos observations of EBs in the Milky Way. We compared the close binary properties ( 20 days) of early-B MS primaries in the three different galaxies. The Milky Way and SMC EB samples are too small to warrant an analysis of period-dependent binary properties. The OGLE-III LMC EB catalog (Graczyk et al., 2011), on the other hand, contains 5 - 40 times more systems, is relatively complete toward shallow eclipse depths, and includes the full I-band and V-band light curves.

In Moe & Di Stefano (2014, hereafter Paper II), we developed a three-stage automated pipeline to analyze EBs with short orbital periods in the OGLE-III LMC database. This pipeline (1) classifies EBs according to their light curve characteristics, (2) measures the intrinsic physical properties of detached EBs, e.g. ages and component masses, based on the observed radii and temperatures, and (3) recovers the intrinsic binary statistics by correcting for selection effects.

In the present study, we utilize EBs in the OGLE-III LMC database to measure the binary fraction, mass-ratio distribution, and eccentricity distribution of B-type MS stars with intermediate orbital periods 20 - 50 days. We organize the rest of this paper as follows. In §2, we define our selection criteria for identifying EBs with B-type MS primaries, intermediate orbital periods, and well-defined eclipse parameters. We next describe an automated procedure we developed to fit detailed physical models to the observed EB light curves, and we present our results for the physical properties of the individual EBs (§3). In §4, we explain the observed trends in the measured EB parameters, paying special attention to the empirical age-extinction and age-eccentricity anticorrelations. We then perform Monte Carlo simulations to quantify selection effects, and present our results for the corrected binary statistics (§5). We summarize our main results and conclusions in §6.

## 2. EB Selection and Classification (Stage I)

In Paper II, we developed a three-stage automated pipeline to fully analyze short-period EBs in the OGLE-III LMC database. In the present study, we adapt our routine to identify intermediate-period EBs with well-defined light curves (Stage I - this section), measure their physical properties (Stage II - §3), and correct for selection effects (Stage III - §5). EBs with intermediate orbital periods exhibit two major differences that must be considered. First, the eclipse widths and , which are expressed as a fraction of the orbital period , become narrower with increasing orbital separation. Given the average number 470 of I-band measurements in the OGLE-III LMC survey (Graczyk et al., 2011), the light curves are not sufficiently sampled if either of the eclipse widths 0.0021 are too narrow. EBs with small MS components and long orbital periods 50 days have narrow eclipses 0.002, and are therefore not Nyquist sampled. This subsampling leads to detection incompleteness, issues with aliasing, and the inability to fully characterize their intrinsic physical properties. Hence, it is the finite cadence of the OGLE-III observations, not geometrical selection effects, that limits our present study of EBs to = 20 - 50 days (see also Söderhjelm & Dischler 2005).

Second, the majority of early-type EBs at 20 days are in eccentric orbits. We must therefore adapt our physical models to simultaneously fit the eccentricity and argument of periastron (§3). In addition, it is possible for an eccentric binary to have a certain combination of eccentricity, periastron angle, and inclination that is sufficiently offset from edge-on (e.g., 86) so that there is only one eclipse per orbit. Indeed, there are many EBs with single eclipses in the OGLE-III LMC database (see below). Unfortunately, we cannot measure the physical properties of these systems. We therefore remove single-eclipse EBs from our well-defined sample, and we account for their removal when we correct for selection effects (see §5). In the following, we review our methods from Paper II, where we pay special attention to the nuances of EBs with intermediate orbital periods.

In this study, we select the 96,000 systems in the OGLE-III LMC catalog (Udalski et al., 2008) with mean magnitudes 16.0 17.6 and observed colors 0.25 0.20. Given the distance modulus = 18.5 to the LMC (Pietrzyński et al., 2013) and typical dust reddenings 0.1 - 0.3 mag toward hot young stars in the LMC (Zaritsky et al., 2004), these stars have luminosities and surface temperatures that correspond to B-type MS primaries. From this sample, we analyze the 221 systems that were identified as EBs with orbital periods = 20 - 50 days (Graczyk et al., 2011). In Table 1, we list the OGLE-III LMC EB identification numbers, observed colors , and numbers of I-band measurements for each of these 221 EBs.

As in Paper II, we measure the intrinsic rms scatter in the -band light curve outside of eclipses for each EB. We then calculate the correction factor 1.0, i.e. the ratio between the actual rms scatter and photometric uncertainties reported in the catalog. For each -band measurement in an EB light curve, we multiply the listed photometric uncertainties by the correction factor to determine the corrected uncertainties.

We classify EBs based on an analytic light curve model of two Gaussians with eight total free parameters. The orbital phase 0 1 is determined by the time of observation and two model parameters: the orbital period (in days) and epoch of primary eclipse minimum (Julian date 2450000). The six remaining analytic model parameters are the average I-band magnitude outside of eclipses , primary and secondary eclipse depths and , primary and secondary eclipse widths and , and the phase of secondary eclipse . The analytic model of Gaussians is:

(1) |

We fit this analytic model to each EB I-band light curve. Specifically, we utilize an automated Levenberg-Marquardt algorithm ( MPFIT, Markwardt, 2009) to minimize the statistic. The MPFIT routine provides robust best-fit solutions and measurement uncertainties for the eight analytic model parameters. Some of the photometric measurements are clear outliers, so we clip up to 2 data points per light curve that exceed 4 from the model. This results in = 8 degrees of freedom. For each EB, we report in Table 1 the eight fitted analytic model parameters and the fit statistics. Excluding the few EBs that exhibit variability or are evolved Roche-lobe filling systems (see below), the goodness-of-fit statistics / = 0.87 - 1.16 indicate the analytic models can adequately describe the EB light curves.

We can measure the physical properties of EBs based solely on the observed photometric light curves (see §3) only if: (1) the binary components are detached from their Roche lobes, (2) the light curves have two well-defined eclipses, and (3) there is no superimposed variability. To be considered well-defined, we require that the 1 uncertainties in the measured eclipse depths and and eclipse widths and are 20% their respective values. These criteria are not satisfied for 91 of the 221 EBs due to a variety of reasons, which we discuss below:

(A) No Secondary Eclipse. For 16 of our EBs, there is no evidence for a secondary eclipse. These EBs may have secondary eclipses that are too shallow and below the sensitivity of the OGLE-III LMC survey, or have eclipse widths that are too narrow and therefore not detected given the cadence of the observations. Most likely, the EBs have a certain combination of , , and as discussed above so that there is only one eclipse per orbit. We list these 16 systems in Category 1 of Table 1, and we show an example in panel A of Fig. 1.

(B) Uncertain. For 32 EBs, both eclipses are observed but one or more of their measured properties are uncertain by more than 20%. This is because one of the eclipses is too shallow and/or too narrow. We group these 32 systems in Category 2 of Table 1. In panel B of Fig. 1, we display an example of a long-period 45 day EB with a secondary eclipse at 0.55 that is too narrow to be accurately measured.

(C) Roche-lobe filling. Three EBs have wide eclipses such that one or both components of the binary must be filling their Roche lobes. We list these three systems in Category 3 of Table 1, and we show an example in panel C of Fig. 1.

(D) Ambiguous Periods. The orbital periods of 23 of our EBs are ambiguous. These 23 EBs can either have twin components 1.0 in nearly circular orbits 0.0 or have half the listed orbital periods and exhibit only one eclipse per very eccentric orbit. Using the orbital periods listed in the OGLE-III LMC catalog, these EBs have primary and secondary eclipses that are nearly identical and separated by almost precisely 50% in orbital phase. Quantitatively, we identify these systems to have values of and uncertainties in eclipse depths, widths, and phases that satisfy:

(2a) | ||||

(2b) | ||||

(2c) |

Given the sensitivity of the data, the observed properties imply the 23 systems have large mass ratios 0.9 with extremely small eccentricities 0.05 (see §3). However, none of the EBs in our sample have eclipse depths that satisfy Eqn. 2a ( 0.9) with secondary eclipse phases 3 10 ( 0.05 - 0.10). Similarly, there is only one EB that satisfies Eqn. 2c ( 0.05) with primary and secondary eclipse depths that are discrepant at the (3 - 10) level ( 0.8 - 0.9). Hence, there are no twin systems in slightly eccentric orbits, and there is only one moderate-mass companion in a nearly circular orbit. The prevalence of 23 twin systems in nearly circular orbits at these moderate orbital periods is therefore highly unlikely. If there is indeed an excess of twins in circular orbits relative to twins in eccentric orbits, our study does not include them. We expect only a few of the 23 EBs that appear to be twins in circular orbits to have the listed orbital periods. The majority of these EBs more likely have orbital periods that are half their listed values, and would therefore exhibit only one eclipse per orbit similar to the systems discussed in (A) above. In panel D of Fig. 1, we show one example where we fold the photometric data with the listed orbital period (in black) and the more plausible scenario that the binary has half the catalog orbital period (in red). We list these 23 EBs in Category 4 of Table 1. We further motivate the removal of these 23 systems in §4 when we show the intrinsic frequency of 0.6 companions with 0.2 is relatively sparse.

(E) Superimposed Variability. Fifteen of the EBs exhibit superimposed variability. Three of these systems are intrinsic variables, two of which (ID-7651 and ID-22929) were already listed as such in the OGLE-III LMC EB catalog. The intrinsic variability is readily apparent in the unfolded light curves. Moreover, the measured intrinsic scatter outside of eclipses is substantially higher than the photometric errors, e.g. 2.8 for ID-3414. We note that a few additional systems with 1.5 - 1.9 may exhibit low-amplitude variations 0.01 mag, but these variations are sufficiently small so as to not to interfere with the light curve modeling. We list the three systems that exhibit definitive intrinsic variability in Category 5 of Table 1. The other 12 EBs exhibit variability in the eclipses themselves, only one of which (ID-17017) was identified as such in the OGLE-III LMC catalog. For these systems, it is possible that more than two bad data points occur near the eclipse. More likely, these 12 EBs display changes in the eclipse depths and/or eclipse phases during the seven years of observations. Apsidal motion due to tidal and relativistic effects are negligible on timescales 7 yrs at these wide orbital separations. Such evolution in the eclipse parameters are most likely caused by orbital motion with a tertiary component (Rappaport et al., 2013). We group these 12 EBs in Category 6 of Table 1, and we display an example in panel E of Fig. 1.

(F) Peculiar. Finally, two EBs have peculiar light curves. ID-343 exhibits a pronounced peak in the folded light curve at = 0.8 between eclipses. This peak may be caused by ellipsoidal modulation in an extremely eccentric orbit. ID-4458, which is shown in panel F of Fig. 1, displays a sinusoidal variation between two eclipses of comparable depth. ID-4458 may contain a hot spot and/or disk, and is similar to the green systems in the top left corner of Fig. 3 in Paper II. We list these two systems in Category 7 of Table 1.

After removing these 91 systems, our well-defined sample contains 130 EBs. We list these 130 systems in Category 8 of Table 1. When necessary, we switch the primary and secondary eclipses to ensure in our well-defined sample. If the epoch of primary eclipse minimum substantially changed from the catalog value in order to satisfy this criterion, we place an asterisk next to our value of in Table 1.

The 130 EBs in our well-defined sample have uncertainties in eclipse depths and and eclipse widths and that are 20% their respective values. The uncertainties in the Gaussian analytic fit parameters have been used only to determine which EBs have detectable and measurable eclipse properties. These uncertainties propagate into our Monte Carlo simulations when we calculate the fraction of binaries that produce detectable eclipses (see §5). The uncertainties in the Gaussian analytic fit parameters are not utilized to calculate the uncertainties in the physical properties of the EBs. Instead, we implement detailed light curve models to measure the values of and uncertainties in the physical model properties, which we now discuss.

## 3. Physical Models (Stage II)

### 3.1. Algorithm

The physical properties of EBs are routinely measured by fitting detailed models to the observational data (Wilson & Devinney, 1971; Prša & Zwitter, 2005; Devor & Charbonneau, 2006; Kallrath & Milone, 2009). Normally, spectroscopic radial velocity observations are required to dynamically measure the component masses and . In the modern era of wide-field photometric surveys, the discovery of EBs is quickly outpacing our ability to obtain follow-up spectra (Devor et al., 2008; Prša et al., 2011a, b). For large EB samples, the physical properties must be inferred based solely on the photometric light curves. MS constraints (Prša & Zwitter, 2005; Kallrath & Milone, 2009) and isochrone fitting (Devor & Charbonneau, 2006) have helped ascertain EB properties from the photometric data. In general, however, these methods lead to large systematic uncertainties and/or solutions that are highly degenerate.

In Paper II, we developed a technique that uniquely and accurately characterizes the intrinsic physical properties of detached EBs with known distances using only the photometric data. The distances to EBs in the LMC are known. In fact, we have already utilized the observed magnitudes and colors to select EBs with B-type MS primaries. For detached EBs with MS primaries, both components are effectively evolving along their respective single-star evolutionary sequences. The photospheric properties of the stellar components, e.g. effective temperatures and , radii and , and luminosities and , therefore depend entirely on the age and component masses and . The systematic uncertainties in the evolutionary tracks are relatively small, e.g. 15% uncertainties in the masses and 30% uncertainties in the ages (see Paper II and §3.3 for further justification and a full assessment of the uncertainties). We can therefore measure the component masses and and ages of detached EBs with known distances based solely on the observed light curve features (Fig. 2).

In our physical models, detached EBs with B-type MS primaries can be uniquely described by nine independent properties. These nine physical model parameters are the orbital period , epoch of primary eclipse minimum , primary mass , secondary mass , age , inclination , eccentricity , argument of periastron , and I-band dust extinction . Given the age and component masses and of the binary, we interpolate the radii and , surface gravities and , effective temperatures and , and luminosities and from pre-MS and MS stellar evolutionary tracks with metallicity = 0.008 (Tognelli et al., 2011; Bertelli et al., 2009). We then use the LMC distance modulus = 18.5 (Pietrzyński et al., 2013), dust reddening law = 0.7 (Cardelli et al., 1989; Fitzpatrick, 1999; Ngeow & Kanbur, 2005), and temperature-dependent color indices and bolometric corrections (Pecaut & Mamajek, 2013) to transform the intrinsic properties of the binary into observed magnitudes and colors. Our physical model parameter space (, , , etc.) of EBs with detached configurations, pre-MS/MS evolutionary constraints, and known distances is quite different than the typical EB parameter space (/, (+)/, etc.) where the distances and evolutionary status of the components are unknown (e.g. Devor & Charbonneau, 2006; Prša et al., 2011a).

Using the physical properties of a binary, e.g., , , , , , , , , , etc., we synthesize photometric light curves with the EB modeling package Nightfall^{1}^{1}1http://www.hs.uni-hamburg.de/DE/Ins/Per/Wichmann/

Nightfall.html. We use the same Nightfall model options adopted in Paper II, e.g. a square-root limb darkening law, default gravity darkening coefficients, model atmospheres, etc., except for three notable distinctions. First, we do not assume circular orbits for our EBs at longer orbital periods, but instead solve for both the eccentricity and periastron angle . Second, we set the albedo of the secondary to = 0.7 and implement one iteration of reflection effects. Considering reflection effects are minuscule for our wider EBs in this study, different treatments of reflection have negligible effects on the synthesized light curves. Finally, we simulate an EB light curve at 1,000 uniformly-spaced discrete orbital phases to ensure narrow eclipses are sufficiently sampled.

Most of our EBs with intermediate orbital periods have eccentric orbits and narrow eclipses. Nightfall and all other EB software packages that account for tidal effects are computationally expensive for eccentric binaries. This is because the three-dimensional photospheric surfaces of the stars need to be recalculated at each of the 1,000 discretely sampled orbital phases. We therefore adapt our algorithm from Paper II to guarantee fast, automated convergence. Namely, we choose initial values for our nine physical model properties that are sufficiently close to the true values to ensure minimization converges quickly to the global solution. The major goal of our algorithm is to synthesize light curves with Nightfall as few times as possible. Our routine can easily be adapted for any population of detached EBs with known distances, and can be used in combination with any EB light curve modeling software.

We decompose our algorithm into three steps.^{2}^{2}2The three steps discussed in this section are not to be confused with the three full stages of our automated pipeline, which classifies EBs (Stage I - §2), fits physical models to the light curves (Stage II - §3), and corrects for selection effects (Stage III - §5). The three steps regarding physical models are all included in Stage II. In Step 1, we select initial values for our nine physical model properties based on the observed light curve features quantified in §2. In Step 2, we make small adjustments in the physical model properties until the analytic model parameters of the synthesized light curve matches those of the observed light curve. In Step 3, we utilize a Levenberg-Marquardt technique, as done in Paper II, to minimize the statistic between the observed and simulated light curves. We elaborate on these three steps below. To help illustrate this procedure, we display in Fig. 3 the light curve of an example EB, ID-2142, and the solutions at the end of each of these three steps.

Step 1. We use the eight analytic model parameters (, , , , , , , ) and observed color from Table 1 to estimate initial solutions for the nine physical model properties. In Fig. 2, we show how the nine observed light curve features can be used to approximate the nine physical properties of the binary. We select the physical parameters and to match the analytic model values. We then estimate and according to the observed phase of the secondary eclipse and the difference in eclipse widths (Kallrath & Milone, 2009, Eqn. 3.1.24 and 3.1.26, see our Fig. 2):

(3a) | ||||

(3b) |

In this study, = 90 if periastron coincides with the observed primary eclipse. For our example ID-2142, and 0.16, indicating 180 and 0.5 - 0.6.

The intrinsic colors of B-type MS stars span a narrow interval 0.3 0.1 (Pecaut & Mamajek, 2013). We therefore initially assume the intrinsic color of an EB to be:

(4) |

where we have accounted for the fact that more luminous B-type MS stars tend to be more massive, hotter, and bluer. The dust extinction is simply estimated from the observed color and our adopted dust reddening law = = 0.7 (see Fig. 2).

We then use the eclipse depths and to approximate the mass ratio = /. For a MS + MS binary in a circular orbit, the ratio of eclipse depths / provides an accurate indicator of the luminosity contrast / (Fig. 2). This luminosity contrast can then be used to infer the mass ratio according to a MS mass-luminosity relation (see Fig. 3 in Paper I). For eccentric orbits, however, the eclipse depth ratio can be modified because the projected distances during primary and secondary eclipses can be different. Nonetheless, deeper eclipses still suggest larger mass ratios. For example, 0.4 mag requires 0.7, regardless of the eccentricity or whether the secondary is a MS or pre-MS star. We use a linear combination of these methods to estimate the mass ratio:

(5) |

where the eclipse depths are in magnitudes.

We next use the observed mean magnitude and sum of eclipse widths to simultaneously measure the primary mass and age . Assuming non-grazing eclipses and standard limb darkening coefficients, the sum of eclipse widths directly provides the relative sum of the radii (+)/. Our EBs occupy a narrow range of magnitudes 16.0 17.6 and therefore span a small interval of total masses = . The orbital separation therefore derives mainly from the known period . We can now use and to determine +. For EBs with B-type MS primaries, we find the following approximation:

(6) |

Given our estimates for and above, we interpolate the stellar evolutionary tracks to determine the primary mass and age that reproduce the sum of the radii according to Eqn. 6 and the observed combined magnitude . Although and both depend on and , they are sufficiently non-degenerate so that we can calculate a unique solution. Namely, the primary mass largely dictates the luminosity and therefore the observed magnitude , while the age primarily determines the radii and therefore the observed eclipse widths (see Fig. 2).

Finally, we select an inclination that approximately reproduces the observed primary eclipse depth (Fig. 2). From our estimates of , = , and , we interpolate the radii and and effective temperatures and from stellar evolutionary tracks. In this step only, we ignore limb darkening and colors of the two stars, and instead assume the stars are uniformly illuminated grey disks (see Paper I). We assume the surface brightnesses of the disks are proportional to the stellar temperatures, i.e. the Rayleigh-Jeans law, because we are observing at relatively long wavelengths in the near-infrared I-band. Using these approximations, we calculate the eclipsed area of the primary at the time of primary eclipse based on the observed primary eclipse depth:

(7) |

Given the eclipsed area and stellar radii and , we then determine the projected distance between the two stars at . The actual physical separation at primary eclipse is already known via (Kallrath & Milone, 2009, Eqn. 3.1.36 evaluated at geometric phase = 0):

(8) |

where derives from our estimates of , , and according to Kepler’s third law, and and are approximated from Eq. 3. Hence, the inclination simply derives from cos = /. We limit our initial approximation of the inclination to the interval = 86.5 - 89.5.

We now have initial estimates for the nine physical model properties. We emphasize that Eqns. 3 - 7 are simple approximations, and that the true values of , , , , , , and may substantially differ from the initial values estimated here. We simply use these estimates as initial parameters in our fitting routine in order to minimize the number of iterations and accelarate convergence toward the final solution (see below and §3.2).

In the top panel of Fig. 3, we compare the I-band and V-band light curves of ID-2142 to a simulated Nightfall model using the values of the nine physical model properties at the end of Step 1. The model matches key features of the observed light curve, but there are three noticeable differences. First, the simulated phase of the secondary eclipse does not match the observations; recall that Eqn. 3 is an approximation. Second, the simulated color is bluer than the observed , indicating we underestimated the dust reddening in our initial step. Finally, the simulated eclipses are slightly deeper than the observed because the more accurate Nightfall model accounts for limb darkening and color effects. This suggests the actual inclination is smaller and/or the mass ratio is slightly different. We correct for these visible discrepancies in the following step.

Step 2. Using the nine physical properties from Step 1, we synthesize an I-band light curve with Nightfall. We then fit the simple analytic model of Gaussians (Eqn. 1) as done in §2 to the simulated Nightfall light curve. In this manner, we measure the analytic parameters of the Nightfall model, e.g. , , , etc.

We adjust the properties in our physical models according to the differences between the simulated and observed analytic model parameters. The adjustments are motivated as follows. If the modeled eclipse widths are wider than the observed , we select a slightly younger age (and vice versa). We increase the dust extinction if the simulated color is too blue. If the modeled primary eclipse is too deep while the modeled secondary eclipse matches observations or is too shallow, we increase the mass ratio and decrease the inclination . However, if both simulated eclipses are too deep (or both too shallow), we only decrease (increase) the inclination . Finally, we adjust according to the position of and differences in the secondary eclipse phases and . In this step, we fix , , and to the values determined in Step 1. Finally, we interpolate from the stellar evolutionary tracks based on the observed mean magnitude and the revised values for , , and .

When adjusting our physical model properties, we choose step sizes that scale with the differences between the observed and simulated analytic model parameters. After making these adjustments, we synthesize another I-band light curve with Nightfall. We iterate this step until all the analytic model parameters of the simulated and observed light curves match within a small tolerance level. In the middle panel of Fig. 3, we show our solution for ID-2142 at the end of Step 2 after five iterations. We therefore required only six Nightfall light curve simulations during this middle step.

Step 3. This final step is essentially the procedure outlined in Paper II. We calculate the photometric correction factors and in both bands. Starting with initial model properties determined at the end of Step 2, we utilize a Levenberg-Marquardt technique (MPFIT, Markwardt, 2009) to minimize the statistic between the simulated and observed light curves. The Levenberg-Marquardt MPFIT algorithm operates by independently varying each of the nine physical model properties from the previous solution. The routine then measures the resulting deviations between the data and models, and then calculates a new solution. This step therefore requires ten Nightfall simulations per iteration. As in Paper II, we simultaneously fit the I-band and V-band light curves. We clip up to 3 data points that exceed 4 from the best-fit model. This results in = 9 degrees of freedom.

In the bottom panel of Fig. 3, we display our final solution for ID-2142 after four iterations of the Levenberg-Marquardt MPFIT routine. We therefore simulated light curves with Nightfall a total of 40 times in Step 3. The physical model properties changed only slightly during this final step. In fact, for ID-2142, the variations were all within the uncertainties of the physical model parameters. We emphasize that Steps 1 and 2 were crucial in guarenteeing rapid convergence toward the final solution in Step 3. Without them, this last step would have required many additional iterations or may have converged to a local minimum.

We utilize this automated procedure for all 130 detached EBs in our well-defined sample. We present our fitted model parameters, physical properties, and fit statistics for these systems in Table 2. For MS binaries in circular orbits, the deeper primary eclipse at time always corresponds to the smaller, cooler, less massive secondary passing in front of the larger, hotter, more massive primary. For eccentric orbits, however, the situation can be reversed depending on the combination of , , and . Indeed, for 18 EBs in our well-defined sample, we determined solutions such that the less massive component was eclipsed at time . To avoid confusion in nomenclature, we list properties in Table 2 according to the primary “p” and secondary “s” eclipse features. Namely, , , and correspond to the component that was eclipsed at the epoch of primary eclipse , and , , and correspond to the component that was eclipsed at the secondary eclipse phase . In the text, we refer to primary mass = max{, }, secondary mass = min{, }, mass ratio = /, etc.

We measure primary masses = 3.6 - 13.9 M, which nearly encompasses the full mass range of B-type MS stars. We determine mass ratios across the interval = 0.20 - 1.00, which confirms the OGLE-III observations are sensitive to EBs with low-mass companions. Our measured dust extinctions cover = 0.10 - 0.58 mag, which is consistent with the range of extinctions found in Paper II. Finally, we determine ages = 0.5 - 190 Myr that span more than two orders of magnitude. We further discuss the EB physical properties, and their interrelations, in §4.

Eleven of the 130 EBs have modest fit statistics / = 1.10 - 1.14, i.e. probabilities to exceed of 0.01 - 0.05 given 530 degrees of freedom. Seven of these EBs are extremely young with estimated ages 0.8 Myr (IDs 5153, 7560, 10422, 13418, 16711, 22691, and 22764). The components in these EBs have small radii, as demonstrated by their narrow eclipses (Eqn. 6), and are therefore consistent with the zero-age MS. The systematic uncertainties in the stellar evolutionary tracks are larger at these young ages, especially considering some of the secondaries may still be pre-MS stars (see Paper II). Three of the 11 EBs with modest fit statistics have primaries at the tip of the MS (IDs 91, 20746, and 21518), as indicated by their wide eclipses. Again, the stellar evolutionary tracks are uncertain at the tip of the MS just prior to the rapid expansion toward the giant phase. The one last EB with a poor physical model fit (ID-17569) has = 1.11 and 0.02. Considering our large sample of 130 EBs, we naturally expect 1 - 3 of these EBs with modest fit statistics. The remaining 119 EBs in our well-defined sample have good fit statistics 0.93 / 1.09. This is testament that the nine independent physical model properties can adequately describe detached EBs with known distances and MS primaries.

###
3.2. Comparison between Initial Estimates

and Final Solutions

In the following, we compare the initial estimates for , , and + in Step 1 according to Eqns. 3, 5, and 6, respectively, to the final solutions in Step 3 from fitting detailed Nightfall light curve models to the data. We can then address the systematic uncertainties in our initial estimates and further justify the mapping between the basic EB light curve parameters to the physical model properties.

For the 130 EBs in our well-defined sample, we compare the initial values of the eccentricities determined from the secondary eclipse phases and eclipse widths and (Eqn. 3) to the final Nightfall solutions (top panel of Fig. 4). The initial estimates agree quite well with the true final values. The rms scatter between the two is only = 0.03. This validates that Eqn. 3 is more than sufficient for starting purposes in our fittng routine. We note that the few systems that change by more than 0.07 between Steps 1 and 3 have narrow, poorly sampled eclipses so that it is more difficult to precisely measure and .

Similarly, in the bottom panel of Fig. 4, we compare the initial estimates of the mass ratios determined from the eclipse depths and (Eqn. 5) to our final values obtained from Nightfall light curve fittings and minimizations. Although the population as a whole shows rough agreement between the solutions at the ends of Steps 1 and 3, individual systems can substantially deviate from the initial estimates. For example, an EB with an initial estimate of 0.6 may actually have a mass ratio anywhere in the interval = 0.3 - 1.0. The rms deviation between the initial and final solutions is = 0.12, or / 20% the respective values. If we had randomly chosen mass ratios in Step 1 while keeping the other initial estimates unchanged, the Nightfall light curve solutions in Step 3 would still converge to the same final values. We simply find that by adopting Eqn. 5 in Step 1 to provide initial estimates for , the number of iterations in Steps 2 and 3 are dramatically reduced.

Finally, we evaluate the discrepancies between the values of + estimated from the sum of eclipse widths + and orbital periods according to Eqn. 6 to the final Nightfall solutions. We measure an rms deviation of (+) = 1.2 R, or (+) / (+) 20% the respective values. The coefficient in Eqn. 6 should therefore be 7.0 1.2 R, valid only for EBs with B-type MS primaries. As with the mass ratios , the approximations for + based on the observed light curve parameters are imprecise but sufficiently accurate to provide initial conditions for our fitting routine.

As mentioned in §3.1, + is primarily an indictor of age of an EB in our sample rather than the component masses and/or . For example, at age = 5 Myr, the sum of the stellar radii must be contained on the interval + 4.4 - 8.7 R given any combination of and 0.2 that satisfies our magnitude limits 16.0 17.6 and 0.25 0.20 and measured range of dust extinctions 0.1 - 0.5 mag. Meanwhile, at age = 100 Myr, the sum of the radii are systematically larger and confined to the interval + 5.3 - 12.3 R given the same photometric requirements. Hence, EBs in our sample with + 5.3 R must be relatively young while those with + 8.7 R must be relatively old. We initially estimated + in Step 1 from the observed sum of eclipse widths + according to Eqn. 6. Although not as accurate as the final solutions, Eqn. 6 provides a model-independent measurement of +. The sum of radii + estimated from Eqn. 6 is therefore a robust and model-independent indicator of age .

In Fig. 5, we display the eccentricities measured in Step 1 from Eqn. 3 as a function of the approximate sum of stellar radii + estimated in Step 1 from Eqn. 6. Both sets of parameters are model independent and based solely on the observed light curve features. According to a Spearman rank correlation test, we find that the approximate values of and + are anticorrelated ( = 0.18) at a statistically significant level ( = 0.04). This suggests that EBs with larger components, which are systematically older, favor smaller eccentricities. We also compare the 19 EBs with approximate + 5.3 R, which must be relatively young, to the 27 EBs with + 8.7 R, which must be relatively old. According to a K-S test, we find these young and old populations of EBs have distributions of eccentricities that are discrepant with each other at the = 0.02 significance level. The anticorrelation between +, which is an indicator of age , and is therefore statistically significant, robust, and model independent. In §4, we further investigate this anticorrelation between and based on the more accurate final solutions obtained from the Nightfall light curve models.

### 3.3. Uncertainties

We now analyze the uncertainties in the final solutions of our Nightfall light curve models (see also Paper II). For each system, we utilize MPFIT (Markwardt, 2009) at the end of Step 3 (§3.1) to calculate the measurement uncertainties. For all 130 well-defined EBs, the nine physical model parameters have unique solutions and finite measurement uncertainties. Some of the model parameters, however, have solutions that are correlated with each other. In addition, uncertainties in the dust reddening law, stellar evolutionary tracks, bolometric corrections, and Nightfall light curve models can lead to large systematic uncertainties in the physical model parameters. In the following, we fully investigate the measurement uncertainties, parameter correlations, and systematic uncertainties in the context of a specific example EB, ID-2142. We then determine the median total uncertainties of each model parameter (Eqns. 9 - 17) for the entire population of 130 well-defined EBs.

For our example EB, ID-2142, we explore the physical parameter space via a Markov chain Monte Carlo (MCMC) technique. Starting with our final solution at the end of Step 3, we implement a Metropolis-Hastings “random walk” MCMC algorithm to generate and select steps in our phase space of nine physical model parameters. At each proposed step, we synthesize a Nightfall light curve model given the proposed nine physical model parameters. The probability of accepting the proposed step is determined by evaluating the difference in the statistic between the proposed step and the current solution. Obviously, if 0, the proposed step is always taken. If the proposed step is rejected, the step length is effectively zero, i.e., the previous solution is counted again. We generate proposed steps according to a Gaussian distribution with a fixed standard deviation for each of the nine physical model parameters. We choose the standard deviation in the step sizes so that approximately one-third of the proposed steps are accepted. We simulate 32,000 proposed steps and light curves with Nightfall, which exceeds the total number of models generated in §3.1 used to fit solutions for all 130 well-defined EBs! It is therefore quite computationally expensive to calculate robust measurement uncertainties and correlations between model parameters for an individual EB with this MCMC algorithm. The distribution of the 12,000 accepted steps and 20,000 repeated solutions provide the nine-dimensional joint probability distribution for the physical models. For each of the nine physical model parameters, we marginalize across the other eight parameters to calculate the one-dimensional probability density function. We also compute the two-dimensional joint probability distributions for each of the = 36 parameter combinations.

In Fig. 6, we display the one-dimensional probability distributions for the nine physical model parameters (diagonal panels) and the two-dimensional joint probability distributions for the 36 parameter combinations (off-diagonal panels). Although some of the parameters are mildly to significantly correlated with each other, the measurement uncertainties are finite for all nine physical model parameters. The MCMC technique confirms the uniqueness and non-degeneracy of the physical model solutions. Moreover, the measurement uncertainties we determined from the robust MCMC algorithm are consistent with the measurement uncertainties we evaluated with the MPFIT routine. We can therefore rely on the MPFIT measurement uncertainties we calculated for all 130 well-defined EBs.

The uncertainties in the orbital parameters and are solely due to the measurement uncertainties and dictated by the sensitivity and cadence of the OGLE-III LMC observations. The solutions for and are therefore independent of the other seven model parameters (note the fairly circular contours in the first and second columns of panels in Fig. 6). For ID-2142, we measure 1 uncertainties of 0.0001 days and 0.004 days. We find the median 1 uncertainties for the entire population of 130 well-defined EBs to be:

(9) |

(10) |

Note that our example ID-2142 has slightly smaller uncertainties than average because it is relatively bright and its eclipses are well sampled.

As with and , the uncertainties in and are primarily determined by the sensitivity and cadence of the OGLE-III LMC observations. The solutions for and are therefore independent of the other parameters, but are slightly correlated with each other (see last two rows in Fig. 6). The eclipses are sufficiently sampled to easily break this degeneracy. For ID-2142, we calculate a 95% confidence interval of = 175 - 180. Note that we measured eclipse widths = 0.0047 = 0.0050, also indicating 180 according to the approximations in Step 1 (Eqn. 3). Based on the Nightfall light curve models, we calculate formal 1 measurement uncertainties of 0.001 and 1.4 for ID-2142.

For ID-2142 and some other EBs in our sample, the measurement uncertainties 0.005 and 1.5 are extremely small. Nightfall treats each stellar component as a three-dimensional polyhedral mesh with a finite number of flat surfaces. We suspect this finite resolution limits the true sensitivity to systematic uncertainties of 0.005 and 1.5. In any case, the measurement uncertainties and increase and dominate the total uncertainties as the eccentricities decrease. We measure median total uncertainties of 0.02 and 4 for 0.5, 0.03 and 10 for 0.3, and 0.05 and 20 for 0.1. Obviously, the periastron angle is not defined, and therefore not constrained, if the orbits are circular. For the entire population of 130 well-defined EBs, we find the following relations adequately describe the median total uncertainties:

(11) |

(12) |

Solutions for the remaining five parameters , , , , and are all correlated with each other (see Fig. 6). Moreover, unlike , , , and , which have relatively symmetric Gaussian errors, the probability density functions of , , , , and are mildly to significantly asymmetric. The three parameters , , and are especially correlated along the observed magnitude . In other words, solutions with more massive primaries require younger ages and higher extinctions to produce the same observed I-band flux. The secondary mass is also anticorrelated with . Finally, the inclination mildly depends on the three parameters , , and that are significantly correlated with each other.

Although , , , and are correlated with each other, there is sufficient information in the observed light curves and our constraints (e.g., distance, evolutionary tracks, dust reddening law) to break the degeneracies and provide unique solutions (see also §3.1). For example, if we were to fix the primary mass at = 11.0 M (i.e., the 2.5 upper limit according to the probability density function in Fig. 6), the other parameters would converge to = 7.1 M, = 3.6 Myr, = 86.35, and = 0.45 mag with a fit statistic that is = 6.7 larger than the best-fit solution. For this larger primary mass, there is no combination of and that can satisfactorally reproduce the observed magnitude and color . Similarly, if we were to fix the primary mass at = 9.2 M (i.e., the 2.5 lower limit), the other parameters would converge to = 6.1 M, = 13.3 Myr, = 86.08, and = 0.39 mag with a fit statistic that is = 5.8 larger than the best-fit solution. In this case, the component masses both decrease by 12% (to maintain the same ratio of eclipse depths ), and so the orbital separation decreases by 4% according to Kepler’s third law. The relative sum of the radii ()/, which derives directly from the sum of eclipse widths , is measured to 1% precision in our Nightfall light curve models. If decreases by 4%, then must also decrease by 4%. According to the MS stellar evolutionary tracks, if and decrease by 12%, then the radii and decrease by 9% given the same age = 7.2 Myr. Hence, the age must increase to = 13.3 Myr so that the sum of radii only decreases by 4%. If the masses decrease, the radii decrease, and the age increases, then the temperatures and both decrease according to the stellar evolutionary constraints. However, if and both decrease, it is difficult to maintain the same values of and with only one free extra paramter . Hence, there is no combination of , , and that can satisfactorially reproduce the observed values of , , , and if = 9.2 M. This line of reasoning holds for all EBs in our sample, and so the physical model parameters will always have unique solutions with finite measurement uncertainties. For ID-2142, we measure formal 1 measurement uncertainties of 0.04 0.5 M, 0.07 0.4 M, 0.25 1.8 Myr, 0.1, and 0.01 0.006 mag. We find similar percentage measurement uncertainties in these parameters for the 130 well-defined EBs in our sample.

The systematic uncertainties in , , , and can be considerably larger and derive from a variety of sources. We first investigate the systematic uncertainties in the adopted bolometric corrections. Our B-type MS primaries and secondaries span a large range of temperatures 10,000 - 30,000 K and therefore a broad interval of bolometric corrections = M M 3.0 - 0.3 mag (Pecaut & Mamajek, 2013). For the hottest stars in our sample with 30,000 K and 3.0, the bolometric corrections are uncertain by 0.2 mag, i.e. / 7% (Bertelli et al., 2009; Pecaut & Mamajek, 2013). To propagate this systematic uncertainty into our solution for ID-2142, we decrease the absolute magnitudes of the bolometric corrections by 7% and repeat our fitting routine from §3.1. The Nightfall light curve models now converge to a final solution of = 9.3 M, = 6.3 M, = 7.5 Myr, = 86.35, and = 0.40 mag. The main effect of decreasing is to decrease the masses and . This is because more of the flux is radiated in the optical and so the component luminosities need to be reduced to maintain the same observed magnitude . Fortunately, the mass ratio = / is not significantly affected by the uncertainties in the bolometric corrections. The decrease in masses lead to slightly longer ages (to maintain the observed eclipse widths), higher inclinations (to maintain the observed eclipse depths), and lower extinctions (to maintain the observed color). The systematic uncertainties in the physical model parameters due to the uncertainties in the bolometric corrections are therefore = 0.11 = 1.1 M, = 0.08 = 0.5 M, = 0.04 = 0.3 Myr, = 0.18, and = 0.1 = 0.04 mag. Because the primary mass is mainly dictated by the observed and the bolometric corrections, we expect similar percentage systematic uncertainties in , , , , and for the other EBs in our sample.

We next propagate the uncertainties in the intrinsic colors, observed colors, and dust reddening law. The uncertainty in the intrinsic colors of B-type MS stars are 0.01 - 0.02 mag (Pecaut & Mamajek, 2013), the color calibrations of stars in the OGLE-III LMC database are also uncertain by 0.01 - 0.02 mag (Udalski et al., 2008), and the coefficient in our adopted dust reddening law = 0.70 has a 10% uncertainty (Cardelli et al., 1989; Schlegel et al., 1998; Fitzpatrick, 1999; Ngeow & Kanbur, 2005). The systematic uncertainty in the dust extinction due to dust/color uncertainties is therefore = max{0.02 mag, 0.1}. To confirm this estimate, we replace the dust reddening law with = 0.63 in our models for ID-2142 and repeat our fitting routine from §3.1. As expected, we measure = 0.48 mag, i.e. the dust extinction increased by = 0.1 = 0.04 mag, while the other parameters do not vary beyond the measurement uncertainties.

We finally investigate the uncertainties in the stellar evolutionary tracks, including the effects of metallicity and rotation. We replace the Z = 0.008 tracks from Bertelli et al. (2009) with the Z = 0.006 non-rotating models from Georgy et al. (2013). We refit ID-2142 and measure = 10.5 M, = 6.8 M, = 9.0 Myr, = 86.38, and = 0.45 mag. Hence, the systematic uncertainties in the stellar evolutionary models, including our ability to interpolate between the tracks, dominates the uncertainty in the age = 0.26 = 1.9 Myr and inclination = 0.21. We then replace the evolutionary tracks with the Z = 0.006 tracks from Georgy et al. (2013) that are rotating on the zero-age MS at = 50% the critical break-up velocity. We note that 80% of B-type MS stars are rotating at 0.5 250 km s (Abt et al., 2002; Levato & Grosso, 2013), and our EBs with intermediate orbital periods = 20 - 50 days may have tidally evolved toward slower rotational velocities. B-type MS stars initially rotating at / = 0.5 have equatorial radii that are only (3 - 4)% larger than their polar radii, but MS lifetimes that are 20% longer (Georgy et al., 2013). It is therefore the differences in the evolutionary tracks of stars with rotation, not the distortions in their shapes, that can significantly affect our model solutions. We refit ID-2142 with the rotating non-synchronized stellar models, and measure = 10.5 M, = 6.9 M, = 8.6 Myr, = 86.30, and = 0.45 mag. For ID-2142, the differences between the non-rotating and rotating tracks from Georgy et al. (2013) are within the measurement uncertainties. This is because the tracks with / = 0.5 do not significantly deviate from their non-rotating counterparts until the ages reach 0.8 the non-rotating MS lifetimes. For ID-2142 and the majority of EBs in our sample with primary ages 0.8, the uncertainties due to the effects of rotation are negligible. For the few systems that are extremely young ( 1 Myr) or old ( 0.8), we expect slightly larger systematic uncertainties in the ages and masses.

By adding the measurement uncertainties and various systematic uncertainties above in quadrature, we estimate the total median 1 uncertainties for the 130 well-defined EBs to be:

(13) |

(14) |

(15) |

(16) |

(17) |

For these five parameters, the uncertainties are dominated by the systematic uncertainties in the bolometric corrections, dust reddening law, and evolutionary tracks. For older EBs with primary ages 80% their MS lifetimes, the total uncertainties in the masses 0.2 and 0.2 and ages 0.45 are slightly larger due to the effects of rotation.

Some of our EB light curve model solutions can be biased due to contamination with a third light source, e.g., a tertiary companion or a background/foreground object along similar lines of sight. In Papers I and II, we estimated that only 10% of our B-type MS EBs in the LMC can be contaminated by a third light source that is bright enough to significantly contribute to the systematic uncertainties. Unlike the previously discussed sources of systematic uncertainties that contribute to all 130 well-defined EBs, contamination by a third light source affects only a small subset of our sample.

In addition to calculating the uncertainties for the nine independent physical model parameters, we also estimate the uncertainties in the dependent physical properties. The total uncertainties in and are 15% but slightly correlated with each other (see above). The total median uncertainty in the mass ratio is therefore = max{0.03, 0.12}. Because the quantity (+)/ is precisely constrained from the observed eclipse widths, the uncertainties in , , and mainly derive from the uncertainties in and according to Kepler’s third law. We measure 0.07 0.3R, 0.07 0.2 R, and 0.06 6 R. Finally, given the (20 - 30)% uncertainties in the luminosities (primarily due to uncertainties in the bolometric corrections) and the 7% uncertainties in the radii, the uncertainties in the temperatures are 8% according to the Stefan-Boltzmann law. Hence, the total median uncertainties are 0.08 1,500 K and 0.08 1,100 K.

## 4. EB Trends

In our sample of 130 EBs, several trends and correlations exist among the nine physical model properties. Most of these trends are caused by geometrical and evolutionary selection effects in our magnitude-limited sample of EBs. We correct for these selection effects in the third stage of our pipeline (§5). Two correlations, however, are intrinsic to the population of binaries with B-type MS primaries. In this section, we first discuss these two empirical relations we uncovered from the data, and then we explain the trends that are caused by selection effects.

In Fig. 7, we display the measured I-band dust extinctions as a function of age for the 130 well-defined EBs. These two parameters are anticorrelated (Spearman rank correlation coefficient = 0.34) at a statistically significant level (probability of independence = 810). We fit a log-linear trend to the total population of 130 EBs (green line in Fig. 7):

(18) |

The slope in the above relation may be biased toward negative values due to a photometric selection effect in our magnitude-limited sample. Specifically, EBs that are intrinsically bluer and more luminous systematically contain younger, short-lived, more massive primaries. These blue, luminous, younger EBs may therefore require larger dust extinctions and reddenings to satisfy our photometric selection criteria (and vice versa). In Fig. 8, we show the measured absolute magnitudes and intrinsic colors as a function of dust extinction for our 130 well-defined EBs. We also display our photometric selection criteria based on the observed magnitudes 16.0 17.6 and observed colors 0.25 0.20 (green lines). Indeed, there are several intrinsically red, low-luminosity, older EBs with 1.2 that are in our sample only because they have small dust extinctions 0.2 mag. If they were to have slightly higher dust extinctions, they would fall below our selection limit of 17.6 (see Fig. 8).

In §5, we account for our photometric selection criteria when analyzing all 130 EBs in our well-defined sample. Here, we correct for photometric selection effects by further culling our sample according to the intrinsic properties of and . To obtain an unbiased subsample, we can choose EBs across any interval of and that also satisfies our selection criteria on observed magnitudes and colors. To retain most of the sample, we select the regions enclosed by 2.63 1.35, 0.341 0.115, and 0.13 (mag) 0.45 (red lines in Fig. 8). The 98 EBs that satisfy these extra selection criteria (filled blue systems in Fig. 8) represent an unbiased sample relatively free from photometric selection effects.

Even within this unbiased sample of 98 EBs, the intrinsic colors and dust extinctions are still anticorrelated ( = 0.25) at a statistically significant level ( = 0.02). As can be observed in the bottom panel of Fig. 8, there are relatively few intrinsically blue systems 0.30 with small dust extinctions 0.2 mag. Similarly, there are few intrinsically redder EBs 0.15 with large dust extinctions 0.4 mag. Intrinsically bluer EBs contain hot primaries that are systematically more massive, short-lived, and younger. Hence, the anticorrelation between age and dust extinction is real.

In Fig. 8, we also display and for the 98 EBs (filled blue) in our unbiased subsample. Although not as prominent, the ages and dust extinctions for the 98 EBs in our unbiased sample are still anticorrelated ( = 0.23) at a statistically significant level ( = 0.02). For example, there is a complete absence of EBs with 0.2 mag at 15 Myr. In contrast, there are many EBs in our unbiased sample with = 0.1 - 0.2 mag at 15 Myr. The intrinsic anticorrelation between and in our unbiased sample demonstrates a relationship between dust content and ages of stellar environments. Young EBs, and young B-type MS stars in general, with 1 Myr are embedded in dusty envelopes and/or molecular clouds with photometric extinctions 0.33 mag. Meanwhile, older EBs with 100 Myr reside in less attenuating environments with 0.26 mag. We fit a log-linear trend to the unbiased sample of 98 EBs:

(19) |

valid for 0.5 Myr 200 Myr (red line in Fig. 8). Even after accounting for selection effects, the value of and measurement uncertainty in the slope 0.029 0.011 is still inconsistent with zero at the 2.6 confidence level. This is a similar probability of significance based on the Spearman rank test above (probability of no correlation = 0.02 between and ). The 30% systematic uncertainty in the ages and 10% systematic uncertainty in the extinctions propagate into Eqn. 19. The values of and total uncertainties are therefore 0.029 0.014 for the slope and 0.33 0.04 mag for the mean dust extinction at 1 Myr. The rms in the measured dust extinctions around the above relation is = 0.08 mag.

It had been previously known that younger early-type stars in the LMC experience slightly higher dust extinctions than late-type stars (Zaritsky, 1999; Zaritsky et al., 2004). In the present study, we have measured the relationship between age and dust extinction . Quantifying age-dependent dust extinctions is crucial when analyzing the spectral energy distributions of unresolved stellar populations in distant galaxies (Panuzzo et al., 2007; Silva et al., 2011). Young O- and B-type stars, which dominate the ultraviolet component in star-forming galaxies, will experience systematically higher dust extinctions than the older, redder stars. To accurately constrain the star-formation histories of these galaxies, it is imperative to account for age-dependent dust extinctions. We note that different galaxies and stellar populations will have slightly different dust extinctions as a function of age. Nonetheless, our empirical age-extinction relation (Eqn. 19) can provide insight when calibrating models of unresolved stellar populations.

Zaritsky et al. (2004) found that the dust extinction distribution toward young, hot stars in the LMC peaks at 0.25 mag with a long tail toward higher values. This is consistent with our total population of 130 EBs with B-type MS primaries (see Fig. 9). By dividing our EB population into young ( 12 Myr) and old ( 12 Myr) subsamples, we find that both subsamples can be fitted with simple Gaussians centered at 0.34 mag and 0.25 mag, respectively. Hence, the non-Gaussian distribution of dust extinction may be simply due to a selection effect with age. The very young EBs, which represent a small fraction of the total population, occupy the long tail toward large dust extinctions. Meanwhile, the long-lived EBs, which comprise the majority of the sample, form the peak in the distribution at 0.25 mag.

We now examine the second physically-genuine trend in our EB population. In Fig. 10, we show the measured eccentricities as a function of age for the 130 EBs in our well-defined sample. The eccentricities and ages are anticorrelated (Spearman rank correlation coefficient = 0.39) at a statistically significant level (probability of no correlation = 510).

This observed anticorrelation is primarily because eccentricities decrease with time due to tidal evolution. The observed trend may be accentuated by a secondary effect, whereby EBs with more massive, short-lived primaries favor larger eccentricities. However, this relation between primary mass and eccentricity cannot fully explain the observed anticorrelation between and . For example, the eccentricities and ages of the 32 EBs with massive primaries 8.5 - 13.9 M are still anticorrelated ( = 0.44) at a statistically significant level ( = 0.01). Similarly, the 98 less massive EBs with 3.6 - 8.5 M have eccentricities and ages that are anticorrelated ( = 0.28) at a statistically significant level ( = 0.005). Although EBs with early-B primaries may be born with systematically larger eccentricities, the anticorrelation between age and eccentricity is dominated by tidal evolution and is observed in both early-B and late-B MS subsamples.

For late-type stars with 1.3 M, orbital energy is most efficiently dissipated into the interior of the stars via convective eddies in the stellar atmospheres (Zahn, 1977; Hut, 1981; Zahn, 1989; Hurley et al., 2002). This equilibrium tide model for convective damping has been tested against observations of late-type binaries in various environments with different ages (Meibom & Mathieu, 2005). For more massive stars 1.3 M with radiative envelopes, such as our B-type MS stars, tides operate dynamically via oscillations in the stellar interiors (Zahn, 1975; Hurley et al., 2002). By estimating the ages of 130 early-type EBs, we have measured the evolution of binary eccentricities due to dynamical tides with radiative damping.

The slope of the observed age-eccentricity anticorrelation provides insight into the tidal evolution of highly eccentric binaries. We fit a log-linear trend to the observations (red line in Fig. 10):

(20) |

The value of and measurement uncertainty in the slope is 0.14 0.03. Hence, the slope is negative at the 5 confidence level, similar to the statistical significance determined from the Spearman correlation test above. Again, systematic uncertainties in the ages and eccentricities contribute to the uncertainties in the coefficients in Eqn. 20. After calculating the total uncertainties, we find the mean eccentricity at 1 Myr is 0.53 0.05 while the slope is 0.14 0.05. The rms scatter in the measured eccentricities around the above relation is 0.16.

The intercept in Eqn. 20 implies a circularization timescale of 5 Gyr for our EBs with B-type MS primaries and moderate orbital periods 20 - 50 days. However, tidal damping is not as efficient when the orbits become less eccentric (Hut, 1981). The true circularization timescale may therefore be longer if the age-eccentricity relation flattens beyond 200 Myr. Conversely, older EBs have systematically larger components (see Fig. 5), and so tidal damping may become more efficient as the primary fills a larger fraction of its Roche lobe. In any case, these short-lived B-type MS primaries will expand beyond 10 R and evolve toward the giant branch long before the orbits are completely circularized.

Our young EBs with generally large eccentricities experience extreme tidal forces. In fact, a few of the EBs with 0.6 in our sample have modest Roche-lobe fill-factors 0.3 at periastron. Tidal evolution of highly eccentric binaries is quite complicated, especially considering second-order effects and non-linear terms can become quite important (Hut, 1981). A full analysis of tidal evolution in our EB sample is therefore not within the scope of the present study. Nonetheless, the observed age-eccentricity anticorrelation provides a constraint for models of tidal evolution in highly eccentric early-type binaries.

In Fig. 11, we display the cumulative distribution function of the eccentricities for the 128 EBs with 0.68 (green). We do not consider the two EBs with 0.71 and 0.77 because highly eccentric binaries are not complete in our EB sample (see below and §5). Moreover, as discussed above, binaries with = 20 - 50 days and = 0.7 - 0.8 nearly fill their Roche lobes at periastron, and are expected to evolve toward smaller eccentricities on rapid timescales. In Fig. 11, we also divide our sample into the 91 old EBs with 10 Myr (red) and 37 young EBs with 10 Myr (blue). Using a maximum likelihood method, we fit a power-law eccentricity probability distribution to the observed EBs. We measure = 0.1 0.2, 0.1 0.2, and 0.8 0.3 for the total, old, and young EB samples, respectively. Our total population of EBs ( = 0.1 0.2) is consistent with the flat distribution ( = 0) observed by Abt (2005) for his sample of binaries with B-type MS primaries and intermediate orbital periods.

If the orbital velocities and energies of a binary population follow a Maxwellian “thermal” probability distribution, then the eccentricity probability distribution = will be weighted toward large eccentricities (Ambartsumian, 1937). Such a population of eccentric and thermalized binaries would suggest the binaries formed through dynamical interactions, either through tidal / disk capture, dynamical perturbations in a dense cluster, three-body exchanges, and/or Kozai cycles with a tertiary companion (Heggie, 1975; Pringle, 1989; Turner et al., 1995; Kroupa, 1995; Kiseleva et al., 1998; Naoz & Fabrycky, 2014). Surprisingly, the observed population of 37 young EBs ( = 0.8 0.3) is fully consistent with a thermal eccentricity probability distribution ( = 1; dashed black line in Fig. 11). This indicates that massive binaries with intermediate orbital periods formed via dynamical interactions on rapid timescales 5 Myr.

Previous observations of spectroscopic (Duquennoy & Mayor, 1991) and visual (Harrington & Miranian, 1977) solar-type binaries have indicated a thermal eccentricity distribution. However, these studies recovered the thermal eccentricity distribution only after applying large and uncertain correction factors for incompleteness. In both the spectroscopic and visual binary surveys, the raw samples were weighted significantly toward smaller eccentricities relative to the thermal distribution. In addition, more recent and complete observations of solar-type (Abt, 2006; Raghavan et al., 2010) and early-type (Abt, 2005) binaries at intermediate orbital periods have revealed a uniform eccentricity distribution that is clearly discrepant with a thermal distribution.

Our raw sample of young early-type EBs is only slightly biased toward small eccentricities. In fact, the small excess of young EBs with 0.1 - 0.3 relative to the thermal distribution in Fig. 11 would be reduced after correcting for selection effects. In other words, we expect even better agreement between our sample of young early-type EBs and the thermal eccentricity distribution after considering observational biases (§5). By choosing only the EBs with young ages, we have probed the initial binary properties of massive stars shortly after their formation. For the first time, we have directly observed the theoretical thermal eccentricity distribution before tides have dramatically reduced the eccentricities.

In the following, we compare other physical model parameters and examine additional trends that could be caused by observational biases. We use these observed distributions to further justify our selection criteria in §2. We also motivate the necessity for incompleteness corrections and Monte Carlo simulations, which we perform in §5.

We display the measured eccentricities as a function of mass ratio in Fig. 12. A Spearman rank test reveals no statistically significant correlation ( = 0.25). The mass ratios of early-type binaries are independent of their eccentricities at intermediate orbital periods = 20 - 50 days.

In §2 (see item D), we removed 23 EBs with nearly identical primary and secondary eclipses separated by 50% in orbital phase. We concluded the majority of these systems have half their listed orbital periods, and therefore exhibit only one eclipse per orbit. If we were to fit physical models to these systems assuming the listed orbital periods, they would all have 0.84 and 0.08 (blue diamonds within red region of Fig. 12). A concentration of 23 EBs in this small corner of the parameter space is highly unlikely considering the density of systems in the surrounding phase space is substantially smaller. We expect only 3 - 5 of the 23 EBs to be twins in nearly circular orbits with the listed orbital periods. Given only the photometric data, however, we cannot easily determine which of the systems are truly twins with small eccentricities and which have half the listed orbital periods. In §2, we simply excluded all 23 EBs with ambiguous periods, and we account for the incidental removal of the 3 - 5 genuine systems in §5. We emphasize that most of the 23 EBs with ambiguous orbital periods have half the listed values, and therefore it was appropriate to remove these systems.

In Fig. 13, we compare the measured eccentricities to the arguments of periastron for the 130 well-defined EBs. Assuming random orientations, the periastron angle should be uniformly distributed across 0 360. However, the observed systems are not evenly concentrated across all and . We notice two observational biases in the data, both of which are due to geometrical selection effects.

First, for modest to large eccentricities 0.4, the EBs cluster near = 0 and = 180. In fact, the two systems with 0.7 - 0.8 have 0. For EBs with 90 or 270, one of the eclipses would occur at periastron while the other at apastron. The eclipse at periaston would be quite narrow according to Kepler’s second law, and may be too narrow to be accurately measured given the cadence of the OGLE-III data (see item B in §2). If the inclination is not sufficiently close to edge-on, e.g. 87, then the eclipse at apastron may be too shallow to be accurately measured (again, see item B in §2). If the inclination was even smaller, e.g. 85, the projected separation at apastron could be large enough so there would be no secondary eclipse. These systems would exhibit only one eclipse per orbit such as those presented in item A of §2. Considering the above, it is extremely difficult to observe and measure highly eccentric EBs with eclipses that occur near periastron and apastron. As the eccentricity increases, well-defined EBs are only detected as the argument of periastron approaches 0 or 180.

Second, there is an overabundance of EBs with 90 relative to those with 270. Quantitatively, there are 90 EBs with 0 180 and only 40 EBs with 180 360. These two values are discrepant at the 4.4 level according to Poisson statistics. This observational bias is due to our definition of the primary eclipse minimum , which determines the reference frame for . Recall the primary eclipse at must be deeper than the secondary eclipse. If 0.2, 89, and the primary is eclipsed closer to apastron, then the eclipse of the most massive luminous component may actually coincide with the secondary eclipse . Indeed, we found 18 EBs in such a configuration whereby is eclipsed at and is eclipsed at (see §3 and Table 2). Sixteen of these 18 EBs have 0 180. If we were to define according to instead of in terms of , then 74 EBs would have 0 180 and 58 EBs would have 180 360. These two values are now consistent with each other, i.e. they only differ at the 1.4 significance level.

As indicated above and discussed in §2, we suspect the majority of the 48 EBs we removed in items A and B of §2 have 0.4 and either 20 160 or 200 340. We test this hypothesis using the statistics of the measured systems in our well-defined sample. Of the 53 EBs with 0.4, 22 have 20, 160 200, or 340. If these 22 systems are complete across the specified intervals of , which total 80, and if the intrinsic distribution of periastron angles is uniform, then we expect 22 360/80 = 99 EBs with 0.4. We detected only 53 EBs with 0.4, implying 99 53 = 46 EBs did not satisfy our selection criteria. These 46 EBs most likley have secondary eclipses that are too narrow, too shallow, or completely absent. This prediction of 46 missing EBs nearly matches the 48 EBs we removed in items A and B of §2. This consistency further demonstrates that geometrical selection effects are understood in our sample and the removal of EBs in §2 were well-motivated.