Revisiting Kepler Transiting Systems: Unvetting Planets and Constraining Relationships among Harmonics in Phase Curves

Space-based photometric missions widely use statistical validation tools for vetting transiting planetary candidates, particularly when other traditional methods of planet confirmation are unviable. In this paper, we refute the planetary nature of three previously validated planets -- Kepler-854~b, Kepler-840~b, and Kepler-699~b -- and possibly a fourth, Kepler-747~b, using updated stellar parameters from Gaia and phase-curve analysis. In all four cases, the inferred physical radii rule out their planetary nature given the stellar radiation the companions receive. For Kepler-854~b, the mass derived from the host star's ellipsoidal variation, which had not been part of the original vetting procedure, similarly points to a non-planetary value. To contextualize our understanding of the phase curve for stellar mass companions in particular and extend our understanding of high-order harmonics, we examine Kepler eclipsing binaries with periods between 1.5 and 10 days. Using a sample of 20 systems, we report a strong power-law relation between the second cosine harmonic of the phase-curve signal and the higher cosine harmonics, which supports the hypothesis that those signals arise from the tidal interaction between the binary components. We find that the ratio between the second and third-harmonic amplitudes is $2.24 \pm 0.48$, in good agreement with the expected value of 2.4 from the classical formalism for the ellipsoidal distortion.


