Joint Analysis of Multicolor Photometry: A New Approach to Constrain the Nature of Multiple-Star Systems Hosting Exoplanet Candidates

Multicolor transit photometry is a technique to observe planetary transits in multiple wavelength passbands (Tingley, 2004; Colón et al., 2012). If another source contaminates in a photometric aperture, the observed light curve is diluted. Due to the stellar color, the dilution ratio depends on each observing passband. Hence, we can obtain the information about contaminants from multicolor transit photometry. This technique is known as one of the planet-validation methods (e.g., Parviainen et al., 2019). The model of normalized light curve in a specific passband ($\lambda$) for a multiple-star system ($F_{\lambda}$) is expressed as

where $f_{\lambda}$ is the undiluted transit light curve in a specific passband whose analytic expression can be found, e.g., in Ohta et al. (2009), and $d_{\lambda}$ is a dilution factor that is the flux ratio of the central star to the total system including the companion(s). Throughout this paper, we assume physically associated and unresolved targets (e.g., eclipsing binary, hierarchical triple binary or planet in a binary system). In other words, we do not assume the contaminating light from background stars. This is because, except for very rare cases of “chance alignment,” those background/foreground sources generally have a good spatial separation (e.g., $>0\farcs 5$) from the target star, which can be easily identified by e.g., adaptive-optics (AO) or speckle observations. Therefore, we can derive $d_{\lambda}$ using the following equation.

where ${\rm Mag}_{*,\lambda}$ denotes the absolute magnitude in a given passband $\lambda$, and $i$ is the index for each star and “cen” means central star, which is the central star of the transiting/eclipsing system. In the calculation of the denominator, “cen” is included in the index $i$. To derive $d_{\lambda}$, we use the MESA Isochrones and Stellar Tracks²²2http://mesa.sourceforge.net/ (MIST: Dotter, 2016; Choi et al., 2016) for the fiducial case, which covers a wide range of stellar effective temperature and metallicity for the synthetic absolute magnitude $\rm Mag_{\lambda}$. Later, we will also use other isochrone models to discuss the isochrone dependence of the fitting results.

The shape of an observed SED also depends on stellar multiplicity. We aim to accurately estimate the properties of the systems showing transit-like signals by adding SED information to multicolor transit light curves. Here, we show our model of apparent magnitudes of the targets. The apparent magnitude in a specific passband is described as

where $\mu$ and $A_{\lambda}$ are the distance and the line-of-sight extinction modulus, respectively. Since we assume unresolved multiple-star systems, ${\rm Mag}_{\lambda}$ is obtained as total brightness of the stars in the system:

where $d_{\lambda}$ is same variable in Equation (2). We use the parallax from $Gaia$ DR2³³3https://gea.esac.esa.int/archive/ (Gaia Collaboration et al., 2016; Andrae et al., 2018; Gaia Collaboration et al., 2018) for the distance modulus. We multiply the color excess ($E(B-V)$) by empirical extinction vectors presented in Yuan et al. (2013) to derive $A_{\lambda}$, where $E(B-V)$ is calculated with ${\rm Mag}_{B}$ and ${\rm Mag}_{V}$, and the observational quantities ${\rm mag}_{B}$ and ${\rm mag}_{V}$. Here, we set the models of ${\rm mag}_{B}$ and ${\rm mag}_{V}$ to the observed values to derive $E(B-V)$. We also apply the MIST isochrones for synthetic absolute magnitudes. The apparent magnitudes in the $B$, $V$ (Bessell, 1990), $g$, $r$, $i$ (Sloan Digital Sky Survey ; York et al., 2000), $J$,$H$,$K$ (Two Micro All Sky Survey ; Skrutskie et al., 2006), W1, W2, W3 and W4 (Wide-field Infrared Survey Explorer ; Wright et al., 2010) bands are collected from catalogues on the VizieR⁴⁴4https://vizier.u-strasbg.fr/viz-bin/VizieR/ website. Finally, we calculate the synthetic apparent magnitudes in the passbands except $W3$ and $W4$, since Yuan et al. (2013) do not provide the extinction coefficients for those bands due to their low sensitivities.

