Orbital and Stellar Parameters for 2M06464003+0109157: A Double-lined Eclipsing Binary of Spotted, Sub-solar Twins

The APOGEE spectrograph (Wilson et al. 2010, 2012; Majewski et al. 2017) covers the wavelength range of the H band (1.51–1.7 μm), with a spectral resolution of R ∼ 22,500. 2M0646 has been observed by APOGEE over the course of 6 epochs, from HJD of 2456312 to 2456673, with distinct spectral reductions in the 14th and 16th data releases (DR14 and DR16) of the Sloan Digital Sky Survey (Blanton et al. 2017; Abolfathi et al. 2018; Ahumada et al. 2020). Olney et al. (2020) report parameters for 2M0646 of Teff = 3870 ± 26 K and $\mathrm{log}g\,=4.416\pm 0.027$ dex, consistent with a low mass main sequence dwarf. Using the DR14 reduction of the data, 2M0646 is reported to have a metallicity of [Fe/H] = −0.302 ± 0.016 dex; however, it drops to −0.6 using the DR16 reduction. We should note, that all of these pipelines are assuming a single star. It is difficult to estimate how an SB2 should affect the stellar parameters, though Olney et al. (2020) note that they do not appear to systematically affect T_eff or $\mathrm{log}g$ distribution compared to single stars in a single age populations. Thus, most likely, they represent the parameters of the flux weighted mean of the two stars, which should usually gravitate to that of the primary.

In examining the cross-correlation function (CCF) of the system's APOGEE spectra (Nidever et al. 2015), we found that there was a signature of two strong components, indicating that this system is a double lined spectroscopic binary. In the DR14 data, components can be resolved in 4 out of 6 epochs, with the other two epochs unresolved due to the closeness of the velocities of the two stars. In the DR16 reduction, both components are resolved in all 6 epochs. The velocities of these components were extracted by fitting a Gaussian profile to the CCF using GaussPy (Lindner et al. 2015). These velocities extracted from each reduction are presented in Table 1; when resolved in both data releases, individual component velocities are typically consistent to within 1σ, though the velocities measured at MJD = 56316.7579 do differ by ∼2σ in each component. At the two epochs where the two components are not resolved in the DR14 CCFs, we report the velocity measured from the system's unresolved CCF in the v₁ column, but note this quantity may better represent the system's bulk heliocentric velocity than the velocity of any individual component at this epoch.

We report velocities measured from both reductions, and we use DR16 reduction for the orbital fit (discussed in Section 3.2).

2M0646 was identified as a P = 1.06493 days eclipsing binary (EB) by Pawlak et al. (2019) using light curves from the ASAS-SN survey (Kochanek et al. 2017). ASAS-SN conducts nightly observations from multiple observing stations around the world. Each station contains four 14 cm aperture telephoto lenses, each equipped with a CCD camera that provides a field of view of 4.5 deg². ASAS-SN has observed 2M06464 over 1878 epochs from 2014 to 2020 December at the time of analysis. Between 2014 and 2018, all epochs were observed in V-band. In 2018, ASAS-SN brought online their first g-band cameras and thus the observations were able to continue in both g-band and V-band. Even with the two different bands, the light curve was still sparsely populated.

These data have been since superseded by the light curves from Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015), which offer better cadencing and precision in the flux. TESS began its survey in 2018 July. It is equipped with four CCD cameras that provide a simultaneous coverage of a 96 × 24 deg² area on the sky. Each sector of the sky is observed for a period of ∼27.4 days with a 30 minute cadence. 2M0646 is located in sector 6, producing a TESS light curve that spans the period from 2018 December 15th to 2019 January 7th.

We used the eleanor pipeline (Feinstein et al. 2019) to download and extract the light curve of 2M0646 from the TESS full frame images archived in STScI's MAST database. This pipeline provides an option for several different calibrations to detrend the light curves and reduce the systematics of the instrument. We adopt the light curve produced by the point-spread function modeling, which performs well for brighter targets such as 2M0646.