INTRODUCTION
Space-based searches for transiting planets, including CoRoT (Auvergne et al. 2009), Kepler , K2 (Howell et al. 2014), and TESS (Ricker et al. 2014), have detected thousands of transiting planet candidates around faint host stars. The observational resources required to confirm or refute the planetary nature of each candidate are prohibitively large, especially in the case of the faintest targets. The necessity to vet the numerous candidates has led to the development of statistical validation techniques (e.g., Torres et al. 2011Torres et al. , 2015Torres et al. , 2017Morton & Johnson 2011;Morton 2012;Díaz et al. 2014; Morton et al. 2016). Such frameworks provide estimates of the probability that a candidate is a * NASA Postdoctoral Program Fellow false positive using a combination of stellar parameters and transit-derived properties. A candidate is declared as a validated planet if the false positive probability is below a certain threshold (usually 0.01 or 0.001, e.g. Crossfield et al. 2016). Since the Kepler era, these methods have been widely used for space-based photometric missions, including K2 (Montet et al. 2015;Crossfield et al. 2016;Mayo et al. 2018;Livingston et al. 2018;de Leon et al. 2021) and TESS (Giacalone et al. 2021).
Despite the popularity and broad applicability of statistical validation methods, there are cases where statistically validated planets have been subsequently shown not to be bonafide planets through more resourceintensive follow-up observations, such as radial velocity monitoring (e.g., Shporer et al. 2017). Similarly, revised stellar parameters, notably those provided by the Gaia mission (Gaia Collaboration et al. 2018, can re-veal misclassified planets as we will be demonstrating in this paper. Beyond radial velocity (RV) follow-up and revised stellar parameters, phase curves, which denote photometric modulations that are synchronous with the orbital period, provide a powerful tool for gaining insights into the mutual gravitational and radiative interactions between objects in binary systems. The three primary astrophysical phenomena driving these modulations are (1) ellipsoidal variation originating from the tidal interaction, (2) beaming (or Doppler boosting) due to relativistic Doppler effects, and (3) mutual illumination due to the mutual heating of the components' atmospheres. In the past, phase curves have been used to validate exoplanets, such as Kepler-41 b (Quintana et al. 2013) and Kepler-78 b (Sanchis-Ojeda et al. 2013). For shortperiod candidates (typically shorter than 5 days), the phase curve can provide an independent estimate of the companion mass and complement other planet characterization techniques without requiring further observations (e.g., Faigler & Mazeh 2011;Shporer 2017).
In this paper, we reassess the planetary nature of Kepler companions in light of updated stellar catalogs and then consider how phase-curve photometry can be used to mitigate the occurrence of misclassified companions. We find that three of them, and possibly a fourth, cannot be planets or brown dwarfs and are thus stellar objects that have been misclassified as planets. Then, we use a sample of Kepler eclipsing stellar binaries to show that the higher-order harmonic seen in the phase curve of one of the misclassified Kepler planets is likely to be a manifestation of the host star's tidal distortion. This paper is organized as follows. We describe the methods used in vetting the planetary nature of companions in Section 2. In Section 3, we apply these techniques to argue against the planetary nature of four previously validated Kepler planets. In Section 4, we extend our study of orbital phase curves to a sample of Kepler eclipsing stellar binaries to investigate the higher-harmonic phase-curve modulations. Finally, a summary of our main findings is given in Section 5.

REVISITING VALIDATED KEPLER PLANETS
Using the Kepler photometric data along with updated stellar parameters can help us identify some of the false positives hidden among the current planetary population. Better constraints on the distance from astrometric missions such as Gaia (Gaia Collaboration et al. 2018 provide better estimates of the physical radius of the stars, thereby leading to more accurate physical radii of the transiting companions. Meanwhile, phase curves, which typically are not used for vetting planetary candidates, enable us to infer the mass of the companion via the characteristic modulations that arise from the mutual gravitational interaction. Usually, the amplitude of the ellipsoidal variation in optical wavelengths does not exceed 200 ppm in a planet-star system (see Millholland & Laughlin 2017), allowing us to rule out some of the obvious false positives. Further analysis of the phase curve, especially the amplitudes at the first harmonic of the orbital frequency, can indicate if the transiting companion is self-luminous (i.e., a star) or a nonluminous planet.

Updated Stellar Parameters
The parameters of transiting companions, such as radius, are often derived from measured quantities that are defined relative to the host star's properties. Therefore, any bias in the stellar parameters will propagate directly as a bias in the companion's properties. We searched for substantial changes in the inferred radii of Kepler planets that are caused by revised stellar parameters. To that end, we first considered the differences in stellar radius of the confirmed Kepler planet hosts between the Kepler Input Catalog (KIC; original catalog) (Brown et al. 2011) and Version 8 of the TESS Input Catalog (TIC v8.2; revised catalog) (Stassun et al. 2019). As seen in the top-left panel of Figure 1, we found that the majority of host stars have comparable stellar radius values in the two catalogs. However, ∼3.5% of the stars show radius discrepancies of larger than a factor of two. These differences are primarily driven by the updated heliocentric distances provided by the Gaia DR2 parallax measurements (Gaia Collaboration et al. 2018) that were used to estimate the TIC stellar radii. The histogram also reveals an asymmetric distribution, with the deviation strongly skewed toward larger ratios (+0.25 versus −0.14). For most objects beyond 1,250 pc, the stellar radii have been revised to larger values owing to the previous underestimation of their distance (see the top-right panel of Figure 1).
For the next step of our assessment of Kepler planetary systems, we took the published planet-star radius ratio for each planet-hosting system and calculated the planet radius, using both KIC or TIC value for the host star's radius. The bottom-left panel of Figure 1 shows the comparison between the TIC-and KIC-based planet radii. Focusing on the upper half of the graph, four putative planets clearly stand out as having significantly larger TIC-based radii, and are therefore suspect as possible misclassified objects. For these systems, we performed independent fits to the host star's stellar energy distribution (SED) using the isochrone package (Morton 2015). In this analysis, we used publicly  Top left: histogram of the ratioR between the TIC-based and KIC-based stellar radius for confirmed Kepler planet-hosting stars. About 3.5% of hosts had their radius over-or underestimated by a factor larger than 2 in the KIC. Top right: TIC-to-KIC radius ratio as a function of the distance listed in the TIC. Cases with larger ratios typically correspond to targets that lie at heliocentric distances greater than 1,250 pc. Bottom left: a direct comparison of the TIC-based planetary radii versus the KIC-based planetary radii. The 1:1 line is shown in green for guidance. The four previously validated candidates with the largest TIC radius values are marked in red. Bottom right: plot of planet radius vs. insolation calculated assuming a circular orbit, which shows the radius inflation trend beyond an insolation flux of around 2 × 10 5 W/m 2 (Demory & Seager 2011). The data were obtained from the NASA Exoplanet Archive in November 2021. It is evident that four of our targets -Kepler-854 b, Kepler-840 b, Kepler-699 b, and Kepler-747 b -have discrepant radii, especially the latter two, given their insolation levels. This indicates that these companions are not planetary in nature. Kepler-447 b, NGTS-1 b, and TOI-1130 c, which lie near Kepler-747 b in this plot, are grazing planets with larger than typical radius uncertainties.
available broadband photometric measurements from the AllWISE source catalog (B, J, H, K s , W 1, W 2, W 3, W 4; Cutri et al. 2013), alongside the Gaia Early Data Release 3 (EDR3) parallax (Gaia Collaboration et al. 2021). As a separate check for consistency, we compared the density of the host star obtained from our SED analysis to the density derived from the Kepler transit fit assuming a circular orbit (described, for example, in Equation 30 of Winn 2010).

Phase-curve Model
Let Φ represent the orbital phase from the mid-transit time (i.e., inferior conjunction). We can mathematically express Φ as follows: Here, T 0 is the time of mid-transit, and P is the orbital period. In our phase-curve analysis, we considered up to five harmonics of the orbital phase, with both sine and cosine terms contribute to the phase curve flux (L PC ): Note that in this formulation all components are independent of each other. Many of these components can be directly tied to the various astrophysical phenomena that occur in binary systems. The first cosine harmonic (B 1 ) originates from the mutual illumination of the stars, whereas the ellipsoidal effect dominates the second (B 2 ) and higher (B 3 , B 4 ) cosine harmonics (Kopal 1959). The beaming effect itself manifests as the first sine term (A 1 ). While the higher-order sine amplitudes are not expected from theoretical modeling of binary-system phase curves, they have been observed for several systems in the past (e.g. Rappaport et al. 2015;Wong et al. 2020a), and therefore we include them for completeness. For binary systems in which both components have significant luminosity, it is important to account for the contributions from both stars. As the orbital phase of the secondary is shifted by 1/2 relative to the primary, all of the odd harmonics for both sine and cosine components (A 1 , A 3 , A 5 , B 1 , B 3 , and B 5 ) are subtractive, whereas the even harmonics (A 2 , A 4 , B 2 , and B 4 ) are additive. As a result, if the contribution from the secondary object dominates, the predicted signals calculated when assuming the primary star is brighter than the secondary may not only show a mismatch with the measurements but could even end up with a reversed sign altogether. In the following subsections, we describe the major phase-curve components.

Ellipsoidal Variation
The ellipsoidal variation arises from the mutual gravitational interaction between the two binary components. The primary physical processes that drive this modulation are the equilibrium tidal bulge raised on the primary star by the companion and the limb-and gravitydarkening of the distended atmosphere. The analytical formulation describing ellipsoidal variation in stellar atmospheres was first put forth by Kopal (1959), with updated mathematical expressions of the corresponding phase-curve amplitudes provided more recently in Morris (1985) and Morris & Naftilan (1993). Generally, the ellipsoidal induced flux variation (L E ) in the phase curve is expressed as a series of cosine terms: where B k are coefficients for the different ellipsoidal harmonics. The classical analytical formulation of ellipsoidal variation provides mathematical expressions for the first four harmonics (e.g., Morris 1985): In the above equations, q denotes the mass ratio between the binary components, i is the orbital inclination, and a/R 1 is the ratio between the orbital semimajor axis and the radius of the primary star. The terms represented by X and G contain the contributions from limband gravity-darkening, respectively. It has been standard practice in the literature to model ellipsoidal variation assuming linear limb-darkening models. In that case, the contribution of limb-and gravity-darkening can be written as: , G 2 = 1 + g, , G 3 = 2 + g, where u is the linear limb-darkening parameter, and g is the gravity-darkening parameter. For both linear limb-darkening and gravity-darkening parameter values, the values can be interpolated using the relevant tabulated parameters provided in Claret & Bloemen (2011). The precision of these coefficients are estimated through Monte Carlo sampling using the uncertainties on the reported stellar parameters. Combining the contributions from both components in the binary system, we arrive at the full mathematical expression of the leading-term amplitude of the ellipsoidal variation (i.e., at the second harmonic of the cosine); here the limb-and gravity-darkening terms are not included (they are set to one), for simplicity: Here, R 2 represents the radius of the secondary, and f 2 denotes the relative luminosity (i.e., flux ratio) of the secondary in the observed bandpass.
For the third-harmonic amplitude, which is typically the second-order term in the description of the photometric ellipsoidal variation, we can combine the two contributions analogously: Now, taking the ratio between these two amplitudes, we obtain The contribution of the secondary, as seen in the expression above, depends both on the radius ratio R 2 /R 1 and the flux ratio f 2 . Our ideal targets contain secondaries that are significantly smaller than the primary (at least by a factor of 2) and contribute negligible flux (in the optical) compared to that of the primary star. In this scenario, both the nominator and the denominator in the parenthesis in Equation 8 can be approximated as q 2 and thus cancel out, making the ratio |B 2 /B 3 | independent of q. Reintroducing the limb-and gravitydarkening terms of the primary star for systems where the secondary is small and faint, we arrive at the approximate expression We note that the limb-and gravity-darkening terms at the second and third harmonics can differ by more than a factor of 2, so careful modeling of their effects is crucial for studying the precise relationship between the ratio of B 2 and B 3 . Besides these terms, all of the parameters can be determined from the fit, and so the above equation allows us to test the classical formulation of Kopal (1959) without needing to know the mass ratio.

Beaming Effect
Beaming, otherwise known as Doppler boosting, is a relativistic effect that is caused by the line-of-sight acceleration of the two binary components around the common center of mass. For circular systems, this effect manifests itself as a photometric modulation at the first harmonic of the sine of the orbital phase and mirrors the shape of the radial velocity curve, albeit with the opposite sign (Loeb & Gaudi 2003): Here, M 1 and M 2 are the masses of the two binary components. The coefficients of the beaming effect α i are given by where F ν is the emission spectrum of the star in frequency units integrated over the Kepler bandpass. In our analysis, we used PHOENIX stellar models (Husser et al. 2013) to estimate the expected amplitude of the beaming effect. If both of the stars are on the main sequence, the luminosity depends strongly on mass (L ∝ M γ , where γ is close to 4). It follows that the theoretical amplitude in Equation (10) should always be positive in the case of a main-sequence binary. Likewise, in a star-planet system, f 2 is negligible, so A 1 is once again expected to be positive.

Mutual Illumination
While stars are poor reflectors overall, the temperature difference between the primary and secondary as well as their large surface areas can lead to nonnegligible mutual illumination effects. This effect is interpreted as heating of the stellar atmosphere, which in turn raises the local brightness and produces photometric modulations as the two components orbit each other. Kopal (1959) showed that the mutual illumination is manifested in the phase curve as a series of cosine terms, with the first harmonic (B 1 ) being the dominant term. Depending on the surface brightness ratio, the value of B 1 stemming from mutual illumination is expected to be negative for typical main-sequence binaries, but maybe positive for systems where the secondary has a significantly higher luminosity than the primary, such as in the case of binaries consisting of a hot white dwarf orbiting a main sequence star (e.g., Wong et al. 2020a).

KEPLER PLANETS RECLASSIFIED AS STELLAR COMPANIONS
Below, we present the results of our full phase-curve analysis for Kepler-854 b, primary transit and secondary eclipse fit for Kepler-840 b, as well as the transit lightcurve fit for Kepler-699 b and Kepler-747 b; the latter two systems have relatively long periods and hence do not show detectable phase-curve modulations. As for Kepler-840 b, we found that stellar activity strongly impacts the phase curve. The detailed data processing and fitting procedures used for the phase-curve and eclipse modeling analysis are described in Section 4.2. The phase-folded light curves for these four systems are shown in Figure 2; the best-fit parameter values are listed in Table 1.

Kepler-854 b
Kepler-854 b (KIC 7532973, TIC 275570329) was validated as a planetary candidate in Morton et al. (2016) (hereinafter M16). However, based on the range of estimated masses derived from the measured ellipsoidal effect in the Kepler photometry, we conclude that this companion is not a planet. In our analysis, we obtained statistically identical phase-curve amplitudes when splitting the Kepler light curve into two halves, indicating a consistent modulation across the time-series. We measured the amplitude of the leading term of the ellipsoidal variation (i.e., the second harmonic of the cosine) to be B 2 = −3519 ± 30 ppm. Given the mass of the host star listed in the TIC (0.84 ± 0.04 M ) we used Equation (4), a linear limb-darkening coefficient value of 0.581±0.040, and a gravity-darkening coefficient value of 0.328 ± 0.017 from Claret & Bloemen (2011), to derive a companion mass of 102 ± 12 M Jup , i.e., likely a late-type dwarf star. This simplified calculation, which ignored the contribution of the companion's tidal distortion to the overall measured amplitude at the second harmonic, was supported by the low fitted relative brightness of the companion (f 2 = 670 ± 120 ppm) and the small radius ratio (R 2 /R 1 = 0.1212 ± 0.0019). The radius ratio derived here is consistent with M16.
The revised stellar radius bumps the host star's radius from 1.25 +0.51 −0.17 R in M16 to 2.56 ± 0.05 R in our isochrone fitting, correspondingly raising the planet radius from 1.50 +0.60 −0.20 R Jup to 3.01 ± 0.07 R Jup (see Table 1). The bottom-right panel of Figure 1 shows the planetary radius values of confirmed planets as reported in NASA Exoplanet Archive against the corresponding insolation flux. There is a strong correlation between planetary radius and insolation flux beyond ∼2×10 5 W/m 2 (e.g., Demory & Seager 2011), but even among these inflated hot Jupiters, the physical radius of Kepler-854 b at 3.01 ± 0.07 R Jup stands out as being unusually large (Sestovic et al. 2018), and therefore it is unlikely that Kepler-854 b is a planet. We note that our derived companion mass is low for an object of this size, which is expected to have a mass of 289 ± 36 M Jup according to the relationships in Chen & Kipping (2017). It is thus possible that the mass obtained from our isochrones analysis could be underestimated. If we instead use the stellar radius from the SED fitting, which is more accurately constrained than the stellar mass, and combine that with the density of 165 ± 18 kg m −3 inferred from the transit fit, we obtain a host-star mass of 1.97 ± 0.25 M , which would yield a companion mass of 239 ± 54 M Jup -consistent with the expected mass of an object with the measured radius. In any case, we can rule out the planetary nature of Kepler-854 b at a high level of confidence.

Kepler-840 b
The strongest evidence against the planetary nature of Kepler-840 b (KIC 11517719, TIC 27847307) comes from the revised physical radius of the planet. We report an updated host-star radius of 1.73 ± 0.11 R from isochrone analysis (up from 1.06 +0.24 −0.12 R in M16), which correspondingly increases the measured planet radius from 1.52 +0.35 −0.17 R Jup to 2.50 ± 0.16 R Jup . The revised radius is too large for a substellar object (see bottom-right panel of Figure 1).
Additional evidence for the non-planetary nature comes from a large secondary eclipse depth of 1060±220 ppm. Assuming reflection as the primary source of the secondary eclipse, this entails a geometric albedo of 1.37 ± 0.31. This inferred geometric albedo is much higher than the values measured for hot Jupiters, which are close to ∼0.1 (e.g. Esteves et al. 2015;Angerhausen et al. 2015;Wong et al. 2020bWong et al. , 2021, or any solar system objects (Dyudina et al. 2016). If we take the secondary eclipse to stem entirely from thermal emission, then the temperature of the secondary is 3455 ± 150 K, which is inconsistent with the equilibrium temperature of 2465 ± 105 K, when assuming zero Bond albedo and a recirculation efficiency ( ) of 0 (Cowan & Agol 2011). Thus, this further supports the conclusion that the companion is a self-luminous object, i.e. a star, and not a planet. We attempted to measure the phase curve but found that it is strongly impacted by stellar activity. Specifically, the phase-curve fits from the two halves of the Kepler time series are not consistent with one another.

Kepler-699 b
Kepler-699 b (KIC 6061119, TIC 168811723) has an orbital period of 27.8 days, and at such a distance the The values presented in the bracket is the error in the last 1-3 digits of the reported values. The reported error may straddle across the decimal point.
phase-curve signal is not strong enough to be robustly detected in the Kepler light curve. We instead fit the transits only. With the revised parallax (Gaia EDR3), our isochrone analysis revealed the stellar radius to be 1.71 ± 0.10 R , up from 0.780 +0.077 −0.067 R reported in M16. This increase in the stellar radius yields a planet radius of 2.76 ± 0.16 R Jup , unphysical for a planetary body, especially given that the companion is not highly irradiated and is therefore unlikely to be inflated to such an extent (see bottom-right panel of Figure 1).

Kepler-747 b
Kepler-747 b (KIC 5780460, TIC 121602473) is a longperiod planet (35.6 days) that requires reexamination in light of the newly available parallax measurement of the host star. The updated stellar radius is 3.20 ± 0.29 R , revised from 0.780 ± 0.042 R in M16, a factor of four increase. This leads to a companion radius of 1.84 ± 0.17 R Jup , an unlikely value given the irradiation level of the companion, as shown in the bottom-right panel of Figure 1. We note that the three planets in the vicinity of Kepler-747 b in the plot are Kepler-447 b (Lillo-Box In the absence of a clear secondary eclipse, as in the case of Kepler-840 b, we cannot completely rule out Kepler-747 b as a planet based on its radius alone.

PHASE-CURVE ANALYSIS OF KEPLER ECLIPSING BINARIES
In the previous section, we found that for Kepler-854 b, the complementary insights provided by the phase curve offer a powerful new perspective into the companions' properties. The phase curve of Kepler-854 shows a significant phase-curve modulation at the third harmonic, with B 3 = −215 ± 32 ppm (Table 1). Such highorder cosine harmonics have been detected in the phase curve of two planetary systems -Kepler-13A and HAT-P-7 (Shporer et al. 2014;Esteves et al. 2015) -though at amplitudes that are at least an order of magnitude smaller. While it has generally been assumed that such higher harmonics are manifestations of the higher-order terms in the host star's ellipsoidal variation (see Section 2.2.1), the astrophysical nature of these signals has not been investigated empirically for a larger sample (for a numerical study, see Gomel et al. 2021).
Due to the lack of such systematic studies and the fact that only a few planetary systems have shown higher than expected amplitudes for the third cosine harmonics, the origins of such harmonics have come under debate. An earlier attempt to look at these higher harmonics among Kepler Objects of Interests by Armstrong & Rein (2015) suffered from both contamination from stellar activities and artifacts from the data analysis (Cowan et al. 2017). This in turn inspired alternative explanations, such as consequences of eccentricity/inclination, pipeline artifacts, and effects of obliquity or weather, although tidal forces remain the most plausible explanation (Cowan et al. 2017;Bielecki & Cowan 2018). The planetary phase curve picture is also further complicated by possible temporal weather patterns (Menou & Rauscher 2009;Armstrong et al. 2016;Cho et al. 2021) or asymmetric cloud coverage (Demory et al. 2013;Shporer & Hu 2015), which can lead to power leakage from lower stronger harmonics (often B 1 for planets) to higher order terms. Given the long observational baseline for Kepler, these variability should however average out. In this section, we supplement our previous analysis of misclassified planetary systems with an investigation of a sample of short-period Kepler eclipsing binaries, particularly those with high signal-to-noise phase-curve modulations. Analyzing these signals may shed new light on the origin of these harmonics and help contextualize similar high-harmonic phase-curve signals in other photometric observations.

Target Selection
We started with the Kepler eclipsing binary (EB) catalog, which contains 2,920 EBs 1 (Prša et al. 2011;Conroy et al. 2014;Kirk et al. 2016). From that list, we selected EBs with orbital periods between 1.5 and 10 days, which are expected to be detached while expecting strong phase-curve signals. We also established a maximum Kepler magnitude of 13.0 to ensure adequate photometric precision. Next, we removed systems that do not show clear eclipses (i.e., contact binaries) and systems with eccentric orbits (i.e., where the secondary eclipse is not half an orbital period away from the primary transit).
To ensure the presence of a coherent phase curve signal in the observed light curve, we divided the light curve into two halves (just as we did for the Kepler-854 and Kepler-840 light curves) and compared the binned phase curve from the two halves, removing any systems which showed observable discrepancies. This test helps us to assess the impact of stellar activity on the phase curve. The Kepler EB catalog was supplemented by other published studies of Kepler eclipsing binaries (Rappaport et al. 2015;Faigler et al. 2015). We note that the target selection criteria were relaxed somewhat for a few targets, such as KIC 4169521, a white dwarf system with a previously studied phase curve (Faigler et al. 2015). Some of the targets show increased scatter in the residuals around the eclipses, which are potentially due to eclipse-timing variations, spot-crossing events, or imperfectly corrected dilution in some Kepler quarters (Stumpe et al. 2012). These targets were also removed.
We performed a preliminary fitting analysis on around 60 EB systems, out of which 20 had both a radius ratio less than 50% and a secondary flux contribution less than 10%. As discussed in Section 2.2, when R 2 /R 1 and f 2 are small, the contribution of the secondary companion to the overall phase-curve signal is negligible, allowing us to focus on the physical properties of the primary. These 20 systems we selected for our analysis are shown in Figure 3 in the context of the ∼200,000 Kepler targets. It is evident that our targets are biased toward hotter host stars. This bias could be introduced due to the likelihood of the cooler binary pairs developing co-rotating spots, which can destroy the long-term coherency of the phase-curve signal (Giles et al. 2017;Rackham et al. 2019). The distribution of our sample of 20 Kepler eclipsing binary systems (red circles) plotted on top of a color map showing the distribution of ∼200,000 targets from the Kepler mission, which primarily looked for planets around Sun-like stars. For the Kepler targets, the plotted parameters were taken from the KIC, whereas for our selected 20 EB systems, we used the more recent values from the TIC.

Model Fitting
We downloaded all of the long cadence photometry available for each target from the Mikulski Archive for Space Telescopes (MAST) and used the latest processed light curves (DR 25; Jenkins et al. 2010Jenkins et al. , 2017. For our fitting analysis, we took the simple aperture photometry (SAP) light curve provided in the data files and discarded all data points with a nonzero quality flag. We corrected the light curves for long-term systematic trends by applying a cubic spline with knots tied every orbital period to the transit-masked photometry (Niraula et al. 2018;Daylan et al. 2021). To clean the photometry of flux ramps and any edge effects incurred by the spline detrending process, we removed the data affected by these spurious features that typically occurred near gaps longer than 6 hours. Finally, to further identify any remaining outliers or significant noise features, we phase-folded the light curve on the orbital period to build a template of a coherent periodic signal. We then used this template to construct the common-mode time series signal throughout the entire light curve and removed all data points that deviated by more than 3.5σ. This step removed around 3% of the time series, on average.
When examining the Lomb-Scargle periodograms of the cubic spline models, we found that they include signals at the first harmonic of the orbital frequency at the level of ∼1% of the total ellipsoidal amplitudes, on average. For targets with large photometric variations, such as KIC 8906039, this power leakage is comparable to or even exceeds the uncertainties on the measured amplitudes. A workaround was therefore developed to address this issue, which we describe later in this section.
To characterize the system, we performed a simultaneous fit of the phase curve, primary transits, and secondary eclipses. The primary transits and secondary eclipses were fit using ellc (Maxted 2016), with a supersampling factor of 10 to account for Kepler's longcadence integration time of 29.4 minutes. Our model accounts for the finite light-travel time across the orbit, which affects the measured secondary eclipse phase relative to the transit phase. We fit for the first five harmonics of both cosine and sine modulations, as described in Equation (2). As the input to our fitting pipeline, we used the median-binned, phase-folded time series, a process that expedites the fitting process relative to fitting the full unbinned, non-phase-folded photometry and naturally removes any significant non-periodic signals.
We used the Markov chain Monte Carlo (MCMC) ensemble sampler emcee (Foreman-Mackey et al. 2013) with 22 parameters, including the radius ratio R 2 /R 1 , mid-transit time T 0 , impact parameter b, scaled semimajor axis a/R 1 , the luminosity ratio S 2 /S 1 (which is directly related to the brightness ratio f 2 ), orbital period P , eccentricity parameters e sin ω and e cos ω (where e is the orbital eccentricity and ω the argument of periastron), quadratic limb-darkening coefficients for both binary components (4 parameters total), and the phasecurve amplitudes (10 parameters). The host stars were assumed to be spherical; any deviation from sphericity due to ellipsoidal tidal distortion is expected to be at the level of a few percent at most, introducing negligible biases to our retrieved parameters. As for the limb darkening, we used the uninformative sampling technique and reparametrization developed by Kipping (2010). We initialized the period and transit time at the values provided by the NASA Exoplanet Archive.
To explore the parameter space, we ran the MCMC with 44 walkers (twice the number of parameters) for about 30,000 steps each. We ensured that the walkers were well-mixed. Once we confirmed that the fits were satisfactorily converged, we removed the burn-in phase. The burn-in phase was identified by visually examining the evolution of the log probability with the number of steps and determining the point when the walkers converged at the global log probability maximum. This typically occurred after about 15,000 steps, and we discarded the chain segments prior to this point before constructing the posterior distributions.
In order to address the aforementioned contamination of the measured phase-curve signal from the spline de-trending process, we tried fitting the uncorrected, nonphase-folded light curves, after masking the transits and eclipses. We used the priors on the transit ephemeris from the previous full phase-curve fits and then ran an additional set of MCMC fits for all 20 systems. The typical uncertainties on the phase-curve amplitudes from fitting the undetrended light curves were larger by a factor of two when compared to the detrended, phasefolded light-curve fits. However, the methods showed broad mutual consistency, except for a few cases such as KIC 8906039. In the Kepler EB phase-curve results presented in this paper, we choose the more conservative approach and report the best-fit amplitudes from the undetrended light-curve fits.

Results
For our 20 selected systems, the phase-folded light curves with best-fit models overplotted are shown in Figures 7 and 8, and the fitted parameters are listed in Tables 2 and 3. As reported in these tables, the first cosine harmonic amplitudes (B 1 ), which is dominated by the mutual illumination, are typically negative. However, some systems show positive B 1 values. These systems all have brightness ratios (S 2 /S 1 > 1) and correspond to the EBs that have hot white dwarf secondaries. This behavior is expected, as the sign of the mutual illumination amplitude is driven primarily by the temperature ratio of two binary components. As for the higher-order cosine amplitudes, we find clear empirical power-law relations between the second (B 2 ) and third (B 3 ) cosine harmonics (left panel of Figure 4) as well as between the second and fourth (B 4 ) cosine harmonics (right panel of Figure 4). Such a positive correlation is expected from Equation (4), as all the cosine harmonics scale proportionately with the mass ratio, and inversely with the scaled distance. Our data set encompasses measured phase-curve amplitudes that span two orders of magnitude, often recommended for claiming robust power law. The scatter can be attributed to the range of inclinations, gravity-and limb-darkening coefficients, as well as the differing contributions from the secondary companions. The observed power law suggests a common provenance for the signals in these systems and strongly implies that the higher cosine harmonics (both B 3 and B 4 ), just as the second cosine harmonic (B 2 ), originate from the mutual tidal interaction between the binary components.
In order to find the best-fit power law, we used the standard distance regression routine in scipy. To estimate the uncertainties in the fit parameters, we performed a fit to half of the data points at a time, carried out the regression on all the possible subsets, and . Left: a comparison of the measured second-(B2) and third-harmonic (B3) cosine amplitudes among our sample of Kepler eclipsing binaries and planet-hosting systems, demonstrating a strong power-law relation between the two. The symbols indicate the method by which the system mass ratio q was determined -from the radial velocity signal or the simultaneously measured beaming amplitude A1. The color scheme corresponds to the effective temperature of the primary star. The error bars are typically smaller than the markers. Right: same, but for the second-and fourth-harmonic cosine amplitudes. The two targets for which the signals are not statistically significant are not included. calculated the median and standard deviation from the resultant distribution of fitted values. In the left panel of Figure 4, we include Kepler-854, which is positioned close to the fitted relation, further highlighting the nonplanetary nature of the phase curve, as planets primarily occur at the bottom-left corner of the plot. Notably, the planetary systems located at the bottom-left corner in Figure 4 deviate from the expected empirical law, which could be due to additional unconsidered phenomena, in-cluding those mentioned previously at the beginning of this section. Figure 5 presents a comparison between the measured and the corresponding expected values given the system parameters for the second (left panel) and third harmonic (right panel), assuming they originate from the ellipsoidal variation as described in Section 2.2.1. For each system, the mass ratio q was estimated based on the published radial velocity orbits (if available) or by Effective Temperature (T 1 ) [K] Figure 6. Observation of the expected linear relation between the ratio of second-and third-harmonic cosine amplitudes |B2/B3| against a term containing the scaled semimajor axis, inclination, limb-darkening coefficients, and gravitydarkening coefficients (see Equation (9)). The best-fit line has a slope 2.24 ± 0.48, shown in blue, which is consistent with the theoretically expected value of 2.4. The shaded region shows the 1σ range of the fitted model. Three targets with statistically insignificant third harmonics are excluded from the fit.
deriving it from the measured beaming amplitude (see Section 2.2.2 and also Equation 1 in Shporer et al. 2010). The ellipsoidal and beaming coefficients were derived using the stellar parameters from the TIC (Stassun et al. 2019). Both panels in Figure 5 show good agreement between the measured and expected values, with consistency at a level close to or better than 1σ for most systems. A few systems do show significant discrepancies (≥ 3σ) between the measured and expected values of B 2 in Figure 5. These include KIC 3735597 and 6449358, where the observed amplitude is larger than the expected value. Within the precision of our data, we do not see more significant discrepancies for the hotter stars, as was predicted by Pfahl et al. (2008). For all of these discrepant systems, the mass ratio was derived from the measured beaming amplitude (A 1 ). Such deviations have been reported for systems such as KIC 4169521 (Faigler et al. 2015) and KIC 9164561 (Rappaport et al. 2015), which exhibit statistically significant differences between the mass based on the beaming amplitude and the mass derived from the radial velocity signal. This is perhaps due to the presence of a phase shift in the mutual illumination signal (Rappaport et al. 2015), which would manifest as a systematic bias in the measured A 1 amplitude and translate into an incorrect beaming-derived mass estimate. For KIC 3735597 and 6449358, a systematic underestimation of the mass from beaming is a plausible explanation, as they simultaneously have larger than expected B 3 measurements.
Additional targets such as KIC 4169521, 6606653, and 9786017 also show disagreement between the observed and measured B 3 amplitudes. It appears that while the classical estimation provides good order-of-magnitude predictions across the board, additional astrophysical phenomena may be contributing nonnegligibly towards the second-and third-harmonic amplitudes in some individual systems. Given that obtaining a precise and accurate mass ratio estimate is difficult, we explored alternative pathways to further test the veracity of the classical formulation of ellipsoidal variation. As outlined in Section 2.2.1, we compared the expected ratio between B 2 and B 3 , given in Equation (9), to our measured value. The advantage of using |B 2 /B 3 | is that it is approximately independent of the mass ratio for systems with relatively small and faint secondaries. Meanwhile, a/R 1 and sin i are fitted parameters from our light-curve analysis. In Figure 6, we plot |B 2 /B 3 | vs. a/R 1 / sin i · (X 2 G 2 )/(X 3 G 3 ), where the expected relation should be linear with a slope of 2.4, following Equation (9). We fitted a linear relation that yielded a slope of 2.24±0.48 and an intercept that is consistent with zero. Despite a large scatter, this value is statistically consistent with 2.4, the value expected from the classical formalism. We did not fix the y-intercept to zero, as it leads to an underestimation of the slope uncertainty.

SUMMARY AND CONCLUSIONS
In this paper, we showed that three statistically validated transiting planets, Kepler-854 b, Kepler-840 b, and Kepler-699 b cannot be planets but are instead small stars. We used updated host-star parameters from the TIC and stellar SED fitting to show that the measured transit depth indicates a companion radius that is too large for a planet. For Kepler-854 b, we also used the Kepler phase curve to show that the transiting object's mass is inconsistent with its being a planet. In addition, we demonstrated that a fourth statistically validated planet, Kepler-747 b, has a somewhat larger radius than other confirmed gas-giant planets with the same stellar irradiation level, suggesting that it too is not a planet, but a small star.
This work demonstrated the use of phase curves as a validation tool for Kepler-854 b, as the phase curve amplitude is sensitive to the transiting object's mass (more directly, the system mass ratio). Meanwhile, the radius is broadly unchanged for objects from about 1 M Jup (∼0.001 M ) to about 100 M Jup (∼0.1 M ) (e.g., Chen & Kipping 2017), making statistical validation methods that primarily rely on the companion radius susceptible to false positive scenarios. Fortunately, orbital phase-  (14) References :  curve modulations can be detected automatically for a large number of objects using the techniques described here and elsewhere (e.g., Shporer et al. 2011;Faigler & Mazeh 2011;Niraula et al. 2018;Wong et al. 2020b). Our work also shows the critical role that accurate host-star parameters play in validating transiting planet candidates (see also Shporer et al. 2017). Here, we have benefited from measurements carried out by the ESA Gaia Mission (Gaia Collaboration et al. 2018), which were not available during the Kepler mission and have since dramatically improved the quality of stellar parameters for planet-hosting stars and, in turn, decreased the number of false positive candidates.
Following the detection of a third-harmonic signal in the phase curve of Kepler-854, we used a sample of Kepler EBs to empirically investigate the nature of such higher-order signals. We found that the amplitude of the measured signal from Kepler-854 is consistent with the behavior of the Kepler EB sample. Moreover, the higher-order harmonic amplitudes measured from the EB sample show excellent agreement with the expected values from theoretical modeling of the tidal distortion of the primary stars. While perhaps somewhat expected, to the best of our knowledge, such an empirical study has not been done before, and our results provide us with a better understanding of not only the Kepler-854 phase curve, but also the contribution of ellipsoidal variation to visible-light phase curves in general.
With the increasing number of transiting planet candidates from the NASA TESS Mission (Guerrero et al. 2021) and the already limited amount of observational resources, statistical validation of candidates is becoming more and more essential. Some of the lessons learned in this work will help improve statistical validation techniques, which will, in turn, support more accurate studies of planet demographics.