Our aim is to constrain the nature of the unresolved candidates “quantitatively” by combining multicolor transit/eclipse photometry, SED, and isochrone fitting. Specifically, we assume the following four scenarios for the systems hosting transiting-planet candidates:

1. 1.

   planet (or dark object) transiting a single star, [PS];

2. 2.

   planet (or dark object) transiting a star in a binary system, [PB];

3. 3.

   eclipsing binary, [EB];

4. 4.

   triple-star system including an eclipsing binary, [ET].

The first two models are computed by setting the flux of the eclipsing object to 0. Additionally, all models have no assumption on the hierarchical structure for the components (objects) within the system; This means that the second model includes both p-type and s-type orbits and the latter two models include “secondary eclipse” for the identified (primary) flux decrease. In this paper, we do not consider systems with four or more stars, because such systems are much less common (e.g., Raghavan et al., 2010) than the above cases.

We model the observed light curves ($F_{\rm\lambda}$) and apparent magnitudes ($\rm mag_{\rm\lambda}$) for each scenario using the values of $\rm Mag_{\rm\lambda}$ for each stellar mass on the isochrones. For each scenario in which two or three objects are involved, we simultaneously fit the observed multicolor light curves and magnitudes. We take a Bayesian approach with the Markov chain Monte Carlo (MCMC) method and estimate the most likely stellar masses on the isochrones and the orbital parameters in $f_{\rm\lambda}$ based on marginalized posterior distributions. In this process, we employ the variable term in the log-likelihood function $\ln{L_{\rm var}}$ for maximization as follows:

where $\lambda_{F}$ and $\lambda_{m}$ indicate the passbands used for the light curve and the apparent magnitude, respectively, the subscripts “${\rm mod}$” and “${\rm obs}$” indicate the modeled and observed data, respectively, $\sigma_{*}$ is the uncertainty for each data point, and $F$ and $\rm mag$ are the models described above. We use a quadratic function for the flux baseline in $F$. In addition, we define $\beta=b\times\frac{R_{\rm cen}}{(R_{\rm cen}+R_{\rm tra})}$, where $R$ is the radius of each object and the subscripts “cen” and “tra” indicate the orbital central star and the transiting/eclipse object, respectively, as an alternative to the impact parameter $b$ in order to account for the case that the transiting object is larger than the central star. In cases where the transiting object is a star (i.e., EB or ET scenarios), we can derive all $R$ values from the effective temperatures and the luminosities on the isochrones by the Stefan-Boltzmann’s law. On the other hand, for the case of planetary system (i.e., PS or PB scenarios), only $R_{\rm tra}$ is independent of the isochrones.

We use the Python package emcee (Foreman-Mackey et al., 2013) ensemble sampler for the MCMC in the fitting process and evaluation of uncertainties in the derived quantities. The fitting parameters in all our models are the time of transit center $t_{c}$, scaled semi-major axis $a/R_{\rm cen}$, reduced transit/eclipse impact parameter $\beta$, quadratic baseline coefficients for the passband $c_{\lambda,0}$, $c_{\lambda,1}$ and $c_{\lambda,2}$, and central stellar mass $M_{\rm cen}$. Additional fitting parameters are the companion star’s mass $M_{\rm com}$ for the PB and ET scenarios, radius of the transiting object $R_{\rm tra}$ for the PS and PB scenarios, and eclipsing star’s mass $M_{\rm tra}$ in the EB and ET scenarios, respectively. We apply uniform distributions as uninformative priors for these fitting parameters. Initial positions of the walkers are set to uniform distributions in their domain ranges, to search the parameter space widely. There are uncertainties about $\rm mag_{mod,{\lambda}}$ due to the parallax and the extinction vector. We thus randomly vary these parameters with Gaussian distributions in each step of the chain. The widths of Gaussians follow the nominal errors in Gaia parallax measurements and in Yuan et al. (2013) for the extinction. After deriving the quantitative physical properties by the MCMC fittings for all scenarios, we compare and evaluate each reduced $\chi^{2}$ (: $\chi^{2}_{\rm red}$) and Bayesian information criteria (BIC ; Schwarz, 1978), and find the most plausible scenario(s).

