Supervised Learning Detection of Sixty Non-Transiting Hot Jupiter Candidates

TL;DR

This paper introduces a supervised learning algorithm to detect non-transiting hot Jupiters in Kepler data by analyzing optical phase curves, identifying 60 high-probability candidates for further confirmation and study.

Contribution

The study develops a novel supervised learning method to find non-transiting hot Jupiters using phase curve data, expanding detection capabilities beyond transiting planets.

Findings

Identified 60 high-probability non-transiting hot Jupiter candidates.

Classified 142,630 Kepler stars for potential phase curve signals.

Derived constraints on albedos and phase curve offsets for candidates.

Abstract

The optical, full-phase photometric variations of a short-period planet provide a unique view of the planet's atmospheric composition and dynamics. The number of planets with optical phase curve detections, however, is currently too small to study them as an aggregate population, motivating an extension of the search to non-transiting planets. Here we present an algorithm for the detection of non-transiting, short-period giant planets in the Kepler field. The procedure uses the phase curves themselves as evidence for the planets' existence. We employ a supervised learning algorithm to recognize the salient time-dependent properties of synthetic phase curves; we then search for detections of signals that match these properties. After demonstrating the algorithm's capabilities, we classify 142,630 FGK Kepler stars without confirmed planets or KOIs and, for each one, assign a probability of a phase curve of a non-transiting planet being present. We identify 60 high-probability non-transiting hot Jupiter candidates. We also derive constraints on the candidates' albedos and offsets of the phase curve maxima. These targets are strong candidates for follow-up radial velocity confirmation and characterization. Once confirmed, the atmospheric information content in the phase curves may be studied in yet greater detail.

Tables

Table 1. Table 1: Parameter Distributions of Synthetic Phase Curve Injections

Parameter	distribution
period, $P$ [days]	$P\sim 𝒰[1,5]$
planet mass, $M_P$ [$M_{\mathrm{Jup}}$]	$M_P\sim 𝒰[0.25,2.5]$
planet radius, $R_P$ [$R_{\mathrm{Jup}}$]	$R_P\sim 𝒰[0.7,1.8]$
inclination, $i$ [∘]	$i\sim 𝒰[45,90]$
mean geometric albedo,
albedo covariance amplitude, $h_A^{}{}_{}{}^2$	$h_A^{}{}_{}{}^2\sim 𝒰[0,0.015]$
mean peak offset,
$\theta$ covariance amplitude, $h_{\theta }^{}{}_{}{}^2$
coherence timescale, $\tau$ [days]	$\tau \sim 𝒰[10,50]$

Table 2. Table 2: A catalog of 60 high probability non-transiting hot Jupiter candidates