The light curve is extracted with the photometry in the units of flux. However, the light curve fitting routine (Section 3.1) expects light curves to be magnitudes. We convert the flux to the T magnitude, ensuring that the average flux of the system is consistent with the magnitude reported in the TESS Input Catalog (T = 11.442 mag, Stassun et al. 2019). The resulting light curve, phase-folded on a 1.0649 days period, is shown in Figure 1, which shows a significant improvement in precision and coverage compared to the ASAS-SN data in Figure 2. In addition to better sampling of the eclipse themselves, the TESS light curve also reveals the presence of out-of-eclipse variations likely due to starspots.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

ellc (Maxted 2016) is a detached binary star modeler that is designed to allow a Markov Chain Monte Carlo (MCMC) sampling mode. We used ellc, as driven by the emcee MCMC package Foreman-Mackey et al. (2013) to fit 2M0646's light curve. Before performing a joint fit of the ASAS-SN, TESS, and APOGEE observations of this system (see next section), we conducted a parameter study to identify the sensitivity of the model results to the number and location of the starspots adopted in the fitting procedure. As the spot modulation is more precisely measured, and less temporally blurred, in the TESS light curve, we focused this parameter study on fitting that light curve in isolation.

For each trial in this parameter study, we initialized 200 MCMC walkers and computed chains for each walker with 20,000 steps. Starting with the broad distribution of parameters that were largely unconstrained, we allowed the model to converge and examined the outputs. Based on these solutions, we put increasingly restrictive priors on the model parameters to exclude unphysical solutions (e.g., negative periods, face-on inclination, unrealistic ratio of the two radii, etc.), without cutting into the physically plausible parameter space to which most of the chains had converged.

This process identified a core set of well motivated priors for the system's stellar and orbital parameters, which we used as consistent initial conditions while varying the number and location of spots in the model. Specifically, we set a Gaussian prior on the period, centered at 1.06 days with a width of 0.025 days. The stellar radii are expressed by ellc as R/a, where R is the physical radius of a given star, and a is the semimajor axis of the system, thus the radius can be inferred from the width of the system's eclipses in phase space. We set a Gaussian prior on R/a centered on R/a ∼ 0.13 with a width of ±0.06. The inclination of the system was set at 83° ± 15°, and the surface brightness ratio as 0.95 ± 0.25. Eccentricity was set as zero. The time of passage of the eclipse T₀ was restricted with a uniform prior from 0 to 1.1 (to prevent the solution wrapping around the phase), and was set centered at 0.75 days.