In order to demonstrate our method and software, we selected several candidate transiting systems identified by the K2 mission, which is an extended mission of the Kepler prime mission (Borucki et al., 2010; Howell et al., 2014). We successfully observed six systems with multicolor transit photometry using the Multicolor Simultaneous Camera for Studying Atmospheres of Transiting Exoplanets (MuSCAT) mounted on the Okayama 188cm telescope (Narita et al., 2015) and/or Simultaneous-color Infra-Red Imager for Unbiased Survey (SIRIUS; Nagashima et al., 1999) atop the Infrared Survey Facility (IRSF 1.4 m) telescope at the South Africa Astronomical Observatory. EPIC ID and stellar parameters of the six targets are summarized in Table 1 (described in Section 3.2). The parallax values were adopted from the Gaia archive for all targets. The apparent magnitudes from VizieR are listed in Table 2. The orbital period, transit depth, and radius ratio of each candidate in the literature are listed in Table 3.

We downloaded the reduced K2 light curves provided by Vanderburg & Johnson (2014). We normalized the light curves by fitting it by a fifth-order polynomial function. We folded the light curves with the orbital periods in Table 3 and removed flux outliers by a $4$ - $\sigma$ clipping for each 0.05 day-long segment. Finally, we folded the light curves with the orbital periods in Table 3 and binned into 0.005 day. Since the six targets listed in Table 1 cover a wide range of the parameter space in terms of the stellar mass, transit depth, and period, etc, those targets become a good set of samples to test and demonstrate our methodology and software to identify the true nature of the systems. Below we give a brief summary of each target before applying our analysis tool.

To better characterize our targets, we determined the stellar properties of the targets using the archived data available. For EPIC 206036749, EPIC 210513446, EPIC 211800191 and EPIC 220621087, we downloaded the high-resolution spectra observed by Keck/HIRES available at ExoFOP-K2. For EPIC 206500801, a spectrum observed with HARPS-N (Cosentino et al., 2012) was downloaded from IA2/TNG archived website⁶⁶6http://archives.ia2.inaf.it/tng/. We made use of SpecMatch-Emp (Yee et al., 2017) to derive stellar properties ($T_{\rm eff}$, $R_{*}$, $M_{*}$ and [Fe/H]) of the targets from high-resolution spectra. These high-resolution spectra were analyzed following Hirano et al. (2018). For EPIC 220696233, for which no high-resolution spectrum was available, we adopted the values of effective temperature, metallicity, and $\log{g}$ listed on ExoFOP-K2 (Huber et al., 2016), but we derived its stellar radius and mass using the empirical relations for low-mass stars (Mann et al., 2015) based on the above parameters and Gaia parallax. The results of these analyzes are summarized in Table 1. Those parameters are helpful to understand the properties of the primary (brightest) star for each system, but in the following analyses, we only use the derived metallicities in generating the isochrones during the fitting process, as there are large uncertainties in spectroscopically derived parameters when the system includes multiple stars. The other stellar parameters in Table 1 are used for comparisons of the derived results.

The most likely scenario for EPIC 206036749 is ET. There are good agreements between the observed quantities and their models for both light curves and magnitudes (Figures 1 and 2). In this scenario, the system consists of an eclipsing star ($0.60\,M_{\odot}$), the central star ($0.48\,M_{\odot}$), and the companion star whose mass is $1.04\,M_{\odot}$. The eclipsing star is estimated to be larger than the central star, suggesting that the reported transit-like events are actually secondary eclipses. The second possible scenario, although much less likely according to the large $\Delta\mathrm{BIC}$, is PB consisting of a star-sized ($0.19\,R_{\odot}$) dark object, a central star with $0.47\,M_{\odot}$, and a companion star with $1.10\,M_{\odot}$. The other two scenarios (PS and EB) are very unlikely, since both $\chi^{2}_{\rm red}$ and BIC values are much worse than the above two scenarios.