KIC	$K_p$	$T_{\mathrm{eff}}\left[K\right]$	$R_{\star }\left[R_{\odot }\right]$	$M_{\star }\left[M_{\odot }\right]$	Fe/H	$P$ [days]	$A_{\mathrm{refl}}$	$A_{\mathrm{ellip}}$	$K\sin i\left[m/s\right]$	Prob
\csvreader[late after line=
\probability

Table 3. Table 4: A catalog of 60 high probability non-transiting hot Jupiter candidates

KIC	$K_p$	$T_{\mathrm{eff}}\left[K\right]$	$R_{\star }\left[R_{\odot }\right]$	$M_{\star }\left[M_{\odot }\right]$	Fe/H	$P$ [days]	$A_{\mathrm{refl}}$	$A_{\mathrm{ellip}}$	$K\sin i\left[m/s\right]$	Prob
\csvreader[late after line=
\probability

Equations

$$d=\frac{\sum _{i=2}^N(r_i-r_{i-1}{)}^2}{\sum _{i=1}^N{r_i}^2}.$$

$$Q_h=N(N+2)\sum _{j=1}^h\frac{{\widehat{{\rho }_j}}^2}{N-j},$$

$$\underset{P\in \{P_1,...,P_4\}}{\min }\left\{\frac{{{\chi }^2}_2}{{{\chi }^2}_1}-\frac{d_2}{d_1}+\frac{m_2}{m_1}\right\}.$$

$$F=F_p(\varphi +\theta )+F_d(\varphi )+F_e(\varphi ).$$

$$\varphi =\frac{t-T_{\mathrm{mid}}}{P},$$

$$F_p=A_p\frac{\sin z+(\pi -z)\cos z}{\pi },$$

$$A_p\approx F_{\mathrm{ecl}}=A_g{\left(\frac{R_p}{a}\right)}^2,$$

$$F_d={\left(\frac{2\pi G}{P}\right)}^{1/3}\frac{M_p{\alpha }_d\sin i}{c{M_{\star }}^{2/3}}\left(\frac{1+e\cos \omega }{\sqrt{1-e^2}}\right)\sin (2\pi \varphi ).$$

$$F_{e}=-\frac{M_{p}\alpha_{2}\sin^{2}i}{M_{\star}}\biggl{[}\left(\frac{a}{R_{\star}}\right)^{-3}\cos(4\pi\phi)\ +\\ \left(\frac{a}{R_{\star}}\right)^{-4}f_{1}\cos(2\pi\phi)+\left(\frac{a}{R_{\star}}\right)^{-4}f_{2}\cos(6\pi\phi)\biggr{]},$$

$$\begin{array}{rl}{f}_1& =3{\alpha }_1\frac{5{\sin }^{2}i-4}{\sin i}\\ f_2& =5{\alpha }_1\sin i\end{array}$$

$$k(t_i,t_j)=h^2\exp \left(-\frac{(t_i-t_j{)}^2}{2{\tau }^2}\right).$$

$$\mathbf{K}=\left[\begin{array}{cccc}k(t_1,t_1)& k(t_1,t_2)& \dots & k(t_1,t_n)\\ k(t_2,t_1)& k(t_2,t_2)& \dots & k(t_2,t_n)\\ ⋮& ⋮& \ddots & ⋮\\ k(t_n,t_1)& k(t_n,t_2)& \dots & k(t_n,t_n)\end{array}\right].$$

$$x_{Ai}=x_A(t_i)=\frac{A_g(t_i)-⟨A_g⟩}{⟨A_g⟩}.$$

$$E[Y_i|x_{1i},x_{2i},...,x_{mi}]=p_i.$$

$$Y_i|x_{1i},x_{2i},...,x_{mi}\sim \mathrm{Bernoulli}(p_i),$$

$$P(Y_i=y_i|x_{1i},x_{2i},...,x_{mi})={p_i}^{y_i}(1-p_i{)}^{1-y_i}.$$

$$\mathrm{logit}(p_i)=\ln \left(\frac{p_i}{1-p_i}\right)=\mathbf{\beta }\cdot {\mathbf{x}}_{\mathbf{i}}.$$

$$\mathbf{\beta }=({\beta }_0,{\beta }_1,{\beta }_2,...,{\beta }_m),$$

$${\mathbf{x}}_{\mathbf{i}}=(1,x_{1i},x_{2i},...,x_{mi}).$$

$$p_i={\mathrm{logit}}^{-1}(\mathbf{\beta }\cdot {\mathbf{x}}_{\mathbf{i}})=\frac{1}{1+e^{-\mathbf{\beta }\cdot {\mathbf{x}}_{\mathbf{i}}}}.$$

$$\begin{array}{rl}l(\mathbf{\beta })& =\ln \left[\prod _{i=1}^{n}p(y_i|{\mathbf{x}}_{\mathbf{i}})\right]\\ & =\sum _{i=1}^n\left[(1-y_i)\ln (1+e^{\mathbf{\beta }\cdot {\mathbf{x}}_{\mathbf{i}}})-y_i\ln (1+e^{-\mathbf{\beta }\cdot {\mathbf{x}}_{\mathbf{i}}})\right].\end{array}$$

$$\begin{array}{r}f(\varphi )=1-A_{\mathrm{refl}}\cos (2\pi (\varphi +{\theta }_{\mathrm{refl}}))\\ -A_{\mathrm{ellip}}\cos (4\pi (\varphi +{\theta }_{\mathrm{ellip}}))\\ +A_{\mathrm{beam}}\sin (2\pi (\varphi +{\theta }_{\mathrm{ellip}})).\end{array}$$

$$R_p={\left(\frac{2A_{\mathrm{refl}}}{A_g}\right)}^{1/2}{\left(\frac{G{M}_{\star }}{4{\pi }^2}\right)}^{1/3}{P}^{2/3}(\sin i{)}^{-1/2}.$$

$${\rho }_p=\frac{{\rho }_{\star }}{2\sqrt{2}}\left(\frac{A_{\mathrm{ellip}}}{{A_{\mathrm{refl}}}^{3/2}}\right)\left(\frac{{A_g}^{3/2}}{{\alpha }_2}\right){\left(\sin i\right)}^{-1/2}.$$

$$S={\iint }_{D}f(R_p,\ln {\rho }_p)d{R}_{p}d\ln {\rho }_p,$$

$$D=\{(R_p,\ln {\rho }_p):A_g\in [0.1,0.3],i\in [15^{\circ },85^{\circ }]\}.$$

$$K\sin i=\frac{A_{\mathrm{ellip}}}{{\alpha }_2}{\left(\frac{G{M}_{\star }}{a}\right)}^{1/2}{\left(\frac{a}{R_{\star }}\right)}^3$$

$$A_g{R_p}^2\sin i=2A_{\mathrm{refl}}{\left(\frac{G{M}_{\star }}{4{\pi }^2}\right)}^{2/3}{P}^{4/3}.$$

$$T_{\mathrm{irr}}=T_{\mathrm{eff}}{\left(\frac{R_{\star }}{a}\right)}^{1/2},$$

$$T_{\mathrm{eq}}=T_{\mathrm{eff}}{\left(\frac{R_{\star }}{a}\right)}^{1/2}{\left[f(1-A_B)\right]}^{1/4},$$

Full text

1 Introduction

While the Kepler prime mission is best known for the detection of thousands of transiting extrasolar planets, notably including some that are Earth-sized and smaller (e.g. Batalha et al., 2013), its long time baseline and high photometric precision allowed for the detection of other small-amplitude signals, such as out-of-transit photometric variations of short-period planets (Borucki et al., 2009). Long before Kepler was launched, the existence and detectability of these optical reflected light phase curves were predicted by Seager et al. (2000) and Jenkins & Doyle (2003).

Optical phase curves, along with their infrared counterparts, provide a spatially-integrated, time-dependent view of the planet’s atmosphere. They yield information about the planet’s atmospheric composition (Rowe et al., 2006; Demory et al., 2013), day-night temperature contrasts and heat redistribution (Knutson et al., 2007; Showman et al., 2009; Cowan & Agol, 2011), cloud existence and reflectivity (Demory et al., 2013; Shporer & Hu, 2015), atmospheric weather variability (Armstrong et al., 2016), and magnetic field strength (Rogers, 2017). A recent review of phase curve methodology and scientific findings is provided by Shporer (2017).

A planet’s optical phase curve is the superposition of four independent components: reflected star light, thermal emission, ellipsoidal variations, and Doppler beaming/boosting (Faigler & Mazeh, 2011). Ellipsoidal variations are brightness changes resulting from tidal deformation to the star (Morris, 1985), while Doppler beaming is a relativistic effect that causes an increase or decrease in light as the star approaches or recedes from the observer (Hills & Dale, 1974; Rybicki & Lightman, 1979). Ellipsoidal variations and Doppler beaming are both sensitive to the planet’s mass, which thereby permits photometric mass constraints (e.g. Lillo-Box et al., 2016). For an average hot Jupiter ($P\sim 3$ days, $R_{p}\sim 1.3\ R_{\mathrm{Jup}}$, $M_{p}\sim M_{\mathrm{Jup}}$) orbiting a Sun-like star, the reflected light component is typically larger than the other components by about an order of magnitude.

There have been numerous detections of optical phase curves of transiting Kepler planets (e.g. Welsh et al., 2010; Shporer et al., 2011; Barclay et al., 2012; Mazeh et al., 2012; Lillo-Box et al., 2014). In addition, comprehensive searches for Kepler planetary phase curves and secondary eclipses have also been published (Coughlin & López-Morales, 2012; Esteves et al., 2013; Angerhausen et al., 2015; Esteves et al., 2015; Lillo-Box et al., 2016). These systematic searches have established several interesting findings, including generally low (but occasionally quite large) hot Jupiter albedos (Demory et al., 2011) and repeated detections of shifts in the maximum of the phase curve away from the sub-stellar point (Faigler et al., 2013; Faigler & Mazeh, 2015; Shporer & Hu, 2015; Hu et al., 2015).

Though it is substantial, the sample size ($\sim$ 15) of Kepler transiting planets with detectable photometric variations is just short of the count necessary to adequately understand these atmospheric properties in a statistical or populational sense. It is therefore worthwhile to increase the number of planets with well-characterized optical phase curves by expanding the search to non-transiting planets and using the phase curve signal itself as a detection mechanism. Given their large expected radial velocity half-amplitudes, candidates identified in this manner can subsequently be confirmed with Doppler velocity measurements.

The prospects of detecting non-transiting planets from their photometric variations alone has been frequently discussed in the literature. The BEaming, Ellipsoidal, and Reflection/heating (BEER) model (Faigler & Mazeh, 2011) was developed in part as a technique to discover non-transiting planets. To date, the BEER model has been very successful in the detection of non-transiting stellar binary companions (Faigler et al., 2012) and the characterization of transiting planets with detectable phase curve variations (Faigler et al., 2013; Faigler & Mazeh, 2015).

Placek et al. (2014) proposed using Bayesian model selection to test a light curve for the presence of significant reflected light, thermal emission, ellipsoidal variation, and/or Doppler beaming phase curve components. Millholland et al. (2016) considered the prospects of discovering non-transiting hot Jupiters in Kepler systems containing known transiting planets with the coupled detection of an optical phase curve and astrometric (stellar wobble-induced) transit timing variations.

In all of these techniques, the primary challenge is that the ellipsoidal and Doppler beaming components are typically significantly smaller than the reflected light component for Jupiter-sized planets. The photometric signature is therefore mostly sinusoidal and difficult to distinguish from stellar or instrumental variability. These prior studies, however, have uniformly considered the phase-folded and time-averaged phase curve signal. When the signal is folded and binned, valuable time-dependent information is lost. This normally discarded information can be helpful in distinguishing planetary signals from stellar ones.

In this paper, we present a systematic procedure to detect non-transiting, short-period giant planets from their phase curves by exploiting the phase curves’ time-dependent properties. The algorithm is rooted in two key ideas:

A planetary phase curve is more temporally consistent in its amplitude and phase than most other stellar or instrumental light curve variability. 2. 2.

The properties of phase curves of known transiting Kepler planets and synthetically generated datasets can be used as a training set on which novel light curves can be compared and classified using supervised learning.

This paper is organized as follows. In §2, we outline the pipeline we developed for detecting candidate phase curves in Kepler light curves. In §3, we inject synthetic phase curves into Kepler light curves and demonstrate the pipeline’s efficiency in recovering them. In §4, we train a logistic regression algorithm to distinguish the properties of planetary phase curves from other signals. We test its classification capabilities on synthetic datasets and the set of Kepler transiting hot Jupiters. In §5, we employ the supervised learning algorithm to search for phase curves around Kepler FGK stars without known planets or planet candidates. In §6, we establish a catalog of candidate planets with phase curves and examine trends in the planets’ albedos and phase curve maximum offsets. We discussion additional properties of the candidates in §7 and conclude in §8.

2 An Automated Detection Pipeline

To detect non-transiting hot Jupiters in the Kepler field, we first construct an algorithm for identifying the single best candidate phase curve within a given light curve. The process consists of three distinct steps:

Search for the periods of potentially significant, approximately sinusoidal signals in the data by calculating Lomb-Scargle periodograms of initially detrended light curves. 2. 2.

At each candidate period, find the best-fit parameters of a phase curve by simultaneously fitting for the phase curve and stellar components. 3. 3.

Quantify the autocorrelations in the residuals of the fits.

The light curves we consider here are the PDC (pre-search data conditioning) photometry (Smith et al., 2012; Stumpe et al., 2012, 2014) provided by the Kepler Science Center and publicly available at Mikulski Archive for Space Telescopes (MAST)111https://archive.stsci.edu/kepler/. The PDC light curves have undergone a removal of systematic errors common to all light curves.

2.1 Step 1- Identifying candidate periods

The first goal is to identify the periods of significant sinusoidal signals in the light curve. To do this, we must first damp the effects of stellar variability or residual instrumental systematics that could contaminate a periodogram. We model our detrending algorithm on the kepflatten routine in the PyKE software package222https://keplerscience.arc.nasa.gov/PyKE.shtml (Still & Barclay, 2012).

We detrend the light curve by fitting it with the mean of polynomials that span a sliding window. In detail, we first split the light curve into sections of length $s$. Then, over a window of size $w=Ns$ (with $N$ an integer), we fit a third degree polynomial. We slide the window through all time sections, stepping in units of $s$. Within each section, we then calculate the mean of the $N$ polynomials that were fit in that section, resulting in a smooth fitted curve. We then detrend the light curve data by dividing by the fit. At this stage, we calculate outliers that are deviant by more than $7\sigma$, remove them, and iterate the process.

We perform two detrending calculations with different window sizes in an attempt to account for stellar variability at a variety of timescales. A step size of $s=2$ days is used in each calculation, but the short-timescale detrending uses $w=6$ days, while the long-timescale detrending uses $w=12$ days. We aim to detect phase curves with orbital periods up to 5 days (though in reality, most detections will be less than 3 days). It is therefore critical that $w$ be large enough to avoid fitting out the phase curve signals we wish to detect. This explains our choice of $w=6$ for the short-timescale detrending.

On each of the two detrended curves associated with a given light curve, we generate a Lomb-Scargle periodogram (Lomb, 1976; Scargle, 1982) within a 1-5 day range. Hence, we assume that we will only be capable of detecting phase curves within this period range. Although there are several known hot Jupiters with sub-day periods (e.g. WASP-18b, Hellier et al., 2009), we found that extending the periodogram much below 1 day leads to significantly more noise and poorer overall results.

We also generate a “local significance” periodogram. At each period in the spectrum, we calculate the number of standard deviations of the periodogram power over the median power within a centered 0.5 day window. The purpose of this additional step is that the peak in a periodogram associated with a true phase curve is sometimes highly significant with respect to nearby powers but not with respect to the entire 1-5 day range. As an illustrative example, Figure 1 displays the Lomb-Scargle periodogram and local significance periodogram for the $R\sim 1.08\ R_{\mathrm{Jup}}$, $P\sim 1.595$ day, transiting hot Jupiter, Kepler-686b. We removed the transits and secondary eclipses in the light curve before the computation.

We extract the highest peak in each of the four periodograms – one Lomb-Scargle and one local significance periodogram for each of the two detrending calculations. We find the unique periods (defined as separate by more than 0.01 days) among these, yielding a collection of 1-4 peak periods. These period(s) constitute candidate periods for a planetary phase curve. The candidate signals are analyzed in more detail in the next section.

2.2 Step 2- Phase curve fitting

The initial detrending calculations described above are sufficient for identifying the periods of significant sinusoidal signals in the data. To robustly characterize the candidate phase curves, however, it is necessary to model the planetary signal and stellar variability simultaneously, rather than attempting to remove the stellar signal a priori. To this end, we use the following steps to find the best-fit sinusoidal signal for each of the candidate period(s) identified.

The best fit is found by chi-square minimization of a multiplicative, two-component (phase curve & stellar variability) light curve model. The phase curve component is modeled by a fixed-period sinusoid with variable phase and amplitude. The stellar variability component, that is, the component remaining after the sinusoid is divided out, is fit using the sliding polynomial technique described in Section 2.1. Here we use three degree polynomials, a step size $s=2$ days, and window size $w=6$ days. The best-fit is characterized by the amplitude and phase of the sinusoid that minimizes the residuals of the joint fit. For the purposes of later comparison, we also perform a one-component fit without the sinusoid. In other words, this just contains the sliding polynomial fit. An example of the one-component (star only) and two-component (star & phase curve) fits are shown for Kepler-13b in Figure 2. Kepler-13b is a $R\sim 1.51\ R_{\mathrm{Jup}}$, $P\sim 1.764\ \mathrm{day}$, transiting hot Jupiter with a large amplitude phase curve.

We note that techniques such as Gaussian Process (GP) regression or autoregressive integrated moving average (ARIMA) modeling might be better suited to capturing the stellar variability component of the light curve (e.g. Haywood et al., 2014; Barclay et al., 2015; Caceres et al., 2015). We did successfully apply a GP regression to take the place of the stellar variability polynomial fit, but we found that any improvements in the fit did not outweigh the larger computational costs.

2.3 Step 3- Quantifying residual autocorrelations

Given the best-fit sinusoidal signal in the data for each of the 1-4 candidate period(s), our goal is to diagnose each signal’s phase and amplitude coherence. This information will be used to select the best of the 1-4 candidate signals and, later in Section 4, to recognize the time-dependent properties of true phase curves.

The phase and amplitude coherence is best examined with the autocorrelation of the residuals of the fit. If the signal experiences large phase and amplitude variations, a sinusoidal fit will be inappropriate, and the residuals will be significantly non-Gaussian. On the other hand, the signal of a true phase curve will have less correlated residuals. We examine two diagnostics of the residual autocorrelation: the Durbin-Watson statistic (Durbin & Watson, 1950) and the Ljung-Box statistic (Ljung & Box, 1978).

The Durbin-Watson test statistic for residuals $\{r_{i}\}$ is given by

[TABLE]

In a Durbin-Watson test, $d$ is compared to critical values at selected significance levels. The test is used to investigate the null hypothesis that the residuals exhibit zero autocorrelation against the alternative hypothesis that the residuals are positively or negatively autocorrelated. Uncorrelated residuals have $d\sim 2$, while positively correlated residuals have $d\ll 2$.

The Ljung-Box statistic is given by

[TABLE]

where $\hat{\rho_{j}}$ is the autocorrelation of the residuals at lag $j$, $N$ is the number of data points, and $h$ is the number of lags being tested. If the residuals are independently distributed, $Q_{h}$ follows a chi-squared distribution with $h$ degrees of freedom, $Q_{h}\sim{\chi_{h}}^{2}$. The residuals show evidence for autocorrelation at significance level $\alpha$ if $Q_{h}>\chi_{1-\alpha,h}^{2}$. Rather than calculate the statistic at one lag, we calculate it at a variety of lags and plot $Q_{h}$ vs. $h$. For example, in Figure 3 we display the $Q_{h}$ vs. $h$ curves for the residuals of the one- and two-component fits of Kepler-13b. (A section of these fits were shown in Figure 2.) $Q_{h}$ exhibits a nearly linear dependence with the lag, with different slopes and y-intercepts between the two fits. Since the “star only” fit has more autocorrelated residuals, the slope and y-intercept of its $Q_{h}$ curve are larger.

Given that the two-component (phase curve & stellar variability) fit described in Section 2.2 will rarely perfectly model a Kepler light curve, there will always be some serial autocorrelation. Therefore to make the Durbin-Watson and Ljung-Box statistics useful, we do not perform the statistical tests directly. Rather, we use the ratio of the test statistics in the two-component fit to their values in the one-component fit. For the Durbin-Watson statistic and Ljung-Box slope and y-intercept, these ratios are denoted $d_{2}/d_{1}$, $m_{2}/m_{1}$, and $b_{2}/b_{1}$, respectively. The ratios for Kepler-13b are displayed in Figures 2 and 3.

As described in Section 2.2, the one-component fit is the “star only” case, that is, the sliding polynomial fit to the data without a sinusoidal component. The ratio of the test statistics in the star & phase curve vs. star only cases thus probes the degree to which the serial autocorrelation is improved or degraded by the inclusion of a sinusoidal signal in the data. Large signal-to-noise phase curves will exhibit a significant reduction in the autocorrelation when the two-component fit is performed, with $d_{2}/d_{1}\gg 1$, $m_{2}/m_{1}\ll 1$, $b_{2}/b_{1}\ll 1$. Similarly, we also examine the degree to which the two-component fit is favored using the ratio of the $\chi^{2}$, denoted by ${\chi^{2}}_{2}/{\chi^{2}}_{1}$. Strongly favored phase curve fits will have ${\chi^{2}}_{2}/{\chi^{2}}_{1}\ll 1$.

Returning to the description of the pipeline, at this stage the autocorrelation diagnostics have been calculated for the best-fit sinusoidal signals associated with each of the 1-4 peak periods. To choose a single best signal, we use the following minimization:

[TABLE]

We thus select the signal with the best fit and lowest autocorrelation. This step concludes the procedure for identifying candidate phase curves in Kepler light curves.

3 Synthetic Phase Curve Injection and Recovery

In Section 2, we built a pipeline for the detection of a candidate phase curve signal in a Kepler light curve. We now construct and inject synthetic phase curves into Kepler light curves and demonstrate the pipeline’s recovery efficiency. These synthetic phase curves will also serve as the training set for the logistic regression in Section 4.

3.1 Phase curve model

Optical phase curve models have been constructed by various authors, including Faigler & Mazeh (2011); Barclay et al. (2012); Placek et al. (2014); Angerhausen et al. (2015) and Esteves et al. (2013, 2015). Here we adopt the notation and model of Esteves et al. (2015). The light curve of a non-transiting planet is composed of three components: the flux coming from the planet (both reflection and thermal), $F_{p}$, and the Doppler beaming and ellipsoidal variation components from the star, $F_{d}$ and $F_{e}$. The composite light curve is given by

[TABLE]

The curve is a function of the phase, $\phi$, which ranges from 0 to 1 and is given by

[TABLE]

where $T_{\mathrm{mid}}$ is the time of inferior conjunction (or mid-transit for a transiting planet).

The planetary brightness is modeled as a Lambert sphere, described by

[TABLE]

where $\cos z=-\sin i\cos(2\pi\left[\phi+\theta\right])$. The planetary brightness is approximately

[TABLE]

and $\theta$ is a phase offset of the peak brightness from the sub-stellar point.

The Doppler beaming flux is given by

[TABLE]

The constant $\alpha_{d}$ is the photon-weighted, bandpass-integrated beaming coefficient.

The ellipsoidal variation component is, using the first three cosine harmonics,

[TABLE]

where

[TABLE]

and $\alpha_{1}$ and $\alpha_{2}$ are functions of the limb darkening and gravity darkening parameters (Esteves et al., 2013).

Though $\alpha_{d}$ is a function of the stellar spectrum and the Kepler transmission function, and $\alpha_{1}$ and $\alpha_{2}$ are complicated functions of limb darkening parameters, they all exhibit nearly linear dependence with $T_{\mathrm{eff}}$ within the range of interest ($4000\ \mathrm{K}<T_{\mathrm{eff}}<7000\ \mathrm{K}$). In Figure 4, we plot the Doppler and ellipsoidal coefficients versus $T_{\mathrm{eff}}$ for 12 host-stars of planets with phase curves from Esteves et al. (2015). We also show best-fit lines. There is a clear linear relationship for each $\alpha$ coefficient with $T_{\mathrm{eff}}$. The trends agree well with an unpictured, additional star with $T_{\mathrm{eff}}=4550\ \mathrm{K}$.

Given the high quality of the fits, we simplify our phase curve model by letting $T_{\mathrm{eff}}$ be a parameter, with the $\alpha$ coefficients calculated using the best-fit lines pictured in Figure 4. Considering equations 6, 8, and 9, there are 9 parameters required to model a phase curve: $P,\ M_{p},\ R_{p},\ i,\ A_{g},\ \theta,\ M_{\star},\ R_{\star},\ T_{\mathrm{eff}}$. We assume that $e\approx 0$ since planets with detectable phase-curves are typically hot Jupiters with low eccentricities.

3.2 Light curve injection

In order to inject simulated phase curves into Kepler datasets, the most thorough procedure would be to insert the signal into pixel-level data, as has been done in measuring the transit signal recovery of the Kepler pipeline (Christiansen et al., 2013, 2015, 2016). However, Christiansen et al. (2013) showed an extremely high fidelity in the preservation of individual injected transit signals from the pixel-level data through image calibration, aperture photometry, and pre-search data conditioning (PDC) systematics removal. Given this precedent, it should be appropriate to inject synthetic phase curves into the PDC light curve data, rather than the pixel-level data.

We picked a random sample of 10,000 Kepler FGK main sequence stars without confirmed planets or KOIs and extracted their PDC light curves. Here we define FGK dwarfs with simple cuts used by Christiansen et al. (2016): $4000\ K<T_{\mathrm{eff}}<7000\ K$, $\mathrm{log}\,g>4.0$. For each target, we generated a single simulated phase curve using the nine physical quantities listed in Table 1. Here $\mathcal{U}$ denotes a uniform distribution. The stellar parameters, $\ M_{\star},\ R_{\star}$, and $T_{\mathrm{eff}}$, are absent from this list. Rather than drawing these parameters randomly, we use the values of the target star. It is important to note that the parameter space spanned by the distributions in Table 1 need not reflect the true distribution of planets; rather, they must encompass all plausible phase curve amplitudes and morphologies capable of being observed.

There are three additional parameters in Table 1 that have not yet been introduced: ${h_{A}}^{2}$, ${h_{\theta}}^{2}$, and $\tau$. These are related to modulations we impose on the planet’s albedo and phase offset. This variability is motivated by the assumption that planets experience surface feature evolution and atmospheric flow that causes their spatially-dependent reflectivity to change slightly over time, which in turn influences the amplitudes and shapes of their observed phase curves. While the magnitude and timescale of this effect is unknown, a fractional albedo variation of up to $\sim 30\%$ seems plausible (Rauscher et al., 2008; Armstrong et al., 2016), and phase offset variations on the order of $\sim 0.15$ have been observed (Armstrong et al., 2016). Moreover, the coherence timescale in the observations of weather on the transiting hot Jupiter, HAT-P-7b is tens to hundreds of days (Armstrong et al., 2016).

To model the albedo and phase offset variability, we use random draws from a Gaussian Process (GP). GPs are used as a non-parametric method of modeling a function in some continuous input space (Gelman et al., 2014). Every realization of a GP is a random variable drawn from a multivariate normal distribution with a mean vector and covariance matrix, where the matrix is constructed from a covariance function that dictates the shrinkage towards the mean and the correlation between pairs of data points. A widely used covariance function is the squared exponential, given by

[TABLE]

In the case of modeling an albedo or phase offset variation, $t$ is time, $h^{2}$ is the maximum covariance of the variation, and $\tau$ is a timescale dictating the smoothness of the variability. $\mathbf{K}$ is the covariance matrix built from the covariance function, given by

[TABLE]

Let $\mathbf{x_{A}}=(x_{A1},x_{A2},...x_{An})$ be a vector representing the fractional deviation of the albedo from its mean at $n$ discrete timesteps. In other words,

[TABLE]

Let $\mathbf{x_{\theta}}$ be the analogous vector for the phase offset. We model $\mathbf{x_{A}}$ and $\mathbf{x_{\theta}}$ as independent multivariate Gaussian random variables: $\mathbf{x_{A}}\sim\mathcal{N}(\mathbf{0},\mathbf{K_{A}})$ and $\mathbf{x_{\theta}}\sim\mathcal{N}(\mathbf{0},\mathbf{K_{\theta}})$. The mean vector is zero because $E[x_{Ai}]=E[x_{\theta i}]=0$.

Figure 5 illustrates several realizations of the GP for $h^{2}=0.0075,\tau=10$ days (top-panel) and $h^{2}=0.015,\tau=20$ days (bottom-panel). In this figure, $\mathbf{x}$ could represent either $\mathbf{x_{A}}$ or $\mathbf{x_{\theta}}$.

We construct a synthetic phase curve for each light curve as follows. We draw a set of parameters from the distributions in Table 1. The covariance matrices, $\mathbf{K_{A}}$ and $\mathbf{K_{\theta}}$, are constructed using the same timescale, $\tau$, but they use their respective covariance amplitudes, ${h_{A}}^{2}$ and ${h_{\theta}}^{2}$. We draw $x_{A}(t)$ and $x_{\theta}(t)$ as GP realizations. For convenience, we make use of the george333http://dan.iel.fm/george/current/ code for Gaussian Process regression (Ambikasaran et al., 2014). Given the random draws for $\left<A_{g}\right>$ and $\left<\theta\right>$, the albedo variation function is calculated as $A_{g}(t)=\left<A_{g}\right>[1+x_{A}(t)]$ and the phase offset $\theta(t)=\left<\theta\right>[1+x_{\theta}(t)]$. These time-dependent quantities and the rest of the parameters are then inserted into the phase curve model outlined in Section 3.1. Finally, the model phase curve is multiplied by the Kepler PDC flux to create the synthetic light curve.

3.3 Recovery efficiency

We now wish to evaluate the efficiency with which these synthetic phase curves are recovered by the pipeline outlined in Section 2. We ran the pipeline on each of the 10,000 synthetic injections. We define the synthetic phase curve to be “recovered” when the period of the detected signal is within 0.01 days of the injected signal. In the upper panel of Figure 6, we plot the histograms of the phase curve semi-amplitudes for all injections and for recovered injections only. The shape of the histogram for the injected data has little physical significance other than to reflect the typical amplitude of the injected phase curves. We show the fraction of recovered injections in the lower panel of the figure. The recovery efficiency is a smooth and increasing function of amplitude, with 20 ppm amplitude phase curves being recovered $\sim 50\%$ of the time.

For the injections that were recovered, we plot in Figure 7 the semi-amplitude of the detected signal versus that of the injected signal. The relationship makes clear that the pipeline attains a reasonably accurate recovery of the phase curve amplitudes. The bias towards underestimation at large amplitude is due to the fact that the pipeline uses sinusoidal fits that average over ellipsoidal features in the light curve.

4 Binary Classification Using Logistic Regression

Having demonstrated that the pipeline is capable of retrieving the signals of synthetic injected phase curves, we proceed to classify light curves by employing a logistic regression algorithm in a supervised machine learning context. We aim to have the algorithm “learn” the properties of planetary phase curves and use this information to determine whether phase curves are present in newly introduced light curves. We offer an introduction to logistic regression in Section 4.1 and apply it to phase curve classification in Section 4.2.

4.1 Logistic regression overview

Logistic regression is a technique used for prediction of a categorical dependent variable based on observed characteristics called explanatory variables or predictors (Gelman et al., 2014). The dependent variable is typically a binary response (e.g. pass/fail, healthy/sick). Like other forms of regression, there may be any number of predictors, and they can be either continuous or discrete.

Consider a dataset of length $n$ containing binary response dependent variables, $Y_{i}$, which may take on the values 0 or 1. Assume each response variable is associated with $m$ explanatory variables: $x_{1i},x_{2i},...,x_{mi}$. The subscript $i$ indicates an index corresponding to each outcome in the dataset. Let the (unknown) probability of success of the outcome be $p_{i}$, where $p_{i}$ is unique to each outcome, but is related to the explanatory variables:

[TABLE]

The outcomes are thus Bernoulli random variables,

[TABLE]

with a probability mass function,

[TABLE]

In a logistic regression model, we assume that the probability of success varies systematically as a function of the explanatory variables:

[TABLE]

Here, $\boldsymbol{\beta}$ is the vector of regression coefficients,

[TABLE]

and $\mathbf{x_{i}}$ is the vector of explanatory variables,

[TABLE]

A value of 1 has been appended as the first element corresponding to the intercept coefficient, $\beta_{0}$. Note that when a given regression coefficient, $\beta_{j}$, is zero, the outcome is independent of the corresponding explanatory variable, $x_{ji}$.

Rewriting equation 17 in terms of $p_{i}$ yields

[TABLE]

Using equations 16 and 20, one can show that the log-likelihood $l(\boldsymbol{\beta})$ can be expressed as

[TABLE]

The regression is solved via maximum likelihood estimation (MLE) by identifying the vector of regression coefficients, $\boldsymbol{\hat{\beta}}$, that maximize $l(\boldsymbol{\beta})$. In a Bayesian context, one can impose a prior distribution on the regression coefficients, $\pi(\boldsymbol{\beta})$. Once the vector of coefficient estimators is found, the unknown probabilities, $p_{i}$, may be calculated with equation 20. More importantly, one may estimate the probabilities of unknown outcomes if provided the observations of the relevant explanatory variables.

In the analysis that follows, we perform logistic regression using the Scikit-learn library444http://scikit-learn.org/ (Pedregosa et al., 2011). We use a regularization strength, $\alpha=1$, and we use Newton-CG optimization for the determination of the regression coefficients.

4.2 Phase curve classification using logistic regression

We now apply the logistic regression model to the problem of classifying phase curve signals in Kepler light curves. The response variable, $Y_{i}$, is whether the detected signal is a phase curve ($Y_{i}=1$) or not ($Y_{i}=0$). The probability, $p_{i}$, is the probability that a given detection is a true planetary phase curve. (Note this is different from the probability that a phase curve is present in the light curve.)

We used a collection of $m=12$ predictors, many of them introduced in Section 2:

Signal period, $P$ 2. 2.

Signal amplitude 3. 3.

Chi-square ratio of the two-component to one-component phase curve fits, ${\chi^{2}}_{2}/{\chi^{2}}_{1}$ 4. 4.

Durbin-Watson statistic ratio, $d_{2}/d_{1}$ 5. 5.

$Q_{h}$ vs. $h$ slope ratio, $m_{2}/m_{1}$ 6. 6.

$Q_{h}$ vs. $h$ y-intercept ratio, $b_{2}/b_{1}$ 7. 7.

Full-width half maximum (FWHM) of the peak at $P$ in the Lomb-Scargle periodogram 8. 8.

Significance of the peak in the Lomb-Scargle periodogram 9. 9.

Significance of the peak in the “local significance” periodogram 10. 10.

Normalized significance of the $P/2$ harmonic peak in the Lomb-Scargle periodogram 11. 11.

Normalized significance of additional peaks in the periodogram 12. 12.

Deviation of the phase-folded light curve from sinusoidality

Predictors 1-9 were discussed in Section 2. Here we introduce predictors 10-12, which were all motivated by inspection of the signals in the synthetic datasets.

Predictor 10: True phase curves should exhibit some power in the periodogram at the $P/2$ harmonic due to the ellipsoidal variation component of the phase curve. We quantify this with predictor 10 by finding the significance of the $P/2$ harmonic peak in the Lomb-Scargle periodogram and taking its ratio with respect to the peak at $P$.

• -

Predictor 11: For true phase curve detections, there should not be significant peaks remaining in the periodogram once the signals at $P$ and $P/2$ are removed. For non-phase curve signals such as astroseismic pulsations (e.g. Grigahcène et al., 2010), significant peaks are accompanied by additional peaks. To calculate predictor 11, we start with the Lomb-Scargle periodogram and remove the peaks at the signal period, $P$, and at the harmonic, $P/2$. We then calculate the local significance periodogram, find the highest peak, and normalize it with respect to the height of the peak at $P$.

• -

Predictor 12: Large-amplitude phase curves should tend to exhibit a deviation from sinusoidality due to the presence of the beaming and ellipsoidal components. We quantify this tendency with predictor 12. We first use the fit in Section 2.2 to detrend the light curve. We phase-fold it at the candidate period, bin it, fit a sinusoid to the phase-folded signal, and calculate the residuals of the fit. The predictor is then the Durbin-Watson statistic of the residuals (see Section 2.3). This quantifies the autocorrelation in the residuals, or the deviation from sinusoidality.

All 12 predictors were shown to exhibit significant differences between phase curve detections and non-phase curve detections in the synthetic datasets, either independently or in combination with other predictors. We show an example of this using a triangle diagram in Figure 8. This figure shows the correlations between four of the predictors: predictor 3, ${\chi^{2}}_{2}/{\chi^{2}}_{1}$; predictor 4, $d_{2}/d_{1}$; predictor 5, $m_{2}/m_{1}$; and predictor 11, $x_{11}$. Plotted in blue are the recovered phase curve injections in the synthetic data. In green are the detections that did not correspond to the injected phase curve.

It is clear to see that there are significant differences in the two populations. Consider, for example, the top left plot of ${\chi^{2}}_{2}/{\chi^{2}}_{1}$ vs. $d_{2}/d_{1}$. Both phase curves and non-phase curves exhibit ${\chi^{2}}_{2}/{\chi^{2}}_{1}<1$ because the two-component fit from section 2.3 has more degrees of freedom than the one-component fit. However $d_{2}/d_{1}$ for phase curves tends to be larger and follow a straight line because the two-component fit does the best job at removing the residual autocorrelations when the signal is a true phase curve, or a near-sinusoid. (Recall that as the amount of positive autocorrelation decreases, $d$ increases and $m$ decreases. Therefore, for strong phase curve detections, we should expect $d_{2}/d_{1}\gg 1$ and $m_{2}/m_{1}\ll 1$.) Generally speaking, the separation between the phase curve and non-phase curve populations continues at smaller scales (where the points are overlapping), but it is less pronounced.

Given the distinct differences in the predictors between signals that are or are not phase curves, the multivariate logistic regression should be capable of picking up on these differences and classifying light curves accordingly. We describe this process next.

4.3 Testing on synthetic phase curves

We applied the logistic regression to the set of 10,000 synthetic phase curves that we produced in Section 3.2. We first randomly split the synthetic dataset into a training set and a test set, with the test set being 5% of the total. We then regressed on the training set using the 12 predictors introduced in the previous section and the outcomes ${Y_{i}}$ indicating whether each detected signal of the training group was the injected phase curve or not. Finally, we ran predictions on the test set, supplying the 12 predictors and calculating the predicted probabilities that each member of the test set was or was not a phase curve. The predictions were then compared to the actual results, that is, the outcomes of whether each signal detected by the pipeline was the injected phase curve or not.

The results of the test are shown in Figure 9. On the x-axis is the semi-amplitude of the signal derived by the two-component fit described in Section 2.2. (Note this amplitude is half the peak-to-peak amplitude.) The y-axis is the probability predicted by the logistic regression that the detection is a true phase curve. Black dots are correct predictions, and red dots are incorrect predictions. More explicitly, black dots with $\mathrm{P(phase\ curve)}>0.5$ are cases where the logistic regression predicted that the signal detected by the pipeline was a phase curve, and the signal did indeed correspond to the injected phase curve. Red dots with $\mathrm{P(phase\ curve)}>0.5$ are cases where the logistic regression predicted that the signal was a phase curve, but the signal did not correspond to the injected phase curve. The situation is reversed for points below the $\mathrm{P(phase\ curve)}=0.5$ line.

Generally speaking, the logistic regression performs remarkably well. In total, correct predictions are made $\sim 90\%$ of the time, and the accuracy increases steeply as $\mathrm{P(phase\ curve)}$ approaches 0 or 1. Larger amplitude detections are assigned more definitive probabilities (i.e. $\mathrm{P(phase\ curve)}$ closer to 0 or 1).

We tested the consistency of the algorithm’s predictive performance using 5,000 realizations of the training and test procedure described in the first paragraph of this section. Each realization used a random partition of the synthetic dataset into the 95% training and 5% test sets. For each trial, we performed the logistic regression and identified the correct and incorrect predictions. Histograms of the fraction of correct predictions are shown in Figure 10. The green histogram shows that the algorithm predicts correctly $\sim 91\%$ of the time. For cases with reported probabilities greater than 0.9 (purple histogram), the prediction is correct $\sim 98\%$ of the time. As indicated by the widths of the histograms, these results are fairly consistent from trial to trial.

We also computed the fraction of correct predictions as a function of the probability reported by the logistic regression. The mean curve over the 5,000 random realizations is displayed in Figure 11. The results indicate that the reported probabilities can indeed be interpreted as the likelihood that a given signal corresponds to a true phase curve.

4.4 Testing on transiting Kepler hot Jupiters

We also tested the detection pipeline and logistic regression on the collection of transiting hot Jupiters in the Kepler field. We first compiled the list of KOIs in the Q1-Q17 DR 24 catalog (Coughlin et al., 2016) from the NASA Exoplanet Archive555http://exoplanetarchive.ipac.caltech.edu (Akeson et al., 2013). We filtered the confirmed planets and planet candidates to those with periods in a 1-5 day range and radii in the range 0.5-3 $R_{\mathrm{Jup}}$. For the unconfirmed planet candidates, we used the KOI false positive probabilities (FPPs) calculated by Morton et al. (2016) to restrict the list to those with a FPP less than 10%. These steps resulted in a collection of 59 total hot Jupiter planets and planet candidates with potentially detectable phase curves.

The pipeline described in Section 2 was executed in the same manner as on the synthetic datasets. The only difference was that the transits and (potential) secondary eclipses were removed from the light curves before processing them. The pipeline recovered planetary phase curves in four cases: HAT-P-7b (also called Kepler-2b), Kepler-13b, Kepler-41b, and Kepler-76b. In these cases, the period of the signal recovered by the pipeline matched the planet period to within 0.01 days. This apparently low recovery rate is not surprising. Esteves et al. (2015) conducted a comprehensive search for planetary phase variations in Kepler transiting planets, finding 14 planets total. Twelve of these have parameters fitting our criteria ($1\ \mathrm{day}<P<5\ \mathrm{days}$, $0.5\ R_{\mathrm{Jup}}<R_{\mathrm{P}}<3\ R_{\mathrm{Jup}}$). In this sense, there are 12 potentially detectable phase curves, and one-third of these have amplitudes large enough to be detected by our pipeline.

Recovery by the pipeline is distinct from recovery by the logistic regression. The next step is to consider predictions of the transiting dataset to confirm that the logistic regression places high probabilities on the four recovered phase curve signals and low probabilities on the rest.

Using the regression that was performed on the synthetic data in the previous section, we made predictions on the light curves of the transiting planets. The results are shown in Figure 12. In all four cases where the pipeline recovered the planet’s phase curve, the logistic regression correctly predicted that the detection was a true phase curve. This is represented by the four annotated black points with $\mathrm{P(phase\ curve)}>0.5$. In all but one of the remaining cases, the logistic regression correctly predicted that the signal detected by the pipeline was not a phase curve. Put another way, all of the black points with $\mathrm{P(phase\ curve)}<0.5$ were cases where the pipeline did not recover the planet’s period, and the logistic regression correctly reported that the signals were not real phase curves.

The one false positive pictured in the upper right of Figure 12 has an amplitude greater than that expected from a planetary phase curve. The signal is from the light curve of KIC 5651104/KOI-840/Kepler-695b. UKIRT images available from the Kepler Exoplanet Follow-up Observing Program (ExoFOP)666https://exofop.ipac.caltech.edu/ show that the star has nearby background or foreground companions, so the signal probably results from blended, non-eclipsing binary stars. We found further evidence for this interpretation in that the phase curve amplitude is inconsistent when the light curve is folded at two times the detected period. Details on the nature of this inconsistency and steps taken to rule out these types of false positives will be discussed in Section 5.1.

In short, the logistic regression performs exceedingly well on both synthetically generated datasets and the light curves of transiting Kepler planets. In the vast majority of cases, the algorithm correctly predicts when a signal is or is not a true planetary phase curve. This gives us confidence in the process of applying the logistic regression to novel light curves and inspecting them for phase curve detections.

5 Application to Kepler FGK stars

We now describe the application of our pipeline and logistic regression algorithm to the search for non-transiting planets around Kepler FGK stars without confirmed planets or KOIs. We first downloaded the Q1-Q17 DR 25 Kepler stellar catalog (Mathur et al., 2017) from the NASA Exoplanet Archive${}^{\ref{NASA Exoplanet Archive}}$ (Akeson et al., 2013). We filtered the target stars for FGK main sequence stars, using the criteria from Christiansen et al. (2016): $4000\ K<T_{\mathrm{eff}}<7000\ K$ and $\mathrm{log}\,g>4.0$. We then removed stars with known planets and planetary candidates, leaving 146,980 targets. For these targets, we retrieved all 17 quarters of Kepler long-cadence, pre-search conditioning (PDC) photometry. These light curves are publicly available at the Mikulski Archive for Space Telescopes (MAST).

We then ran the pipeline outlined in Section 2 on the set of light curves, with successful convergence on 142,630 of them. This procedure identified the most likely phase curve signal in each light curve, if a phase curve were to be present. Using the logistic regression that was performed on the synthetic data in Section 4.3, we then made predictions on this novel set of light curves, evaluating the likelihood that each detected signal could correspond to a true phase curve.

The results are displayed in Figure 13. This plot is analogous to Figures 9 and 12, with the x-axis being the semi-amplitude of the signal derived by the two-component fit described in Section 2.2 and the y-axis being the probability output by the logistic regression described in Section 4.2. The points have been colored by their density, and the panel at the right is a histogram of the probabilities on a log scale. As expected and required, the vast majority of targets have $P(\mathrm{phase\ curve})<0.5$.

5.1 Binary star filtering

Among the potential sources for astrophysical false positives, short period binary stars are one of the most significant concerns. Binary stars have an occurrence rate of $\gtrsim 10\%$ at periods of 1-5 days (Kirk et al., 2016), about two orders of magnitude greater than the occurrence rate of hot Jupiters with phase variations detectable by our pipeline. Their phase curves are dominated by ellipsoidal variations (Faigler et al., 2012), meaning that our pipeline would detect them at half their orbital period. Although the amplitudes of their ellipsoidal variations are typically 100-1000 ppm, the amplitude is given by $A_{\mathrm{ellip}}\propto\sin^{2}i$ (see equation 9), so inclined systems could appear with amplitudes less than 100 ppm. Moreover, potential background or foreground stars could contaminate the photometric aperture and make the ellipsoidal variations appear with small amplitudes.

Fortunately, there are means of discriminating between a binary system and a giant planet phase curve, even at small amplitudes. Because ellipsoidal variations dominate the binary stars’ light curve, they would be detected at half their orbital period. However, the Doppler beaming amplitude is non-negligible (Faigler et al., 2012), meaning that the periodogram of a binary light curve should show significant power at the true period, or twice the detected period. Planetary phase curves should not show any signal at $2P$. Furthermore, when the light curve is folded at twice the detected period and split at the midpoint, the beaming amplitude of a binary system will cause the two halves to be distinct, but a planetary signal will have two identical curves. To address these features and eliminate binary stars, we constructed the following three filters. The values of the thresholds were determined by inspection of the synthetic datasets.

Require $\mathrm{power}_{2P}/\mathrm{power}_{P}<0.4$, where $P$ is the detected period (or, in the case of binaries, half the true period). 2. 2.

If the light curve is folded and binned at $2P$, and if $y_{i}$ represents the data in the first half and $z_{i}$ represents the second half, then require $0.8<r<1.2$, where $r=\frac{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-z_{i})^{2}}{{\sigma_{y}}^{2}+{\sigma_{z}}^{2}}$, ${\sigma_{y}}^{2}=\sum_{i=1}^{N}{\sigma_{yi}}^{2}$ and ${\sigma_{z}}^{2}=\sum_{i=1}^{N}{\sigma_{zi}}^{2}$. We take $N=200$, such that 400 points span the $2P$ folded light curve. 3. 3.

Reqiure $\frac{\mathrm{Amp}}{300\ \mathrm{ppm}}<\left(\frac{P}{1\ \mathrm{day}}\right)^{-4/3}$, where Amp is the semi-amplitude.

To apply these filters to the sample, we eliminated all points with $P>0.5$ that did not adhere to these limits. The results of these filters are shown in Figure 16, where they are discussed in conjunction with the additional vetting in the next section.

5.2 Further vetting via physical plausibility considerations

In addition to binary stars, there are likely additional types of false positives that are not being taken into account (Shporer, 2017). It is therefore worthwhile to take further steps to address these.

Another powerful method for false positive vetting is to consider whether the signal morphologies are physically plausible given the expectations for planetary phase curves. Similar to the binary star filtering in the previous section, the additional vetting that we will now describe was only applied to cases with predicted probabilities, $P>0.5$ (hereafter called “candidates”).

We started with a fitting procedure for assessing the physical plausibility of the candidate phase curves. For each candidate, we calculated a one-component sinusoidal fit, representing a fit for reflected light only. We then calculated a three-component fit modeled by

[TABLE]

This is a simple way to assess the phase curve without applying the full phase curve model described in Section 3.1. The $A_{\mathrm{refl}}$ term represents the reflected light component with one phase, $\theta_{\mathrm{refl}}$, and the $A_{\mathrm{ellip}}$ and $A_{\mathrm{beam}}$ terms are the ellipsoidal and beaming components with a separate phase, $\theta_{\mathrm{ellip}}$. The thermal component is incorporated into the $A_{\mathrm{refl}}$ term, since the two cannot be separated without more advanced modeling. An example of the one-component and three-component phase curve fits for one of the candidates, KIC 5716330, is shown in Figure 14.

We then used these light curve fits to examine the candidates for false positives. We first calculated the Durbin-Watson statistic, $d$, of the residuals of the three-component fits (see equation 1). If the statistic is far from 2, this indicates autocorrelated residuals and implies that the light curve is not well-described by the fit. We eliminated all candidates with $d<1.85$ or $d>2.15$, where this range was determined by inspection of the synthetic phase curves.

Next, we examined the candidate phase curves’ physical plausibility by considering $A_{\mathrm{refl}}$ and $A_{\mathrm{ellip}}$. If, for example, $A_{\mathrm{refl}}$ is too large, the planet would have an abnormally large radius or albedo. If $A_{\mathrm{ellip}}/A_{\mathrm{refl}}$ is too large, the density would be abnormally high.

To convert $A_{\mathrm{refl}}$ and $A_{\mathrm{ellip}}$ into physical quantities, we first solve for the planet’s radius in terms of $A_{\mathrm{refl}}$ and other observables:

[TABLE]

The factor of 2 in front of $A_{\mathrm{refl}}$ comes from the fact that $A_{\mathrm{refl}}$ is a semi-amplitude but equation 7 considers a peak-to-peak amplitude.

Next, we solve for the planet’s mean density in terms of $A_{\mathrm{ellip}}/{A_{\mathrm{refl}}}^{3/2}$:

[TABLE]

Here $\alpha_{2}$ is a constant in the ellipsoidal variation model (see equation 9). It varies nearly linearly with $T_{\mathrm{eff}}$ (see Figure 4).

For each candidate, we calculated estimates of $R_{p}$ and $\rho_{p}$ by using the Kepler stellar catalog properties and by assuming $A_{g}\in[0.1-0.3]$ and $i\in[15^{\circ},85^{\circ}]$. This domain thus corresponds to an allowable region in $\rho_{p}-R_{p}$ space. The candidate’s plausibility may be considered by comparing this region to the radii and density of known hot Jupiters. In Figure 15, we plot $\ln\rho_{p}$ vs. $R_{p}$ for all hot Jupiters with measured masses and radii, as obtained by the Exoplanets Data Explorer777exoplanets.org (Han et al., 2014). (We considered planets with $P<5\ \mathrm{days}$ and $0.2\ M_{\mathrm{Jup}}<M_{p}\sin i<5\ M_{\mathrm{Jup}}$.) We also plot the $\rho_{p}-R_{p}$ regions for two candidates, one with strong physical plausibility and one with weak plausibility.

A quantitative metric is necessary for comparing a given candidate’s allowed $\rho_{p}-R_{p}$ domain to the region occupied by known hot Jupiters. We first calculated a kernel density estimation of the values for known hot Jupiters in $\ln\rho_{p}$-$R_{p}$ space, denoted $f(R_{p},\ln\rho_{p})$ and illustrated with the gray gradient in Figure 15. We then calculated the integral, $S$, of the kernel density over each candidate’s allowed domain,

[TABLE]

where

[TABLE]

We calculated the $S$ integrals for all candidates with probabilities, $P>0.5$. We then removed all candidates with $S<5$. Though the choice of this threshold is somewhat arbitrary, it corresponds to $\gtrsim 50\%$ of a candidate’s $\rho_{p}$-$R_{p}$ domain overlapping with that of the known planet detections.888It is important to note that this is a quite conservative vetting procedure, since the albedo and inclination ranges are not as wide as possible, and since there could feasibly be new planet detections on the outskirts of the known domain. However, in the process of selecting the best possible candidates, conservative vetting is a logical solution. We also verified that the four recovered transiting hot Jupiters from Section 4.4 had $S$ integral measures greater than this threshold.

Figure 16 shows the results of the binary star filtering (Section 5.1) and the physical plausibility vetting (Section 5.2) on the aggregate of candidates with $P>0.5$. These filters clearly have a dramatic effect on the number of plausible planetary phase curve candidates.

6 Results: A Catalog of Candidate Planets with Phase Curves

In the remainder of the paper, we focus our analysis on the subset of highest probability candidates. Figures 9 and 12 show that true phase curves should be reported with probabilities very close to 1. To this end, we selected all candidates with $P>0.97$ and $\mathrm{Amp}<70$ ppm. Lastly, we established a final candidate list through a small amount of visual inspection. In particular, we inspected the Lomb-Scargle periodograms and removed any that showed signs of significant peaks apart from the primary peak. We also examined the $2P$-folded phase curves and checked for any subtle indications of amplitude inconsistency in each half.

Following this final vetting procedure, we present 60 remaining high probability candidates. In Table 5, we list their parameters. The parameters $K_{p}$, $T_{\mathrm{eff}}$, $R_{\star}$, $M_{\star}$, and Fe/H are from the Q1-Q17 DR 25 Kepler stellar catalog (Mathur et al., 2017). The uncertainties on the periods were calculated using a Gaussian fit to the peak in the Lomb-Scargle periodogram. The units of $A_{\mathrm{refl}}$ and $A_{\mathrm{ellip}}$ are ppm.

We also report estimates of the candidates’ minimum RV semi-amplitudes, $K\sin i$. To calculate this, we note that the amplitude of the ellipsoidal variation component, $A_{\mathrm{ellip}}$, obtained via equation 9, is directly related to $K$ through the relation,

[TABLE]

The errorbars on $K\sin i$ were obtained via propagation of the uncertainties on the KIC stellar parameters and the $A_{\mathrm{ellip}}$ estimates obtained from the least-squares fit introduced in the first paragraph of Section 5.2. Unfortunately these uncertainties are quite large, mostly due to the lack of constraint on the stellar parameters.

An online repository of the candidates may be found at https://smillholland.github.io/Non-transiting_HJs/. The repository contains a downloadable candidate catalog and a variety of diagnostic plots, including the candidate phase curves, periodograms, and $\rho_{p}-R_{p}$ diagrams. We also provide information relevant to follow-up observations, including estimates of the orbital ephemerides.

6.1 Albedo constraints

Under the assumption that the majority of remaining high probability candidates correspond to true planets, we are presented with a collection of phase curves that may be examined in greater detail. One quantity of interest is the planetary albedo. Without additional information, it is impossible to disentangle the degeneracy between the albedo, planet radius, and inclination. However, we can consider the composite quantity, $A_{g}{R_{p}}^{2}\sin i$, which can be calculated directly from the observable quantities as follows:

[TABLE]

Here, $P$ and $A_{\mathrm{refl}}$ are the period and reflection semi-amplitude of the candidate phase curve, and $M_{\star}$ is taken from the Kepler stellar catalog.

In Figure 17, we display the results of the $A_{g}{R_{p}}^{2}\sin i$ calculations. In the top panel, we show the histogram of $A_{g}{R_{p}}^{2}\sin i$ for the candidates in blue. We also show the histograms of all synthetic injections and the recovered synthetic injections from Section 3. The ratio of these two gives a recovery efficiency curve, and in the bottom panel, we use this efficiency curve to construct the bias-corrected histogram for the candidates. We compare these results to 12 known transiting hot Jupiters with phase curves from Esteves et al. (2015). It is clear from the comparison of the bias-corrected candidate CDF to the transiting hot Jupiter CDF that the candidates are a good match to known planets. The candidate CDF is slightly shifted towards lower values, which is an expected outcome of the bias correction. Furthermore, the bias-corrected histogram peaks at $A_{g}{R_{p}}^{2}\sin i\sim 0.1$, providing further evidence that hot Jupiters are generally fairly dark.

It is interesting to check whether the candidates show any systematic trends in albedo with other quantities. In Figure 18, we plot $A_{g}{R_{p}}^{2}\sin i$ vs. $T_{\mathrm{irr}}$, where $T_{\mathrm{irr}}$ is defined as

[TABLE]

such that the incident stellar irradiation is $F_{0}=\sigma T_{\mathrm{irr}}^{4}$ (Perna et al., 2012). $T_{\mathrm{irr}}$ is not the planet’s equilibrium temperature, which is rather given by

[TABLE]

where $A_{B}$ is the Bond albedo and $f$ is a re-radiation factor ranging between 1/4 (homogeneous redistribution) and 2/3 (instant re-radiation) (Esteves et al., 2013). However, $T_{\mathrm{irr}}$ is a close proxy to $T_{\mathrm{eq}}$.

The gray gradient in Figure 18 illustrates regions with stronger detection bias. This was calculated by finding the point density of recovered synthetic injections and dividing it by the point density of all synthetic injections.

There is little obvious trend in $A_{g}{R_{p}}^{2}\sin i$ vs. $T_{\mathrm{irr}}$. There is a weak indication that the average $A_{g}{R_{p}}^{2}\sin i$ decreases as a function of $T_{\mathrm{irr}}$ up to $T_{\mathrm{irr}}\sim 2700\ \mathrm{K}$, at which point the phase curve would start involving substantial thermal radiation, and $A_{g}$ would not just include reflected light. This may suggest that cooler planets are brighter on average, with a greater proportion of reflective clouds. It is not surprising to see the lack of a strong relationship, however. We are using Kepler stellar catalog parameters, so $T_{\mathrm{irr}}$ is not known to high accuracy. It would be valuable to revisit this after more precise stellar parameters have been obtained. Data from the Gaia spacecraft (Gaia Collaboration et al., 2016) may be a near-term resource for these improved estimates.

6.2 Offsets of the phase curve maxima

A phase offset between the phase curve maximum and sub-stellar point has been observed for most planetary phase curves. The maximum is shifted eastward for thermal phase curves (Showman & Guillot, 2002; Knutson et al., 2007), as superrotating jets advect the hottest region downwind. Both eastward and westward shifts have been observed for optical phase curves (Shporer & Hu, 2015; Esteves et al., 2015).

For phase curves of transiting planets, the phase offset is easily determined by comparing the phase curve maximum to the secondary eclipse. Determining $\Delta\theta$ is not easy for non-transiting planets, but it is possible in some cases. If beaming and/or ellipsoidal components of the light curve are detected, the phase offset can be obtained from the phase relationship between the reflected light and beaming/ellipsoidal components.

Towards this goal, we used the three-component phase curve fits (equation 22) of the candidates from Section 5.2. We calculated the phase offset between the reflection maximum and sub-stellar point using $\Delta\theta=\theta_{\mathrm{refl}}-\theta_{\mathrm{ellip}}$. For the candidate pictured in Figure 14, the phase offset is consistent with an eastward shift of the phase curve maximum.

For each of the candidate phase curve fits, we also calculated the differences in the Aikaike Information Criterion (AIC) between the one-component and three-component fits. If the three-component fit is favored, then $\Delta\mathrm{AIC}$ is large and the ellipsoidal/beaming components are more significant. In Figure 19, we plot $\Delta\theta$ versus $\Delta\mathrm{AIC}$ and $T_{\mathrm{irr}}$. $\Delta\theta>0$ corresponds to eastward shifts in the phase curve maxima (i.e. the phase curve peaks before the substellar point), while $\Delta\theta<0$ corresponds to westward shifts.

From the left panel in Figure 19, we see that candidates with smaller $\Delta\mathrm{AIC}$ favor westward shifts in the phase curve maxima, whereas the maxima of candidates with larger $\Delta\mathrm{AIC}$ tend to be eastward shifted. Moreover, the right panel provides weak evidence for a positive correlation of $\Delta\theta$ with $T_{\mathrm{irr}}$. Both of these correlations strengthen appreciably when considering candidates with higher probability thresholds (e.g. $P>0.985$, rather than the 0.97 used here). We emphasize again that we are using Kepler stellar catalog parameters in the calculation of $T_{\mathrm{irr}}$, which likely explains much of the scatter in Figure 19. It will be valuable to see how the results improve with more reliable spectroscopic and astroseismic stellar parameter estimates.

It is logical that $\Delta\mathrm{AIC}$ and $T_{\mathrm{irr}}$ should show a similarly signed correlation with $\Delta\theta$, as $\Delta\mathrm{AIC}$ is closely related to $T_{\mathrm{irr}}$. As $\Delta\mathrm{AIC}$ increases, the ellipsoidal component of the light curve becomes more significant. Among other dependencies, this correlates with a small $a$ or large $R_{\star}$ (see equation 9), both of which would act to increase $T_{\mathrm{irr}}$ (equation 29).

Several previous studies have also observed a positive trend in $\Delta\theta$ vs. $T_{\mathrm{irr}}$ or $T_{\mathrm{eq}}$ (Angerhausen et al., 2015; Shporer & Hu, 2015; Esteves et al., 2015; Shporer, 2017). Optical phase curves contain contributions from both reflected light and thermal radiation. For the hottest planets, the phase curve contains a significant thermal contribution, and the eastward phase shift matches that which is observed in IR phase curves. This phase shift is thought to be due to a superrotating equatorial jet advecting the hot spot eastward of the substellar point (e.g. Showman & Guillot, 2002). For the coolest planets, the phase curve is reflection-dominated, and the westward phase shift has been attributed to the presence of reflective clouds condensing in the cooler regions of the planet, westward of the substellar point (e.g. Shporer, 2017). The linear correlation in $\Delta\theta$ vs. $T_{\mathrm{irr}}$ could then just arise from the relative contribution of each phase curve component.

The presence of both thermal and reflection components in the phase curve makes it challenging to draw physical interpretations of the atmospheric dynamics. Considering the hotspot alone, an increase in $T_{\mathrm{irr}}$ should result in a decrease in the magnitude of the phase offset, as heat redistribution becomes less efficient (Perna et al., 2012). This arises from the balance between radiative and advective timescales, $t_{\mathrm{rad}}$ and $t_{\mathrm{adv}}$. With

[TABLE]

and

[TABLE]

$t_{\mathrm{adv}}/t_{\mathrm{rad}}$ increases rapidly with $T_{\mathrm{irr}}$, such that the heat redistribution efficiency and the hotspot offset both decrease with $T_{\mathrm{irr}}$ (Perna et al., 2012).

Ohmic dissipation, which becomes significant at large $T_{\mathrm{irr}}$ due to an increase in thermal ionization and the resulting electrical conductivity, should also act to slow the zonal winds and reduce the magnitude of the hotspot offset (Perna et al., 2010a, b; Batygin et al., 2013). This effect works constructively with the balance between $t_{\mathrm{adv}}$ and $t_{\mathrm{rad}}$ described above.

In order to see these theorized trends of $\Delta\theta$ vs. $T_{\mathrm{irr}}$ in the data, it would be necessary to isolate the thermal component from the reflection component of the phase curve (e.g. Armstrong et al., 2016). This would be a valuable undertaking once the candidates presented in this paper are confirmed and the planets’ irradiation temperatures may be determined more accurately.

7 Discussion

7.1 Caveats & Future Work

There is an important caveat related to the identification and characterization of non-transiting planet candidates presented in this paper. The concern relates to the effects of nearby stellar companions. Using high-resolution optical and near-IR imaging, Furlan et al. (2017) found that $\sim 30\%$ of KOI host stars have at least one companion star within 4″. Since the Kepler detector has $\sim$4″wide pixels and the aperture photometry is obtained from at least a few pixels, this implies that $\gtrsim 30\%$ of the candidates presented here will have signal amplitudes that are smaller than they would otherwise be without nearby companions.

In addition to influencing the constraints on the candidates’ albedo, planet radius, and inclination (see Section 6.1), this photometric contamination can also result in false positive phase curve detections in the form of background or foreground binary star phase curves, analogous to the dominant transit false positive sources (e.g. Morton et al., 2016). Although the binary star filtering and physical plausibility vetting outlined in Sections 5.1 and 5.2 have certainly helped reduce some of these false positives, we do not claim to have fully accounted for the problem.

It might be possible to further reduce the effects of background or foreground binary stars by including them in the synthetic dataset that was first produced in Section 3.2. To do this, one could first assume a distribution of angular separations and magnitude differences of nearby stellar companions using the results of Furlan et al. (2017). Then, the distribution of short-period binary stars could be empirically modeled after the findings of Kirk et al. (2016) and other papers in this series. After injecting signals appropriate to these distributions, one could search for predictors to add to the logistic regression that would help distinguish the binary signals from the planetary phase curve signals.

7.2 Candidate follow-up observations

The main result of this paper is a catalog of 60 non-transiting, short-period giant planet candidates with phase curves (presented in Section 6). These candidates can be confirmed with follow-up radial velocity observations. Even though the target stars are dim, hot Jupiters are readily detectable even with low RV precision. Moreover, the known orbital ephemerides from the phase curves will significantly aid the RV measurements. Once the RVs confirm a planet candidate’s existence, the comparison between the RV and photometric ephemerides will allow for a more precise measurement of the shift between the phase curve maximum and sub-stellar point. Moreover, the combined RV-photometric modeling could permit a break in the mass/inclination degeneracy.

After RV confirmation, other valuable follow-up observations include Spitzer infrared photometry of the brightest targets. There are very few planets with both optical and infrared phase curve observations (Christiansen et al., 2010; Demory et al., 2013). The combination of these observations, and particularly a comparison of the offset of the phase curve maxima, would be strongly constraining on the planets’ atmospheric dynamics.

7.3 Candidate stellar metallicities and effective temperatures

It is well-known that the occurrence rate of giant planets correlates strongly with stellar metallicity (Santos et al., 2004; Fischer & Valenti, 2005; Johnson et al., 2010). Consideration of the stellar metallicity of the highest probability candidates thus provides a strong independent check on their validity; the distribution of candidate metallicities should be systematically higher than the distribution of metallicities of the rest of the population.

The metallicities reported in the Kepler stellar catalog have been primarily photometrically-derived (Brown et al., 2011) and have also been shown to systematically underestimate the true metallicity and scatter (Dong et al., 2014; Petigura et al., 2017). Dong et al. (2014) used LAMOST spectroscopic data to show that the Kepler metallicities are best fit by the relation,

[TABLE]

with the relation most secure in the range $-0.3<\mathrm{[Fe/H]_{LAMOST}}<0.4$. It is thus necessary to interpret any analysis of Kepler metallicities with caution. Nevertheless, it should still be meaningful to compare the statistical distribution of candidate metallicities to the rest of the population.

In Figure 20, we display the metallicity histograms and cumulative distribution functions of all FGK stars in gray and of the candidate stars above certain probability thresholds in blue and green. It is clear to see that the metallicities of the candidates are systematically higher than the rest of the FGK stars. The candidates above the $P>0.985$ threshold are shifted further, indicating a smaller amount of false positive contamination as $P$ increases.

An identical analysis may be performed for $T_{\mathrm{eff}}$, since the occurrence of hot Jupiters drops off sharply for $T_{\mathrm{eff}}\lesssim 5000$ K. As shown in Figure 20, $T_{\mathrm{eff}}$ is systematically high for the candidates, which should be expected if the candidates are true giant planets.

7.4 KOI host stars

Our analysis has focused on the collection of Kepler FGK stars without known planets or planet candidates. Although the search for non-transiting hot Jupiters around Kepler transiting planet hosts was the original motivation of this work (Millholland et al., 2016), here we have intentionally excluded them. In KOI systems, the detection of astrometrically-induced transit timing variations for the transiting planet may serve as an additional predictor for the existence of a non-transiting giant planet. These targets are thus best handled separately, and we will be addressing them in future work.

8 Conclusion

The prospects of using phase curves to detect non-transiting, Jupiter-mass planets has been an outstanding problem for some time (Faigler & Mazeh, 2011; Placek et al., 2014; Shporer & Hu, 2015; Shporer, 2017). In this paper, we constructed a supervised learning algorithm for the detection of non-transiting Kepler planets with optical phase curves. Our algorithm relies on exploiting the time-dependent properties of planetary phase curves, namely the fact that phase curves are more temporally consistent in their amplitudes and phase compared to other types of light curve variability.

We first developed a pipeline to identify candidate phase curve signals in Kepler light curves. We demonstrated the pipeline’s recovery efficiency using light curves containing synthetically injected phase curves. The phase curves recovered by the pipeline exhibited significantly different properties from “non-phase curves”, that is, signals where the pipeline’s detection did not match the injected phase curve. We then developed a logistic regression algorithm in a supervised learning context to classify phase curves and non-phase curves. The algorithm performs exceedingly well in its ability to predict phase curves in synthetic datasets and the set of Kepler transiting hot Jupiters.

We applied our algorithm to the full set of Kepler FGK stars without known planets or planetary candidates and identified 60 high-probability planet candidates. We examined trends in the candidates’ albedos and phase offsets and discussed physical explanations for the observed trends. All of our candidates are available for inspection at https://smillholland.github.io/Non-transiting_HJs/.

For more than two decades, the characterization of hot Jupiters has been consistently improved via careful observational study from both the ground and space. Yet many zeroth-order aspects of these bizarre worlds – including their energetics, their interior structures, and their origins – remain poorly understood. Aided by the supervised machine learning detection of optical phase curves, a substantial augmentation to the known population of hot Jupiters can potentially be obtained. Planetary candidates identified in this manner can be readily verified with sparsely sampled Doppler velocity observations. These measurements will not only enable the confirmation of new planets, but will also enhance the analysis of the phase curves and aid in the understanding of the planets’ atmospheric dynamics. We are particularly optimistic that the techniques outlined here can be refined and extended to forthcoming space-based photometric surveys, including, notably, the TESS and Plato Missions.

9 Acknowledgements

We are indebted to the referee, Jon Jenkins, whose expertise in the subject and thorough review of the manuscript helped facilitate significant improvements. S.M. is supported by the National Science Foundation Graduate Research Fellowship Program under Grant Number DGE-1122492. This material is also based upon work supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under Cooperative Agreement Notice NNH13ZDA017C issued through the Science Mission Directorate. We acknowledge support from the NASA Astrobiology Institute through a cooperative agreement between NASA Ames Research Center and Yale University. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work has also made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org.