Using these priors, we then attempted to fit the TESS light curve with models assuming no star spots, one spot (on the system's primary), one spot per star, and two spots on the primary component. The best fitting model from each of these spot configurations is shown in Figure 1, and system and spot parameters are given in Tables 2 and 3, respectively. As Figure 1 shows, the best fit from the no-spot model is unable to reproduce the system's out-of-eclipse variations. Adding a single spot to the model significantly improves the fit, but still leaves systematic residuals in phase.

Ultimately, we find two spots are necessary to fully fit the out-of-eclipse light curve. In the two spots, one-per-component model, the spots were initialized with a size of 40° ± 20°, a longitude (l) of 180° ± 180°, a latitude (b) of 30° ± 30°, and a brightness factor of 0.7 ± 0.3. At the final initialization of the system, only one uniform prior was placed on l from 0° to 360° to keep it physically plausible. However, the model is not sensitive to the distribution of spots—the fit is comparable if both spots are located on the same stars versus on one star each. Similarly, the model is not sensitive to discriminating between dark spots and bright faculae. The resulting parameters across different models are shown in Table 2. The preferred spot sizes (∼35°) and contrasts (∼70%) are consistent across both dual-spot configurations.

Although the average l and b show a wide distribution in Table 3, this average is a superposition of various discrete locations that are equally favored. If faculae are used, they tend to be positioned facing toward the other star in the binary, whereas darker spots tend to be located on the opposite side. We interpret this as two stars that are tidally locked irradiating each other's photosphere, and consider the model where both spots are equally distributed across the stars as most likely.

The addition of the star spots to the solution does shrink the radii of the individual stars and increase the surface brightness ratio compared to the spotless model. This is due to ellc suppressing the average brightness of the spotted primary. Both single and double spot models produce a comparable effect. The spotted models are more consistent with the near-equal mass ratio system (see Section 3.2). Although both solutions that have two star spots fit the LC well, whether both of the spots are on the primary, or there is one on each of the stars, there is an 6% radii difference between the two which represents the systematic uncertainty in the solution. We chose, in our final fit, to adopt the "one-spot-per-component" model, as this seemed most consistent with the underlying symmetry in the system's parameters/mass ratio.

Initially, we explored extracting stellar and orbital parameters by using ellc to fit the system's light curve and The Joker (Price-Whelan et al. 2017, 2018) to fit the system's radial velocity measurements. Given the disjoint epochs of the TESS light curve and the APOGEE radial velocities, and the longer but sparser ASAS-SN monitoring, however, fitting the light curves and radial velocities separately led to non-negligible uncertainties in the system's period and T₀ (∼0.02 days).

To best constrain the system's orbital period and ephemeris, we utilized ellc's light curve and radial velocity modeling capabilities to jointly fit the system's TESS and ASAS-SN (V- and g-band) light curves, as well as the APOGEE radial velocities. In this joint fit, we utilize the ASAS-SN light curves mainly for the constraints their in-eclipse data points provide for the system's period and ephemeris; as 2M0646's starspots may evolve substantially over the baseline of the ASAS-SN monitoring, we artificially down-weight the ASAS-SN out-of-eclipse data points (which we define as any data point with phase between 0.05–0.45, and 0.55–0.95) by inflating their errors by a factor of 100. We then use emcee to identify 2M0646's maximum likelihood stellar and orbital parameters by using ellc to compute self-consistent light curve and radial velocity solutions for a given set of model parameters, and then compute a joint chi square value from the residuals of the fits to the TESS, ASAS-SN and APOGEE data. Given the relatively small number of radial velocity data points, and their orthogonal information content relative to the light curve, we apply a 20× multiple to the APOGEE residuals in computing the joint chi square value. In computing the light curve models, we adopt limb darkening coefficients computed by Claret & Bloemen (2011) for the V- and g-bands (for the ASAS-SN light curves), and for a white-light/Kepler bandpass (for the TESS light curve). In this joint fit, we applied priors sparingly, primarily to limit the solution to physically plausible solutions (i.e., restricting eccentricity to 0 < e < 1, the mass ratio to 0 < q < 1, etc.). We initialized the chain with Gaussian distributions around parameters which were found to well fit the light curve in the parameter study described earlier. Specifically, we initialized the walkers using the mean values and standard deviations listed in Table 4.

The remaining parameters that are customizable by ellc (e.g., third light, irradiation, apsidal motion rate, asynchronous rotation factor, surface gravity, [M/H], reflection coefficient, external noise, ${T}_{\exp }$) were not modified from their default values, as appropriate for a detached system of solar-metallicity main sequence dwarfs, or set to zero in the initial run.

We initialize a model with 15 walkers, each of which takes 5000 steps. We used the MCMC chain to estimate parameters for r₁, r₂, surface brightness ratio, inclination, T₀, and period. Before measuring parameters from the chain, we reject the first 10% of the chain to ensure that it has fully burned in (i.e., become relatively independent of where it was initialized). We then calculate the median and standard deviation of each parameter in the remainder of the chain, and adopt these values to characterize the system.

To investigate the effect of the two RV pairs of points near the eclipses (using DR16 data) on the derived orbit, we repeat the joint ellc fit and exclude them. The two RV fits are shown in Figure 3, and have negligible difference, either visually in examining the fitted orbit, or in the derived parameters.

Zoom In Zoom Out Reset image size

Best fitting parameters for the model containing all the data are given in Table 5. The derived masses and radii are show in Figure 4. The full corner plots of star+orbit and spot parameters are shown in Figures 5 and 6 respectively. To convert the radii from units of semimajor axis to solar radii, we used a variation of Kepler's 3rd Law to calculate the semimajor axis of the system. We then multiply the resultant semimajor axis by the relative radii (r/a) produced by ellc. Performing this conversion and propagating the errors in r/a, P, i, and v₁ + v₂, we find R₁ = 0.66 ± 0.05 R_☉ and R₂ = 0.57 ± 0.06 R_☉.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size