In both of the results for DSEP and PARSEC, the ET scenario is selected as the minimum BIC model as shown in Table 5. The derived properties are consistent within $\approx 1\sigma$ among the three different isochrone models. The uncertainties for the eclipsing star’s mass and the companion star’s mass in DSEP are large due to the bimodal posterior distributions of the fitting parameters. The highest peak is compatible with the results for the other two isochrones, and the second peak is located around $1.0\,M_{\odot}$ - $0.4\,M_{\odot}$ - $0.6\,M_{\odot}$ for $M_{\rm tra}$ - $M_{\rm cen}$ - $M_{\rm com}$, respectively. This disagreement may be due to the fact that the combination of these three stellar masses degenerates with respect to the SED, although the reason why the bimodal distribution was only seen for DSEP is unclear.

Figure 5 indicates that the secondary eclipse (phase near 0.5) is much shallower than the primary transit in the K2 light curve. This suggests either that the orbital inclination and eccentricity are large (grazing eclipses) or that the transiting/eclipsing object is much cooler than the eclipsed star. The former case is possible in the ET scenario, but the possibility of a large eccentricity has to be low due to the short orbital period (1.13 days) and the observed timing of the secondary eclipse near phase $\approx 0.5$. On the other hand, the latter possibility is compatible with the PB scenario in Table 5; in the PB scenario, the transiting object is assumed to be “dark”, meaning that it does not affect the color-dependence of the transit light curves and SED, but this scenario also allows for a very late-M dwarf or a brown dwarf as the transiting object, which would give negligible impacts on the observable quantities. The detection of the secondary eclipse is indicative of the possibility that the transiting/eclipsing object is actually a very low-mass star rather than a giant planet.

The folded light curve also exhibits a phase-shifted sinusoidal modulation. This modulation is unlikely to be caused by ellipsoidal variations, in which case it would show two peaks in one orbit (Morris & Naftilan, 1993; Mazeh & Faigler, 2010). The flux peak/bottom phase is not consistent with the Doppler boosting, either. It may be ascribed to interactions of close-in binary stars such as the O’Connell effect (Wilsey & Beaky, 2009). As we noted in Section 3.1.1, the contamination star is not resolved within $0\farcs 5$, corresponding to $\sim 200$ au in the AO image. This is a tight constraint for a triple-star systems. In Crossfield et al. (2016), EPIC 206036749 is reported to have a Neptune-sized planetary candidate, but our analysis implies a compact multiple-star system (i.e., ET or PB).

The most likely scenario for EPIC 206500801 is PS, consisting of a $0.22\,R_{\odot}$ transiting/eclipsing object and the central star with $1.33\,M_{\odot}$. The $\chi_{\rm red}^{2}$ and BIC in the other scenarios are much worse than those for PS. The minimum BIC scenario with DSEP is also PS, in which the derived parameters in Table 5 are consistent within $\approx 1\sigma$ with those of MIST. On the other hand, the analysis with PARSEC favors PB, consisting of a $0.18\,R_{\odot}$ transiting object, central star with $1.09\,M_{\odot}$, and companion with $0.93\,M_{\odot}$. This disagreement in the scenario selections may be due to the difference in SED models for different isochrones.

When we naively interpret the observed transit/eclipse depths for EPIC 206500801, an EB scenario with e.g, $\approx 0.2\,M_{\odot}$ - $\approx 1.3\,M_{\odot}$ for $M_{\rm tra}$ - $M_{\rm cen}$ is also expected as a plausible solution, but this possibility was ruled out in the above analysis. To further investigate this, we reanalyzed the same data set imposing a prior for the eclipsing star being less massive than $0.2\,M_{\odot}$. In this case, the solution for the EB scenario was found to be $0.16\,M_{\odot}$ - $1.25\,M_{\odot}$ for $M_{\rm tra}$ - $M_{\rm cen}$, where the $\chi^{2}_{\rm red}$ and BIC values were 1.43 and 3078.53, respectively, which are also listed in Table 4 as $\rm EB^{*}$ in addition to the other scenarios. While the BIC value is the lowest than those for all the scenarios without a prior, the goodness of fit in $\chi_{2}^{2}$ ($\chi^{2}$ for the apparent magnitudes) is the worst among all scenarios. The disagreement between the fitting results with and without the prior suggests that there are multiple local mimima in $\chi^{2}$, leading to a degeneracy among the fitting parameters. The MCMC analysis without a prior explored a wider range of the parameter space, but this degeneracy of the fitting parameters prohibited a robust estimation.

Consequently, we were unable to conclude if EPIC 206500801 is a PS or EB. This is because the flux from a low-mass star less massive than $0.2\,M_{\odot}$ is negligible compared to the solar-type stars in the total SED, and we cannot distinguish between a giant transiting planet and a very low-mass star. Fortunately, it is relatively straightforward to identify such a close-in binary, by e.g., making a few radial-velocity measurements. Note that binary systems whose mass ratio is lower than 0.2 and orbital period is shorter than 10 days are not so common (see Figure 17 in Raghavan et al., 2010). In Figure 5, the phase curve looks sinusoidal, but its exact mechanism is not known as in the case of EPIC 206036749. A weak secondary eclipse-like signal ($\sim 0.1\%$) can be seen with the blue arrow in Figure 5, making EB a more likely scenario for this system. Further follow-up observations would be required to make a definitive conclusion.

The best-fit scenarios for EPIC 210513446 are EB and ET, which are equally likely based on the BIC values in our analysis. In the EB scenario, masses of the eclipsing star and the central star are $0.82\,M_{\odot}$ and $0.54\,M_{\odot}$, respectively, and the observed eclipse is actually a secondary eclipse. In the ET scenario, the system consists of an eclipsing star with $0.56\,M_{\odot}$, the central star with $0.48\,M_{\odot}$, and the companion star with $0.76\,M_{\odot}$. For both DSEP and PARSEC models, the ET scenario was selected as the minimum BIC model as in Table 5. The derived physical parameters are consistent within $\approx 1\sigma$ with those for MIST.

As we noted in the Section 3.1.3, a contaminant source was resolved in $\approx 0^{\prime\prime}.5$, corresponding to $\sim$ 160 au with $\Delta{\rm mag_{K}}\approx 1$ in the AO image. This magnitude difference is almost consistent with the difference between a $0.6\,M_{\odot}$ - $0.5\,M_{\odot}$ binary and a $0.8\,M_{\odot}$ companion star in our ET scenario. If the system includes similar-mass eclipsing binary stars, we would observe secondary eclipses whose depths are similar to the primary, but the secondary eclipse as well as the photometric variations synchronized with the eclipses was not observed with the 1.15-day period in Figure 5 (black dots). We did not find phase variations in the light curve folded with the double period (2.30 days; grey dots in Figure 5), either. On the other hand, we detected a large periodic variation at 7.54 days based on the generalized Lomb-Scargle periodogram (Zechmeister & Kürster, 2009), which is shown with blue dots in Figure 5. Given the large amplitude of this flux modulation, this periodicity at 7.54 days most likely comes from the companion star, emitting the largest flux among the three stars in the system. Since the 2.30-day orbital period more consistently explains the observed data (i.e., the secondary eclipse), we conclude that EPIC 210513446 is likely to be an ET including an eclipsing binary ($0.6\,M_{\odot}$ - $0.5\,M_{\odot}$) with a period of 2.30 days.

The minimum BIC scenario is PB consisting of a $0.16\,R_{\odot}$ transiting object, the central star with $0.30\,M_{\odot}$ and the companion star with $1.04\,M_{\odot}$, but derived values of $\chi^{2}_{\rm red}$ and BIC are similar between PB and ET scenarios, and we cannot select the most likely scenario from this analysis. The ET scenario was selected as the minimum BIC model (Table 5) for both DSEP and PARSEC models. The derived values of the physical parameters are consistent within $\approx 2\sigma$ with those of MIST.

For the EB and ET scenarios, secondary eclipses would be seen in the K2 light curve, but were not detected. Thus, we suspected that the orbital period in Table 3 was incorrectly estimated. We folded the K2 light curve with the double period, and found the synchronized photometric phase curve and two eclipses as shown by grey dots in Figure 5. We checked the ephemeris and reanalyzed the K2 light curve using the same data set and conditions as above, but with the double period. As a result, we derived the ET scenario as the most plausible one with $0.50-0.39-1.01\,M_{\odot}$ for the MIST isochrone, which is almost the same as the original result in Table 4 within $1\,\sigma$. If the EB or PS scenarios were true, then the secondary eclipse would be much shallower than the primary transit, which was not observed. We thus conclude that ET is the most reasonable scenario for the EPIC 211800191 system. As we noted in Section 3.1.4, no apparent contaminant was resolved within $0\farcs 5$, corresponding to $\sim 200$ au in the AO image.

The minimum BIC scenario for MIST is ET, consisting of the eclipsing star with $0.24\,M_{\odot}$, the central star with $0.23\,M_{\odot}$, and the companion star with $0.33\,M_{\odot}$. The result for DSEP selects the ET scenario as the minimum BIC model as in Table 5, which is consistent within $\approx 1\sigma$ with that for MIST. On the other hand, the result with PARSEC favors PB scenario, consisting of a planet with a size of $0.02\,R_{\odot}$, the central star with $0.20\,M_{\odot}$, and the companion star with $0.44\,M_{\odot}$.

As we noted, EPIC 220621087 was already validated as a true planet host (K2-151: Hirano et al., 2018). In our assessment, however, EPIC 220621087 is a possible false positive, which consists of similar-sized M dwarfs as in Table 4. If the false-positive scenarios were true, then we would see the multiple features in the cross-correlation (line profile) of its high-resolution spectrum, since the masses of the two components are similar to each other in the ET scenario. In Hirano et al. (2018), a single-lined profile was shown for EPIC 220621087, ruling out the false positive scenario that our analysis software has returned. Similarly, a PB scenario can be ruled out by AO imaging and cross-correlation profile. In addition, phase-curve variations are not seen for both single and double periods in Figure 5, suggesting a low probability with the system being a binary.

The inconsistency of our analysis results made us consider that it is difficult to distinguish a true scenario from other possibilities, when the transit depths for all passbands are similar and relatively shallow ($\lesssim 0.1\%$) in comparison with the precision of the ground-based photometry. In such cases, the scenario with a single star and ones including multiple stars of similar masses produce similar observational signals, leading to a mis-classification in the output. The same is true in the case of the scenario with small transiting objects and the scenario hosting grazing object. Another possible reason for the inconsistency is the incompleteness of the isochrone models; the low-mass range ($<\,0.3\,M_{\odot}$) of the isochrones is known to show systematic deviations in theoretical values from the empirical ones (Mann et al., 2015), suggesting that the analysis for cool targets may have large systematic uncertainties in the derived parameters.

The most likely scenarios are PS and PB for the MIST isochrones. In the former case, the central star has a mass of $0.61\,M_{\odot}$, while for the latter case, the central star is a $0.54\,M_{\odot}$ M-dwarf and the companion star has $0.38\,M_{\odot}$. In both scenarios, the radius of the transiting planet is $\approx 0.06\,R_{\odot}$. The flat bottom of the transit curve in Figure 1 is also suggestive of the planetary nature of the candidate. Both analyses with DSEP and PARSEC resulted in PS as the minimum BIC scenario. The properties in the PS scenarios are consistent within $\approx 2\sigma$ among these analyses.

As we noted in Section 3.1.6, no contaminant is resolved within $0^{\prime\prime}.5$ by AO imaging, corresponding to $\approx 130$ au. In Figure 5, phase-curve variations are not seen at both single and double periods. In the PS scenario, the planet candidate is a super-Neptune-sized planet orbiting a $0.6\,M_{\odot}$ star (i.e., early M dwarf). In order to rule out the possibility of EB completely, high resolution spectroscopy would be required. Although this target is optically faint ($\sim 16$ mag in the V band), the detection of a $0.5\,M_{\odot}$ - $0.5\,M_{\odot}$ binary (i.e., the secondary feature in the cross-correlation profile) would be relatively easy by e.g., near-infrared high-resolution spectroscopy.

There are two types of systematic effects due to the isochrones in the fitting result as shown in Section 4. One pertains to the selection of the most plausible scenario with the minimum BIC value. This systematic effect solely depends on the goodness of fit, because the degree of freedom is the same (fixed) for the different isochrones. In the fitting processes, the light curve model (the first term in Equation 5) is relatively insensitive to the selected isochrone model, since only the dilution factor $d_{\lambda}$ depends on the color returned by the isochrone model, and the other orbital parameters are independent of isochrones. However, the SED model (the second term in Equation 5) heavily relies on the isochrone-dependent color difference, and therefore the plausible-model selection by BIC is very sensitive to the isochrones adopted in the analysis. The other impact of isochrone models on the fit result is the difference in the derived values for each scenario. In most of the cases presented in this work, the derived values are consistent within $\approx 2\sigma$ among the three isochrones. The bimodal posterior distribution shown by EPIC 206036749 is sometimes seen for ET scenarios, suggesting that there is a degeneracy in the expected combined contrast in cases of more than two stars.

We estimated the expected error propagation from the metallicity uncertainty, adopting EPIC 206036749 as a representative sample. In doing so, we used the MIST isochrones with [Fe/H] = 0.27 and 0.43 and assessed the variations in the derived parameters. The best-fit stellar masses in ET are 0.57 - 0.43 -1.06 $\,M_{\odot}$ and 0.60 - 0.46 - 1.06 $M_{\odot}$ for $M_{\rm tra}$ - $M_{\rm cen}$ - $M_{\rm com}$, respectively. Thus, the systematic error for the estimated masses arising from fixing the metallicity of the system is no more than $0.03\,M_{\odot}$, which is generally smaller than statistical errors in the derived masses (Table 4). We also investigated age-dependency by applying MIST isochrones with 5.0 Gyr, finding 0.58 - 0.44 - 1.02 $M_{\odot}$ for $M_{\rm tra}$ - $M_{\rm cen}$ - $M_{\rm com}$. The small deviations of those values from the fiducial result (Table 4) suggest that our analysis is insusceptible to relatively small uncertainties in the system age and metallicity, as long as the system is in the main sequence of the evolutionary track.

In following up candidate transiting systems, it is not always possible to conduct multicolor transit photometry from the optical (e.g., by MuSCAT) to near-infrared (e.g., by IRSF) wavelengths as we did for EPIC 206036749 and EPIC 220621087. In particular, for the ET scenario containing more than two stars, completely different fitting results may be derived, depending on the passbands used in the analysis.

To evaluate the dependence of the fitting result on the light-curve passbands used in the analysis, we performed an additional analysis for the ET scenario of EPIC 206036749, which shows significant transit-depth variations against wavelength. We estimated the properties for the two different subsets of light curves: K2 plus SIRIUS light curves, and K2 plus MuSCAT light curves. The result of this additional test is shown in Table 6. The analysis of K2 plus SIRIUS light curves resulted in bimodal posterior distributions as in Section 4.1, but the highest-peak position was consistent with the K2+MuSCAT result within $1\sigma$. With regard to the median values, the result for the K2 plus SIRIUS analysis is closer to the result in Table 4, in which all the light curves were used. We conclude that, while using a subset of multicolor light curves produces consistent results in most cases, limited light curve data sometimes result in a worse convergence in the posterior distributions of fitting parameters.