#  The population of long-period transiting exoplanets

Daniel Foreman-Mackey11affiliation<EMAIL_ADDRESS>Sagan Fellow 22affiliation:
Astronomy Department, University of Washington, Seattle, WA, 98195, USA ,
Timothy D. Morton33affiliation: Department of Astrophysics, Princeton
University, Princeton, NJ, 08544, USA , David W. Hogg44affiliation: Simons
Center for Data Analysis, 160 Fifth Avenue, 7th floor, New York, NY 10010, USA
55affiliation: Center for Cosmology and Particle Physics, New York University,
4 Washington Place, New York, NY, 10003, USA 66affiliation: Max-Planck-
Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany
77affiliation: Center for Data Science, New York University, 726 Broadway, 7th
Floor, New York, NY, 10003, USA , Eric Agol22affiliation: Astronomy
Department, University of Washington, Seattle, WA, 98195, USA , and Bernhard
Schölkopf88affiliation: Max Planck Institute for Intelligent Systems
Spemannstrasse 38, 72076 Tübingen, Germany

###### Abstract

The Kepler Mission has discovered thousands of exoplanets and revolutionized
our understanding of their population. This large, homogeneous catalog of
discoveries has enabled rigorous studies of the occurrence rate of exoplanets
and planetary systems as a function of their physical properties. However,
transit surveys like Kepler are most sensitive to planets with orbital periods
much shorter than the orbital periods of Jupiter and Saturn, the most massive
planets in our Solar System. To address this deficiency, we perform a fully
automated search for long-period exoplanets with only one or two transits in
the archival Kepler light curves. When applied to the $\sim 40,000$ brightest
Sun-like target stars, this search produces 16 long-period exoplanet
candidates. Of these candidates, 6 are novel discoveries and 5 are in systems
with inner short-period transiting planets. Since our method involves no human
intervention, we empirically characterize the detection efficiency of our
search. Based on these results, we measure the average occurrence rate of
exoplanets smaller than Jupiter with orbital periods in the range 2–25 years
to be $2.0\pm 0.7$ planets per Sun-like star.

methods: data analysis — methods: statistical — catalogs — planetary systems —
stars: statistics

††software: batman (Kreidberg, 2015), ceres (Agarwal et al., 2016), corner.py
(Foreman-Mackey, 2016), emcee (Foreman-Mackey et al., 2013), exosyspop (Morton
& Foreman-Mackey, 2016), george (Ambikasaran et al., 2016), isochrones
(Morton, 2015a), matplotlib (Hunter et al., 2007), numpy (Van Der Walt et al.,
2011), scipy (Jones et al., 2001), transit (Foreman-Mackey & Morton, 2016),
vespa (Morton, 2015b)

##  Section 1 Introduction

Data from the Kepler Mission (Borucki et al., 2011) have been used to discover
thousands of transiting exoplanets. The systematic nature of these discoveries
and careful quantification of survey selection effects, search completeness,
and catalog reliability has enabled many diverse studies of the detailed
frequency and distribution of exoplanets (for example, Howard et al., 2012;
Petigura et al., 2013; Foreman-Mackey et al., 2014; Dressing & Charbonneau,
2015; Burke et al., 2015; Mulders et al., 2015). So far, these results have
been limited to relatively short orbital periods because existing transit
search methods impose the requirement of the detection of at least three
transits within the baseline of the data. For Kepler, with a baseline of about
four years, this sets an absolute upper limit of about two years on the range
of detectable periods. In the Solar System, Jupiter – with a period of 12
years – dominates the planetary dynamics and, since it would only exhibit at
most one transit in the Kepler data, an exo-Jupiter would be missed by most
existing transit search procedures.

Before the launch of the Kepler Mission, it was predicted that the nominal
mission would discover at least 10 exoplanets with only one or two observed
transits (Yee & Gaudi, 2008), yet subsequent searches for these signals have
already been more fruitful than expected (Wang et al., 2015; Uehara et al.,
2016). However, the systematic study of the population of long-period
exoplanets found in the Kepler data to date has been hampered due to the
substantial technical challenge of implementing a search, as well as the
subtleties involved in interpreting the results. For example, false alarms in
the form of uncorrected systematics in the data and background eclipsing
binaries can make single-transit detections ambiguous.

Any single transit events discovered in the Kepler light curves are
interesting in their own right, but the development of a general and
systematic method for the discovery of planets with orbital periods longer
than the survey baseline is also crucial for the future of exoplanet research
with the transit method. All future transit surveys have shorter observational
baselines than the Kepler Mission (K2, Howell et al. 2014; TESS, Ricker et al.
2015; PLATO, Rauer et al. 2014) and given suitable techniques, single transit
events will be plentiful and easily discovered. The methodological framework
presented in the following pages is a candidate for this task.

A study of the population of long-period transiting planets complements other
planet detection and characterization techniques, such as radial velocity (for
example Cumming et al., 2008; Knutson et al., 2014; Bryan et al., 2016),
microlensing (for example Gould et al., 2010; Cassan et al., 2012; Clanton &
Gaudi, 2014; Shvartzvald et al., 2016), direct imaging (for example Bowler,
2016), and transmission spectroscopy (for example Dalba et al., 2015). The
marriage of the radial velocity and transit techniques is particularly
powerful as exoplanets with both mass and radius measurements can be used to
study planetary compositions and the formation of planetary systems (for
example Weiss & Marcy, 2014; Rogers, 2015; Wolfgang et al., 2016).
Unfortunately the existing catalog of exoplanets with measured densities is
sparsely populated at long orbital periods; this makes discoveries with the
transit method at long orbital period compelling targets for follow-up
observations. Furthermore, even at long orbital periods, the Kepler light
curves should be sensitive to planets at the detection limits of the current
state-of-the-art radial velocity surveys.

There are two main technical barriers to a systematic search for single
transit events. The first is that the transit probability for long-period
planets is very low; scaling as $\propto P^{-5/3}$ for orbital periods $P$
longer than the baseline of contiguous observations. Therefore, even if long-
period planets are intrinsically common, they will be underrepresented in a
transiting sample. The second challenge is that there are many signals in the
observed light curves caused by stochastic processes – both instrumental and
astrophysical – that can masquerade as transits. Even when the most
sophisticated methods for removing this variability are used, false signals
far outnumber the true transit events in any traditional search.

At the heart of all periodic transit search procedures is a filtering step
based on “box least squares” (BLS; Kovács et al., 2002). This step produces a
list of candidate transit times that is then vetted to remove the substantial
fraction of false signals using some combination of automated heuristics and
expert curation. In practice, the fraction of false signals can be
substantially reduced by requiring that at least three self-consistent
transits be observed (Petigura et al., 2013; Burke et al., 2014; Rowe et al.,
2015; Coughlin et al., 2016).

Relaxing the requirement of three transits requires a higher signal-to-noise
threshold per transit for validating candidate planets that display only one
to two transits. Higher signal-to-noise allows matching the candidate transit
to the expected shape of a limb-darkened light curve, as well as ruling out
various false alarms. This is analagous to microlensing surveys, for which a
planet can only be detected once, thus requiring high signal-to-noise for a
reliable detection (Gould et al., 2004).

Recent work has yielded discoveries of long-period transiting planets with
only one or two high signal-to-noise transits identified in archival Kepler
and K2 light curves by visual inspection (Wang et al., 2013; Kipping et al.,
2014b; Wang et al., 2015; Osborn et al., 2016; Kipping et al., 2016; Uehara et
al., 2016). These discoveries have already yielded some tantalizing insight
into the population of long-period transiting planets but, since these
previous results rely on human interaction, it is prohibitively expensive to
reliably measure the completeness of these catalogs. As a result, the existing
catalogs of long-period transiting planets cannot be used to rigorously
constrain the occurrence rate of long-period planets.

In this paper, we develop a systematic method for reliably discovering the
transits of large, long-period companions in photometric time series _without
human intervention_. The method is similar in character to the recently
published fully automated method used to generate the official DR24 exoplanet
candidate catalog from Kepler (Mullally et al., 2016; Coughlin et al., 2016).
Since the search methodology is fully automated, we can robustly measure the
search completeness – using injection and recovery tests – and use these
products to place probabilistic constraints on the occurrence rate of long-
period planets. We apply this method to a subset of the archival data from the
Kepler Mission, present a catalog of exoplanet candidates, and estimate the
occurrence rate of long-period exoplanets. We finish by discussing the
potential effects of false positives, evaluating the prospects for follow-up,
and comparing our results to other studies based on different planet discovery
methods.

##  Section 2 A fully automated search method

To find long-period exoplanets in the Kepler light curves, we search for
individual, high signal-to-noise transit signals using a fully automated
procedure that can be broken into three main steps:

  1. 1.

an initial candidate search using a box-shaped matched filter,

  2. 2.

light curve-level vetting (using automated model comparison) to remove signals
that don’t have a convincing transit shape, and

  3. 3.

pixel-level vetting to remove some astrophysical false positives.

The following sections describe each of these steps in more detail.

The model comparison step (step 2) is the key component of our method that
enables robust automation but it is also computationally expensive because we
must estimate the marginalized likelihoods of several different models
describing a transit and other processes that “look” like transits but are
actually caused by noise. This step is conservative: unless a signal is a very
convincing transit, it won’t pass the test. In practice, this means that all
but the highest signal-to-noise events will be rejected at this step.
Therefore, in the inexpensive first step – the initial candidate search – we
can restrict the candidate list to high signal-to-noise events without a
substantial loss in detection efficiency.

###  2.1 Step 1 – Initial candidate events

It is not computationally feasible to run a full model comparison at every
conceivable transit time in the light curve so we must first find potentially
interesting events. For our purposes, “interesting” means high signal-to-noise
and previously unknown.

To generate this list, we use a method much like the standard “box least
squares” (BLS; Kovács et al., 2002) procedure with a single (non-periodic)
box. After masking any known transits, we filter the PDC light curves (Stumpe
et al., 2012; Smith et al., 2012) using a running windowed median with a half-
width of 2 days to remove stellar variability. We then compute the signal-to-
noise of the depth of a 0.6 day-long top hat on a grid of times spanning the
full baseline of observations.

In detail, at each proposal time $t_{0}$, we hypothesize a box-shaped transit
with duration $\tau$

$\displaystyle m(t)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\mu-\delta,&\mbox{\,if
$|t-t_{0}|<\tau/2$}\\\ \mu,&\mbox{\,otherwise}\end{array}\right.\quad.$ (3)

Assuming that the uncertainties on the observed fluxes $f(t)$ are Gaussian
with known variance ${\sigma_{f}}^{2}$, the likelihood function for the mean
flux $\mu$ and transit depth $\delta$ can be analytically computed to be a
two-dimensional Gaussian with mean and covariance given by linear least-
squares. This likelihood function provides a natural scalar objective: the
signal-to-noise of the measured depth computed as a function of time. In
principle this scalar is also a function of duration but we only use a single
transit duration because the following steps in this procedure are only
sensitive to transits with very high signal-to-noise, and in practice, the
final results are insensitive to the specific choice of duration.

To avoid edge effects, we apodize this detection scalar near any large gaps in
the time series using a logistic function with width equal to one transit
duration. Finally, we estimate the background noise level in the detection
scalar time series using a robust running windowed variance estimate of the
detection scalar. We accept peaks that are more than 25-times this background
noise level as candidates.

For the Kepler light curves, this procedure yields at least one candidate
event in about 1 percent of targets. For these targets, we investigate the
three highest signal-to-noise events in the following step.

###  2.2 Step 2 – Light curve vetting

In this step of the method, the goal is to discard any signals that are not
sufficiently “transit-like” in shape. This step is similar to the method
independently developed and recently published by the Kepler team (Mullally et
al., 2016). To quantify the quality of a candidate, we perform a model
comparison between a physical transit model and a set of other parameterized
models for systematics. In order for a candidate to pass this vetting step,
the transit model must be “preferred” to any other model as measured using the
Bayesian Information Criterion (BIC). The BIC is not the optimal choice for
this model comparison, but it is more computationally tractable than the
alternatives, such as computing thousands of precise marginalized likelihoods
or expected utilities for each model. The BIC can be efficiently computed and
it exhibits the desired behavior – decreasing with increasing likelihood but
flexible models are penalized – and we find that it performs sufficiently well
in practice.

For up to three candidate transit times per light curve, we select a
contiguous chunk of PDC light curve approximately centered on the proposed
transit with no more than 500 cadences (about 10 days) and compute the BIC of
each model for this data set. The BIC for a model $k$ in the set of $K$ models
is given by

$\displaystyle\mathrm{BIC}_{k}$ $\displaystyle=$
$\displaystyle-2\,\ln\mathcal{L}^{*}+J\,\ln N$ (4)

where the likelihood function $\mathcal{L}$ is evaluated at its maximum, $J$
is the number of free parameters in the model, and $N$ is the number of data
points in the data set.

For each model, we describe the data using a Gaussian Process (GP; Rasmussen &
Williams, 2006) with a Matérn-3/2 covariance and mean given by the chosen
model $m_{k}(t;\,{{\boldsymbol{\theta}}})$ parameterized by the parameter
vector ${\boldsymbol{\theta}}$.

We consider the following mean models (this list provides a qualitative
justification for each model):

  * •

transit – a limb-darkened, exposure-time integrated transit light curve,

  * •

variability – a pure GP model to capture stellar variability,

  * •

outlier – a single outlier to account for a bad data point,

  * •

step – a step function to describe sudden pixel sensitivity dropouts (SPSDs;
for example Christiansen et al., 2013), and

  * •

box – a box to catch signals that are well fit by the search scalar but
insufficiently transit-like to be convincing.

The functional forms of these models are given in Appendix A and the details
of the technical methodology of GP fitting are described in Appendix B.

Figure 1 shows representative events that fall into different classes and the
corresponding maximum likelihood model. For each candidate event, the BIC of
each of these models is computed and the event is only passed as a candidate
if the transit model is preferred to all the other models. The box model is
the most restrictive comparison, vetoing about half of the candidate events in
the Kepler light curves, followed by the variability model. To further
restrict to non-grazing transits, we also reject events where the maximum
likelihood impact parameter is greater than $1-R_{\mathrm{P}}/R_{\star}$.

Since the search procedure described here was tuned to discover transit
signals, we do not consider the distribution or potential astrophysical nature
of any models besides the transit model. In the future, it would be
interesting to relax this goal and investigate the other model classes; in
particular, the box model is sensitive to astrophysical phenomena, notably
occultations of white dwarfs. In a cursory investigation it is clear that the
majority of signals labeled box in our analysis are noise; however, a subset
are likely to be astrophysical in nature.

The reliability of this method of automated vetting is limited by the specific
models selected in this step. We find that these are sufficient for the
targets discussed below but different target lists or data sets might require
additional models to be included for robust selection.

Figure 1: Representative examples of candidate events flagged by the initial
search. Each example falls into a different model category and the figure
shows the data as black points and the best fit mean model prediction. The
examples represent the following model categories: _(a)_ step _(b)_
variability, _(c)_ box, and _(d)_ transit.

###  2.3 Step 3 – Pixel-level vetting

To minimize contamination from background eclipsing binary systems, we require
candidate events to pass a centroid shift test similar to the one used in the
official Kepler transit search pipeline (Bryson et al., 2013). To measure the
centroid shift, we model the flux-weighted centroid traces independently in
each coordinate as a multiple of the best-fit transit model and a GP noise
model. By properly normalizing the transit model, we measure the in-transit
centroid shift $\Delta_{\mathrm{centroid}}$ in pixels. We reject any candidate
event where the estimated transit location is more than half a pixel from the
out-of-transit centroid

$\displaystyle\Delta_{\mathrm{centroid}}\,\left(\frac{1}{\delta}-1\right)$
$\displaystyle>$ $\displaystyle 0.5$ (5)

where $\delta$ is the observed transit depth (Bryson et al., 2013).

##  Section 3 Results: a catalog of long-period transiting exoplanet
candidates

To limit the scope of this paper while still demonstrating the applicability
of our method, we search the Kepler archival light curves of the brightest and
quietest Sun-like stars for long-period transiting exoplanets. In this
section, we describe the target selection process and the parameter estimation
procedure.

###  3.1 Target selection

We select the $\sim 40,000$ brightest and quietest G and K dwarfs from the
Kepler catalog using the most recent catalog of stellar
parameters111Parameters from the q1_q17_dr24_stellar table from the NASA
Exoplanet Archive (Huber et al., 2014, with updates). and the cuts used by
Burke et al. (2015):

  * •

$4200{\,\mathrm{K}}\leq T_{\mathrm{eff}}\leq 6100{\,\mathrm{K}}$,

  * •

$R_{\star}\leq 1.15\,R_{\odot}$,

  * •

$\mathrm{data\,span}\geq 2\,\mathrm{years}$,

  * •

$\mathrm{duty\,cycle}\geq 0.6$,

  * •

$K_{p}\leq 15{\,\mathrm{mag}}$, and

  * •

$\mathrm{CDPP}_{7.5{\,\mathrm{hrs}}}\leq 1000{\,\mathrm{ppm}}$.

We continue by excluding the light curves of known eclipsing
binaries222http://keplerebs.villanova.edu/ (Kirk et al., 2016), other known
false positives (Coughlin et al., 2016), a planet with known transit timing
variations (Kepler-9), and four especially noisy stars (KIC 4482348, KIC
4450472, KIC 5438845, and KIC 10068041). The final catalog contains 39,036
targets and the parameter distribution is shown in Figure 2.

Since these data have already been searched for short-period planets, we
assume that all high signal-to-noise candidates with three or more transits
have been previously found (Coughlin et al., 2016). To remove these candidates
from consideration, we mask the cadences within two transit durations of the
time when a short-period planet candidate is known to transit333We
specifically use the q1_q17_dr24_koi from the NASA Exoplanet Archive
http://exoplanetarchive.ipac.caltech.edu/..

Figure 2: The distribution of stellar parameters for Kepler targets selected
for this search (orange) compared to the distribution of the full Kepler
target catalog (black).

###  3.2 Parameter estimation

For each transit candidate, we constrain the physical parameters of the system
by fitting a section of light curve around each transit using an exposure-time
integrated Keplerian orbit with a quadratic limb darkening law for the central
body Foreman-Mackey & Morton (2016). It has previously been established that
the orbital period of a transiting planet with only one transit can still be
constrained given a measurement of the stellar density and an assumption about
the orbital eccentricity (for example Wang et al., 2015; Osborn et al., 2016).
Qualitatively this works because the transit of a bound body cannot have an
arbitrary period for a given duration. This is the same argument used to
justify the “photoeccentric effect” (Dawson & Johnson, 2012) and the method of
“asterodensity profiling” (Kipping et al., 2014a). In particular, this
suggests that the periods of single transits in systems with multiple inner
planets will be especially well constrained (Kipping et al., 2012). In this
paper, we do not take advantage of the extra constraints provided by the inner
planets, instead treating each long-period transiting system in isolation, but
this would be a good follow-up project.

In the following paragraphs, we describe the components of the probabilistic
model used to infer the planet candidates’ properties. To perform parameter
estimation under this model, we use the Markov Chain Monte Carlo (MCMC)
package emcee444http://dfm.io/emcee (Foreman-Mackey et al., 2013) with an
ensemble of 40 walkers. We run each chain until at least 750 independent
samples – in most cases, we actually produce thousands of independent samples
– are obtained555The integrated autocorrelation time is estimated using a
robust iterative method as suggested by Alan Sokal:
http://www.stat.unc.edu/faculty/cji/Sokal.pdf. and discard the first third of
the chain as burn-in. The posterior constraints on a few physical parameters
for the single transit candidate in the light curve of KIC 8505215 are shown
in Figure 3 and all the chains are made available
online666http://dx.doi.org/10.5281/zenodo.58273.

Figure 3: The posterior constraints on the physical parameters for the single
transit candidate found in the light curve of KIC 8505215. The contour plots
show estimates of the two-dimensional marginalized probability densities and
the histograms show the marginalized density for each parameter. This figure
was generated using corner.py (Foreman-Mackey, 2016).

#### Priors

For each candidate in our sample, we take the constraints on the stellar
parameters from the Kepler DR24 stellar properties catalog and assume an
empirical beta function prior on the eccentricities based on the observed
eccentricity distribution of long-period planets discovered using radial
velocities (Kipping, 2013a). Table 1 lists all the fit parameters and their
prior distributions. Besides these listed priors, we add the extra constraint
that no other transits can occur in the baseline of the Kepler observations.
This constraint is overly conservative because there is some probability that
a second transit could occur in a data gap but we find that, in practice, most
of the posterior mass is at longer periods and the period inferences are not
significantly affected.

#### Likelihood function

As above, we model the light curve as a Gaussian Process (GP) with a physical
transit model as the mean, and a covariance matrix described by a Matérn-3/2
kernel function. The full likelihood function and some details of GP
regression are given in Appendix B. For computational efficiency, we first
perform a joint optimization of the physical parameters and GP hyperparameters
to find the maximum _a posteriori_ model then keep the hyperparameters fixed
and run MCMC sampling for the 11 physical parameters alone.

Table 1: The inferred parameters and priors used in the inference

name | symbol | units | prior
---|---|---|---
mean flux | $\log f_{\star}$ |  $\cdots$ | $\log f_{\star}\sim\mathcal{U}(-1,\,1)$
stellar massaaStellar parameters and uncertainties taken from the Kepler catalog (Huber et al., 2014) | $M_{\star}$ | $M_{\odot}$ | $M_{\star}\sim\mathcal{N}(M_{\star,\mathrm{cat}},\,\sigma_{M,\star,\mathrm{cat}})$
stellar radiusaaStellar parameters and uncertainties taken from the Kepler catalog (Huber et al., 2014) | $R_{\star}$ | $R_{\odot}$ | $R_{\star}\sim\mathcal{N}(R_{\star,\mathrm{cat}},\,\sigma_{R,\star,\mathrm{cat}})$
limb darkening | $q_{1}$ |  $\cdots$ | $q_{1}\sim\mathcal{U}(0,\,1)$
$q_{2}$ |  $\cdots$ | $q_{2}\sim\mathcal{U}(0,\,1)$
planet radius | $\log R_{\mathrm{P}}$ | $R_{\odot}$ | $\log R_{\mathrm{P}}\sim\mathcal{U}(-10,\,2)$
reference time | $t_{0}$ | days |  $t_{0}\sim\mathcal{U}(t_{\mathrm{cand}}-0.5,\,t_{\mathrm{cand}}+0.5)$bbThe reference time is constrained to be within half a day of the candidate transit time
semi-major axis & inclination | $\sqrt{a}\sin i$ | ${R_{\odot}}^{1/2}$ | $\sqrt{a}\sin i\sim\mathcal{U}(-10^{3},\,10^{3})/\sqrt{a}$
$\sqrt{a}\cos i$ | ${R_{\odot}}^{1/2}$ | $\sqrt{a}\cos i\sim\mathcal{U}(0,\,10^{3})/\sqrt{a}$
eccentricity | $\sqrt{e}\sin\omega$ |  $\cdots$ |  $e\sim\beta(1.12,\,3.09)$ccKipping (2013b)
$\sqrt{e}\cos\omega$ |  $\cdots$ | $\omega\sim\mathcal{U}(-\pi,\,\pi)$

Note. — There is one further constraint that complicates these priors: the
period of the orbit must be longer than some minimum period $P_{\mathrm{min}}$
set by the transit time and the full baseline of Kepler observations.

##  Section 4 Catalog of transit candidates

Applying the search procedure described in Section 2 to the Kepler light
curves of the 39,036 targets selected in Section 3.1, we find 16 convincing
transit candidates. Visual inspection of each candidate confirms the
reliability of the classification and no candidates are manually removed from
the catalog. Of these, three candidates have two transits in the Kepler
baseline and the remainder have only one observable transit. The candidates
and their inferred physical parameters are listed in Table 2 and the light
curves are plotted in Figure 4. The inferred radius and orbital periods of the
candidates are compared to the short-period Kepler sample and the Solar System
in Figure 5.

Two of the shortest period candidates – both with two observed transits – have
previously been studied in detail (KIC 8800954 and KIC 3239945; Kipping et
al., 2014b, 2016). Table 2 indicates the candidates that were also discovered
by earlier searches for long-period transiting systems using visual inspection
(Wang et al., 2015; Uehara et al., 2016). The consistency between our results
and the earlier catalogs is reassuring. In the light curves of targets with
previously known short-period planets, our automated search did not find any
candidates that weren’t previously detected by visual inspection (Uehara et
al., 2016) and one candidate (KIC 3230491) reported by the human analysis was
discarded as grazing by our search. The Planet Hunters citizen science project
(Fischer et al., 2012) reported five long-period candidates with one or two
observed transits in our target list (Wang et al., 2015). Of these, we also
find two (KIC 8410697 and KIC 10842718) although we find a second transit in
the KIC 8410697 system that was previously missed. We do not recover the three
other candidates reported by Wang et al. (2015): KIC 5536555, KIC 9662267, and
KIC 12454613. The transits of these candidates are all low signal-to-noise and
they do not pass our initial signal-to-noise threshold. Six of the candidates
in Table 2 have not been previously published.

Of the 16 candidates, 5 have known inner planets with three or more observable
transits (Coughlin et al., 2016). Given the fact that only 844 of the 39,036
targets had previously known planets, this means that systems with short-
period transiting planets are nearly a factor of 20 more likely to host long-
period planets accessible by our method than systems with no known inner
transiting planets. This difference cannot be accounted for by differences in
completeness between targets with known planets and without because the
average detection efficiency for both populations is consistent within
sampling uncertainty. Qualitatively, this suggests that these long-period
planets occur with a higher frequency in multi-planet systems or are
preferentially aligned with the plane of any inner planets but a more detailed
analysis would be needed to make a quantitative statement (see, for example,
Tremaine & Dong, 2012; Fang & Margot, 2012; Ballard & Johnson, 2016; Moriarty
& Ballard, 2015).

The candidate in the light curve of KIC 4754460 is an individual transit
candidate but another deeper eclipse can be found at a Kepler Barycentric
Julian Date (KBJD) of 1587.13; right at the beginning of Quarter 17. This
eclipse was missed by the automated search because only the second half of the
eclipse is observed. The most likely explanation of this system is that the
listed candidate is the secondary eclipse of a binary system but we will keep
the candidate in the list and treat this effect statistically in Section 7.

Five candidate transit events in the light curves of four targets were
rejected because of a significant centroid shift or a large impact parameter.
These events are probably astrophysical eclipses from binary star systems that
were not found by previous studies of long-period eclipsing binary systems. We
do not consider these events further in the following analysis but Table 3
lists these events and their properties for posterity.

Table 2: The inferred parameters for the long-period transiting exoplanet
candidates

kic id | $T_{\mathrm{eff}}$ | $R_{\star}$ | Kp | period | $t_{0}$ | radius | duration | impact |  $T_{\mathrm{eq}}$**The equilibrium temperature is computed assuming zero albedo. | $\mathrm{Pr}_{\mathrm{planet}}$ | KOI/Kepler$\dagger$$\dagger$footnotemark:
---|---|---|---|---|---|---|---|---|---|---|---
| K | $R_{\odot}$ |  | years | KBJD | $R_{\mathrm{J}}$ | hours |  | K |  |
3218908b | $5513_{-139}^{+172}$ | $0.75_{-0.05}^{+0.26}$ | $14.6$ | $7.0_{-3.4}^{+9.5}$ | $766.6722_{-0.0114}^{+0.0096}$ | $0.514_{-0.093}^{+0.092}$ | $21.45_{-0.61}^{+0.72}$ | $0.24_{-0.17}^{+0.21}$ | $129_{-39}^{+41}$ | $0.96$ | 1108 / 770
3239945c | $4786_{-88}^{+106}$ | $0.71_{-0.03}^{+0.05}$ | $14.0$ | $2.9328721_{-0.0000026}^{+0.0000026}$ | $420.28714_{-0.00068}^{+0.00069}$ | $0.876_{-0.039}^{+0.039}$ | $16.202_{-0.071}^{+0.077}$ | $0.207_{-0.110}^{+0.071}$ | $142.8_{-4.3}^{+4.3}$ | $0.73$ | 490 / 167
4754460 | $5766_{-158}^{+172}$ | $1.13_{-0.22}^{+0.38}$ | $14.9$ | $5.9_{-3.0}^{+11.8}$ | $826.8369_{-0.0046}^{+0.0046}$ | $0.67_{-0.15}^{+0.16}$ | $15.92_{-0.54}^{+0.55}$ | $0.893_{-0.037}^{+0.018}$ | $171_{-64}^{+60}$ | $0.95$ |
6551440 | $6050_{-182}^{+155}$ | $1.10_{-0.15}^{+0.43}$ | $13.6$ | $4.0_{-1.2}^{+4.2}$ | $1039.0589_{-0.0037}^{+0.0037}$ | $0.282_{-0.083}^{+0.093}$ | $10.85_{-0.30}^{+0.37}$ | $0.60_{-0.37}^{+0.20}$ | $170_{-45}^{+38}$ | $0.97$ |
8410697c,a | $5918_{-152}^{+157}$ | $1.00_{-0.16}^{+0.35}$ | $13.4$ | $2.8688097_{-0.0000054}^{+0.0000053}$ | $542.1231_{-0.0013}^{+0.0013}$ | $0.698_{-0.078}^{+0.107}$ | $19.77_{-0.10}^{+0.12}$ | $0.15_{-0.11}^{+0.12}$ | $206_{-13}^{+15}$ | $0.95$ |
8426957 | $5992_{-183}^{+153}$ | $1.08_{-0.13}^{+0.46}$ | $13.6$ | $54_{-36}^{+88}$ | $784.677_{-0.013}^{+0.013}$ | $1.04_{-0.25}^{+0.30}$ | $39.4_{-1.4}^{+1.6}$ | $0.889_{-0.059}^{+0.037}$ | $85_{-28}^{+46}$ | $0.80$ |
8505215b | $5087_{-98}^{+102}$ | $0.71_{-0.03}^{+0.06}$ | $13.0$ | $9.1_{-3.4}^{+9.5}$ | $140.0492_{-0.0018}^{+0.0017}$ | $0.277_{-0.017}^{+0.017}$ | $20.06_{-0.16}^{+0.18}$ | $0.28_{-0.19}^{+0.20}$ | $103_{-23}^{+19}$ | $0.96$ | 99 / none
8738735b | $6000_{-129}^{+101}$ | $1.04_{-0.06}^{+0.21}$ | $13.9$ | $9.9_{-5.0}^{+14.9}$ | $697.8538_{-0.0049}^{+0.0059}$ | $0.355_{-0.044}^{+0.045}$ | $27.44_{-0.37}^{+0.62}$ | $0.28_{-0.20}^{+0.30}$ | $137_{-39}^{+43}$ | $0.97$ | 693 / 214
8800954c | $5286_{-101}^{+110}$ | $0.76_{-0.04}^{+0.06}$ | $13.4$ | $1.9279957_{-0.0000091}^{+0.0000092}$ | $492.7652_{-0.0024}^{+0.0024}$ | $0.386_{-0.025}^{+0.025}$ | $15.76_{-0.13}^{+0.14}$ | $0.18_{-0.13}^{+0.17}$ | $189.4_{-7.4}^{+7.2}$ | $0.95$ | 1274 / 421
9306307 | $5762_{-160}^{+209}$ | $0.92_{-0.15}^{+0.44}$ | $14.0$ | $4.3_{-1.1}^{+3.3}$ | $1191.35648_{-0.00018}^{+0.00018}$ | $1.22_{-0.36}^{+0.49}$ | $8.499_{-0.042}^{+0.042}$ | $0.6399_{-0.0072}^{+0.0068}$ | $126_{-32}^{+33}$ | $8.7\times 10^{-6}$ |
10187159b | $5185_{-143}^{+179}$ | $0.91_{-0.11}^{+0.80}$ | $14.4$ | $4.9_{-1.8}^{+7.6}$ | $604.1102_{-0.0031}^{+0.0023}$ | $0.43_{-0.13}^{+0.21}$ | $11.81_{-0.20}^{+0.27}$ | $0.23_{-0.16}^{+0.21}$ | $119_{-43}^{+34}$ | $0.91$ | 1870 / 989
10287723b | $4500_{-119}^{+153}$ | $0.73_{-0.05}^{+0.04}$ | $13.4$ | $4.9_{-1.3}^{+4.3}$ | $393.5976_{-0.0029}^{+0.0031}$ | $0.266_{-0.024}^{+0.027}$ | $9.49_{-0.30}^{+0.33}$ | $0.68_{-0.36}^{+0.13}$ | $114_{-22}^{+13}$ | $0.95$ | 1174 / none
10321319 | $5749_{-133}^{+154}$ | $0.94_{-0.13}^{+0.38}$ | $11.9$ | $5.5_{-2.1}^{+8.1}$ | $554.3562_{-0.0063}^{+0.0064}$ | $0.163_{-0.037}^{+0.046}$ | $16.84_{-0.38}^{+0.42}$ | $0.47_{-0.31}^{+0.26}$ | $153_{-47}^{+37}$ | $1.00$ |
10602068 | $5628_{-157}^{+165}$ | $0.91_{-0.08}^{+0.36}$ | $14.9$ | $3.16_{-0.83}^{+2.65}$ | $830.80892_{-0.00015}^{+0.00015}$ | $2.00_{-0.35}^{+0.66}$ | $12.804_{-0.092}^{+0.068}$ | $0.6027_{-0.0078}^{+0.0072}$ | $159_{-36}^{+28}$ | $3.9\times 10^{-9}$ |
10842718a | $5754_{-156}^{+159}$ | $1.04_{-0.14}^{+0.38}$ | $14.6$ | $12.7_{-6.6}^{+20.2}$ | $226.2344_{-0.0047}^{+0.0047}$ | $0.74_{-0.16}^{+0.16}$ | $35.92_{-0.38}^{+0.51}$ | $0.26_{-0.17}^{+0.19}$ | $128_{-43}^{+47}$ | $0.90$ |
11709124b | $5688_{-101}^{+108}$ | $0.97_{-0.09}^{+0.19}$ | $14.5$ | $4.3_{-1.3}^{+4.7}$ | $657.2674_{-0.0016}^{+0.0018}$ | $0.83_{-0.11}^{+0.12}$ | $17.75_{-0.44}^{+0.54}$ | $0.51_{-0.22}^{+0.12}$ | $166_{-39}^{+28}$ | $0.94$ | 435 / 154

††footnotetext: The KOI number and, if applicable, the Kepler number for the
target.

ccfootnotetext: Candidate has two observed transits.

Note. — The values and uncertainties indicate the 16-th, 50-th, and 84-th
percentiles of the posterior samples for each parameter.

Figure 4: Sections of PDC light curve centered on each candidate (black) with
the posterior-median transit model over-plotted (orange). The y-axis shows the
relative apparent flux of the light curve in parts per thousand (ppt).
Candidates with two transits are folded on the posterior-median period. The
plots are ordered by increasing planetary radius from the top-left to the
bottom-right.
Figure 5: The catalog of long-period transiting exoplanet candidates (green
points with error bars) compared to the Kepler candidates (blue points) and
confirmed planets (black points; Morton et al., 2016) found in our target
sample, and the Solar System (orange squares). The thin black error bars to
the left of each candidate indicate the minimum period allowed for each
candidate by the prior assumption that no other transit occurred during the
baseline of Kepler observations of the target. The vertical solid line shows
the absolute maximum period accessible to transit searches that require at
least three transits in the Kepler data.

Table 3: The signals rejected with a centroid shift or large impact parameter  kic id | time | depth | duration | reason
---|---|---|---|---
| KBJD | ppt | hours |
3230491 | 315.3 | 9.0 | 7.4 | impact
6342758 | 553.9 | 10.3 | 9.9 | impact
8463272 | 641.0 | 35.5 | 4.8 | impact
8463272 | 1206.7 | 35.5 | 4.8 | impact
10668646 | 1449.2 | 5.7 | 12.4 | centroid

##  Section 5 Empirical search completeness

To measure the completeness of the search procedure described in Section 2, we
exploit the fact that transit signals are sparse and rare. Therefore, most
light curves contain no transits and we can reliably measure the recovery rate
of our method on synthetic transit signals – with known properties – injected
into real light curves. This procedure is standard practice in the transit
literature and it has been used to determine the completeness of the KOI
catalog (Christiansen et al., 2013, 2015) and other independent transit
searches (Petigura et al., 2013; Dressing & Charbonneau, 2015; Foreman-Mackey
et al., 2015).

To reliably capture the full structure of the search completeness function,
the simulations must sample the (high-dimensional) space of all properties
that affect the probability of detecting a transit: the stellar properties
(including variability amplitudes and time scales), the planet’s physical
properties and orbital elements, and any observational effects (noise,
spacecraft pointing variations, _etc._). For the modest goals of this paper,
we only need a robust constraint on the transit detection efficiency
_integrated_ across the target sample but, even so, many simulations per star
are required.

The procedure for measuring the recovery rate of simulated transits is as
follows:

  1. 1.

First, a star is randomly selected from the target list, and the PDC light
curve and stellar properties for that star are loaded.

  2. 2.

Planetary properties are sampled from the distributions listed in Table 4 with
phase uniformly distributed across the baseline of observations. These
properties are re-sampled until the transit is visible in at least one non-
flagged cadence.

  3. 3.

The transit signal induced by this planet is computed and multiplied into the
PDC light curve.

  4. 4.

The transit search method described in Section 2 – including de-trending and
_all automated vetting_ – is applied to this light curve with the injected
transit signal.

  5. 5.

This candidate is flagged as recovered if at least one transit within one
transit duration passes all the cuts imposed by the automated vetting.

The fraction of recovered simulations as a function of the relevant parameters
gives an estimate of the probability of detecting an exoplanet transit with a
given set of parameters, _conditioned on the fact that it transits the star
during a time when the star was being observed by Kepler_. We will call this
function $Q_{\mathrm{det},k}({{\boldsymbol{w}}})$ where ${\boldsymbol{w}}$ is
the set of all parameters affecting the transit detectability and $k$ is an
index running over target stars.

Figure 6 shows the fraction of recovered simulations as a function of planet
radius and orbital period based on 819,752 injected signals. This figure shows
the transit detection efficiency falling with decreasing planet radius. This
is the expected behavior because the depth (and signal strength) of a transit
scales with the area ratio between the planet and the star. There is also a
slight decrease in the completeness to larger planet radius. This trend is
introduced in steps 2 and 3 of the search procedure where the tuning
parameters were chosen to maximize the yield of convincing small transit
discoveries. The decreasing completeness with orbital period is less intuitive
because, on average, the signal strength should increase as the duration of
the transit increases. In this case, this simplistic treatment misses two
important factors. First, in step 1 of the search procedure (Section 2.1) only
a single transit duration is used and second, longer transits are less easily
distinguished from stellar variability and they will, therefore, be discarded
in the conservative light curve vetting step (Section 2.2).

This detection efficiency must then be combined with the geometric transit
probability function and the window function. For the star $k$, the geometric
transit probability is given by (Winn, 2010)

$\displaystyle Q_{\mathrm{geom},k}({{\boldsymbol{w}}})$ $\displaystyle=$
$\displaystyle\frac{R_{\star,k}+R}{a_{k}}\,\frac{1+e\,\sin\omega}{1-e^{2}}$
(6) $\displaystyle=$
$\displaystyle\left[\frac{4\,\pi^{2}}{G\,M_{\star,k}}\right]^{1/3}\,\left[\frac{1+e\,\sin\omega}{1-e^{2}}\right]\,(R_{\star,k}+R)\,P^{-2/3}$
(7)

where $R$ is the planet radius, $P$ is the orbital period, $e$ is the orbital
eccentricity, $\omega$ is the argument of periastron, $R_{\star,k}$ is the
radius of star $k$, and $M_{\star,k}$ is the star’s mass. All of these
parameters are included in ${\boldsymbol{w}}$.

In Equation (6), the term $(R_{\star,k}+R)$ takes grazing transits into
account. This might seem counter intuitive because, as part of the search
procedure, we rejected candidates where the maximum likelihood model had a
grazing transit. However, since the measurement of $Q_{\mathrm{det},k}$
included a cut on the measured impact parameter, the $Q_{\mathrm{det},k}$ term
already takes this effect into account. In other words, $Q_{\mathrm{det},k}$
quantifies the probability that a transit of a given shape will be detected
_given that it transits at all_ and $Q_{\mathrm{geom},k}$ – the way it is
written in Equation (6) – is the marginalized probability that the system will
transit given its physical parameters.

Approximating the window function using a binomial probability of observing at
least one transit, we find (following Burke & McCullough, 2014)

$\displaystyle Q_{\mathrm{win},k}({{\boldsymbol{w}}})$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}1-(1-f_{\mathrm{duty},k})^{T_{k}/P}&\quad\mbox{if
$P\leq T_{k}$}\\\
T_{k}\,f_{\mathrm{duty},k}/P&\quad\mbox{otherwise}\end{array}\right.$ (10)

where $f_{\mathrm{duty},k}$ is the duty cycle and $T_{k}$ is the full
observation baseline for target $k$.

Combining these detection effects, the total detection efficiency is given by

$\displaystyle Q_{k}({{\boldsymbol{w}}})$ $\displaystyle=$ $\displaystyle
Q_{\mathrm{det},k}({{\boldsymbol{w}}})\,Q_{\mathrm{win},k}({{\boldsymbol{w}}})\,Q_{\mathrm{geom},k}({{\boldsymbol{w}}})\quad.$
(11)

So that our planet candidate catalog can be easily used for other projects, we
also provide an analytic approximation to the relevant integrated detection
efficiency function

$\displaystyle Q_{\mathrm{det}}(P,\,R)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{K}\int
Q_{\mathrm{det,k}}({{\boldsymbol{w}}})\,p({{\boldsymbol{w}}}_{\\{P,\,R\\}})\,\mathrm{d}{{\boldsymbol{w}}}_{\\{P,\,R\\}}$
(12)

where $p({{\boldsymbol{w}}}_{\\{P,\,R\\}})$ is the prior distribution of all
the parameters except the period and radius. We find that a good fit to this
integrated completeness is given by the function

$\displaystyle Q_{\mathrm{det}}(P,\,R)$ $\displaystyle\approx$
$\displaystyle\frac{\mathrm{min}[\mathrm{max}[a(P)\,b(R),\,0],\,1]}{1+\exp\left[-k(P)\,(\ln
R/R_{\mathrm{J}}-x(P))\right]}$ (13)

where

$\displaystyle a(P)=a_{1}\,\ln P/{\,\mathrm{yr}}+a_{2}\,,\quad b(R)=b_{1}\,\ln
R/R_{\mathrm{J}}+b_{2}\,,\quad$ (14) $\displaystyle k(P)=k_{1}\,\ln
P/{\,\mathrm{yr}}+k_{2}\,,\,\mathrm{and}\quad x(P)=x_{1}\,\ln
P/{\,\mathrm{yr}}+x_{2}\,.$ (15)

When fit to the set of 819,752 injected transits, the best fit parameters are
given in Table 5 and the approximation is plotted in Figure 7. Note that we do
not use this approximation in the following analysis but instead compute the
relevant integrals using the injection results directly.

Table 4: Distributions of physical parameters for transit simulations  name | distribution
---|---
period | $\log P\sim\mathcal{U}(\log 2{\,\mathrm{yr}},\,\log 25{\,\mathrm{yr}})$
radius ratio | $\log R_{\mathrm{P}}/R_{\star}\sim\mathcal{U}(\log 0.02,\,\log 0.2)$
impact parameter | $b\sim\mathcal{U}(0,\,1+R_{\mathrm{P}}/R_{\star})$
eccentricity |  $e\sim\beta(1.12,\,3.09)$aaIncluded in the Wang et al. (2015) catalog.
$\omega\sim\mathcal{U}(-\pi,\,\pi)$
limb darkening | $q_{1}\sim\mathcal{U}(0,\,1)$
$q_{2}\sim\mathcal{U}(0,\,1)$

Figure 6: An empirical estimate of the search completeness as a function of
planet radius and orbital period. In each bin, the completeness is estimated
by the fraction of recovered simulations. The projected histograms show the
integrated completeness as independent functions of period and radius.
Figure 7: An analytic approximation to Figure 6 with the same color scale. The
contours indicate the 0.1, 0.3, 0.5, and 0.7 levels.

Table 5: The fit parameters for the analytic approximation to the completeness function  parameter | value | parameter | value
---|---|---|---
$a_{1}$ | $-0.13$ | $k_{1}$ | $0.70$
$a_{2}$ | $0.95$ | $k_{2}$ | $3.06$
$b_{1}$ | $-0.20$ | $x_{1}$ | $-0.07$
$b_{2}$ | $0.90$ | $x_{2}$ | $-0.91$

##  Section 6 The occurrence rate of long-period exoplanets

Using the catalog of exoplanet discoveries (Section 4) and the measurement of
the search completeness (Section 5), we can now estimate the occurrence rate
of long-period exoplanets. To simplify the analysis, we will make the strong
assumption that none of the candidates are astrophysical false positives (the
eclipse or occultation of a stellar mass companion, around either the target
star or a faint background star). We revisit this assumption and discuss its
validity in the following section. As a further simplification, we also
neglect the measurement uncertainties on the planet parameters (including
orbital period). This assumption is justified because we are only making high-
level measurements of the mean occurrence rate in bins larger than the
uncertainties.

Assuming a Poisson likelihood, the occurrence rate density in a volume $V$ –
defined as $P_{\mathrm{min}}\leq P<P_{\mathrm{max}}$ and $R_{\mathrm{min}}\leq
R<R_{\mathrm{max}}$ – is (see, for example, the Appendix of Foreman-Mackey et
al., 2014)

$\displaystyle\Gamma_{V}\equiv\frac{\,\mathrm{d}^{2}N}{\,\mathrm{d}\ln
P\,\mathrm{d}\ln R}$ $\displaystyle=$
$\displaystyle\frac{C(P_{\mathrm{min}},\,P_{\mathrm{max}};\,R_{\mathrm{min}},\,R_{\mathrm{max}})}{Z(P_{\mathrm{min}},\,P_{\mathrm{max}};\,R_{\mathrm{min}},\,R_{\mathrm{max}})}$
(16)

where $N$ is the expected number of planets per G/K dwarf, $C(\cdots)$ is the
number of detected planets in the volume, and

$\displaystyle
Z(P_{\mathrm{min}},\,P_{\mathrm{max}};\,R_{\mathrm{min}},\,R_{\mathrm{max}})$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{K}\int
p({{\boldsymbol{w}}}_{\\{P,\,R\\}})\,Q_{k}({{\boldsymbol{w}}})\,{\boldsymbol{1}}[P,\,R\in
V]\,\mathrm{d}{{\boldsymbol{w}}}$ (17)

where $p({{\boldsymbol{w}}}_{\\{P,\,R\\}})$ is the prior distribution of all
the parameters except the period and radius and ${\boldsymbol{1}}[\cdot]$ is 1
if the argument is satisfied and 0 otherwise. Using the $J$ injections sampled
uniformly in period and radius and other parameters from $p(w_{\\{P,\,R\\}})$,

$\displaystyle
Z(P_{\mathrm{min}},\,P_{\mathrm{max}};\,R_{\mathrm{min}},\,R_{\mathrm{max}})$
$\displaystyle\approx$
$\displaystyle\frac{K\,V}{J}\,\sum_{j=1}^{J}Q_{k_{j}}(w^{(j)})$ (18)

where the sum is over all injections in the volume $V$.

Using the injection results from Section 5 and the catalog of discoveries from
Section 4, we compute the occurrence rate in the period range 2 to 25 years
and in two radius bins between 0.1 and $1.0\,R_{\mathrm{J}}$. The calculated
occurrence rates are listed in Table 6. Integrating the two bin model in this
range, we find an expected occurrence rate of

$\displaystyle
N_{0.1\,R_{\mathrm{J}}-1\,R_{\mathrm{J}},\,2{\,\mathrm{yr}}-25{\,\mathrm{yr}}}$
$\displaystyle=$ $\displaystyle 2.00\pm 0.72{\,\mathrm{planets}}$ (19)

per G/K dwarf with radii in the range $0.1\,R_{\mathrm{J}}-1\,R_{\mathrm{J}}$
and periods in the range $2{\,\mathrm{yr}}-25{\,\mathrm{yr}}$. This result is
qualitatively consistent with the Solar System where there is one planet –
Jupiter – in this parameter range and Saturn is just outside the range with an
orbital period of 29 years. In Section 8, we compare with similar occurrence
rate estimates from the literature.

The occurrence rates given here should be interpreted with a few caveats in
mind. First, when we inferred the periods of the planets with only one
transit, we assumed that the period was long enough that no other transit
occurred during the Kepler lifetime. This neglects the small but non-
negligible posterior probability – less than one percent for the typical
candidate – that a second transit might have occurred in a data gap. All of
the candidates in our catalog are consistent with having periods this long but
the geometric transit probability decreases quickly with orbital period. For
the purposes of this paper, we neglect this effect because its rigorous
treatment is subtle, but comment that this would only ever decrease the
occurrence rate estimate. Second, we assume that each planet candidate
transits the star that is characterized by Huber et al. (2014); we assert that
each planet does not transit a fainter companion star or a background star. If
the planet does transit a companion star, then the companion star must be
fainter, and hence denser, causing the period to be underestimated. If the
planet transits a background star, it is more likely to be a giant star due to
Malmquist bias, hence the density of the star and period of the planet would
be overestimated. Either of these scenarios has a small probability, so we
expect that our population estimates will stand, while the parameter estimates
for individual candidates should be taken as provisional until more detailed
follow-up is carried out, including high-contrast imaging, high-resolution
spectroscopy, and parallax measurements. Third, we assume that the Huber et
al. (2014) parameters are accurate for each star that is transited by a planet
candidate, and that each transit is unaffected by blending. Malmquist bias,
Eddington bias, and metallicity bias may affect the stellar parameters (Gaidos
& Mann, 2012), and so we again caution that the individual parameter estimates
should be taken as provisional until more detailed follow-up is completed.

Table 6: The occurrence rate density in two radius bins

$R_{\mathrm{min}}\,[R_{\mathrm{J}}]$ | $R_{\mathrm{max}}\,[R_{\mathrm{J}}]$ | rate densityaaKipping (2013b) | integrated ratebbIncluded in the Uehara et al. (2016) catalog.
---|---|---|---
$0.1$ | $0.4$ | $0.45\pm 0.20\,(0.36\pm 0.16)$ | $1.57\pm 0.70\,(1.26\pm 0.56)$
$0.4$ | $1.0$ | $0.18\pm 0.07\,(0.16\pm 0.06)$ | $0.42\pm 0.16\,(0.36\pm 0.14)$
$0.1$ | $1.0$ | $0.24\pm 0.07\,(0.22\pm 0.06)$ | $1.41\pm 0.41\,(1.29\pm 0.37)$

bbfootnotetext: The integrated rate is computed by integrating the rate
density over the bin. Note that the first two rows do not sum to the last row
because each row is computed assuming that the rate density is uniform across
the bin.

Note. — These values are computed in the period range 2–25 years.

##  Section 7 Astrophysical false positives

Various configurations of eclipsing binary stars can mimic the signal of a
transiting planet. However, the occurrence rate calculation presented in the
preceding section assumes no astrophysical false positives among the
candidates identified in this work. In this Section we explore the validity of
this assumption.

While an eclipsing binary (EB) typically produces a photometric dip much
deeper than a transiting planet, the depth of the signal may be comparable to
that of a planet if the eclipse is grazing, or if the EB only comprises a
small fraction of the total light in the photometric aperture – a so-called
blended eclipsing binary (BEB). Additionally, if a binary star has an
eccentric orbit, it may be oriented so as to present only a secondary
occultation and not a primary eclipse, causing a shallow and potentially flat-
bottomed photometric dip without an accompanying tell-tale deep primary
signal.

To determine to what extent the catalog of detections presented in this work
may contain such false positives, we simulate populations of detected signals
to predict how many we should expect. To accomplish this, we use the Python
package exosyspop (Morton & Foreman-Mackey, 2016), which we developed for this
purpose and utilizes the isochrones, vespa, and batman packages (Morton,
2015a, b; Kreidberg, 2015) for simulations of stellar populations and their
eclipses.

With exosyspop one can define the parameters of a population model and
generate synthetic catalogs according to the model (and the parameters of a
survey) very efficiently. For example, a population of EBs may be defined by a
binary fraction, power-law distributions in mass ratio and period (within
given bounds), and a beta distribution for eccentricity. This population,
initialized with a catalog of target stars (each of which has a duty cycle and
total span of observation), may then be “observed,” returning a catalog of
objects detectable via either primary or secondary eclipse (according to
randomly oriented orbital geometries and accounting for observation duty cycle
and data span). This synthetic catalog includes signal-to-noise estimates of
both the primary and secondary eclipses, the number of detected primary and
secondary eclipses, and the trapezoidal shape parameters of each detection
(depth, duration, and ingress-to-duration ratio, as defined in Morton 2012).

In order to predict how many EBs or BEBs we might expect to detect in this
particular search of Kepler data, we first need to choose reasonable
parameters for the binary star population. To do this, we calibrate the
population parameters using the catalog of detected Kepler eclipsing binaries.
We find that a binary fraction of 25% between periods of 20 days and 25 years
with a log-flat period distribution and eccentrities distributed according to
$\beta(0.8,2.0)$ is able to reproduce well both the number and period
distribution of observed Kepler EBs between 20 and 1000 days. We thus fix
these binary star population parameters for our subsequent EB and BEB
simulations.

To simulate synthetic populations of EB detections, we assign binary stars to
the Kepler target list described in Section 3.1 according to the above EB
population parameters. We consider an EB to be detected if it presents fewer
than three eclipses (either primary or secondary, but not both), if the
signal-to-noise ratio is $>$15, and the duration of the detected eclipse is
$<$2.5 d. In 100 realizations of these synthetic observations, we see $7.2\pm
2.5$ single- or double-eclipsing EB signals.

To simulate BEBs, we assume an exponentially varying background field star
density across the Kepler field, from 0.005 arcsec-2 at a Galactic latitude
$b=20$ to 0.05 arcsec-2 at $b=5$ (matching up well with the simulations of
Morton & Johnson (2011) at many different Galactic latitudes). Each Kepler
target star is then assigned a number of background stars drawn from a Poisson
distribution with mean given by the expected number of stars to be found
within a circle of 4 arcsec radius, given the appropriate density at its
Galactic lattitude. We draw the specific background stars from a TRILEGAL
(Girardi et al., 2005) field star simulation toward the center of the Kepler
field. Binary companions are then assigned to these background stars according
to the same stellar binary population distribution as the EB population above,
and synthetic detected populations are “observed” according to the same rules
(accounting appropriately for the diluted eclipse depths in the Kepler
bandpass). In 100 synthetic observations, we see an average of only 0.41
detected BEBs.

These prediction results suggest that we should indeed expect to see some
astrophysical false positives in our search. However, this does not mean that
we should fear that $\sim$7 of the planet candidates might be EBs. In
particular, we note that these simulations do not include the full vetting
procedure described in Section 4, and it is likely that the three impact-
parameter-rejected candidates are EBs and that the centroid-rejected candidate
is a BEB. Thus, we might expect maybe two or three additional false positives
among our planet candidates.

In order to more precisely quantify which of the candidates might indeed be
false positives, we can inspect the synthetic observation simulations in more
detail. In particular, we can analyze the shape distribution of the different
scenarios and compare them with the observed shapes of the actual Kepler
detections in order to quantify the probability that each of them may be a
false positive. These distributions and the observed shape parameters are
plotted in Figure 8.

Following the method of Morton (2012) used to compute false positive
probabilities for the regular Kepler KOI catalogs (Morton et al., 2016), we
can calculate the posterior probability for each of our candidates to belong
to each of the three scenarios we consider (EB, BEB, or planet) as follows:

$\displaystyle{\rm
Pr}_{i}=\frac{\pi_{i}\mathcal{L}_{i}}{\sum_{j}\pi_{j}\mathcal{L}_{j}},$ (20)

where the $\pi$ factors are the hypothesis priors, $\mathcal{L}$ are the
hypothesis likelihoods, and the sum over $j$ is over all the hypotheses. In
this case, we determine the relative hypothesis priors from the synthetic
observations, using the mean numbers of “observed” EBs (7.2) and BEBs (0.4),
and choosing the expected number of planets to be 12. We calculate the
hypothesis likelihoods using the depth/duration distributions of synthetic
populations of each scenario and evaluating these distributions at the
observed depths and durations of each candidate signal. To estimate the
expected shape distribution of the planet scenario, we define a custom
exosyspop population of planets according to the two-bin population model
described by the median posterior values in Table 6 and generate a population
of 1000 detected signals.

We list the probability that each candidate is a planet in Table 2. We find
that the two deepest signals in our candidate catalog (9306307 and 10602068)
are very likely to be EBs, though we note that this result may be dominated by
the fact that our planet population is fixed to have a maximum radius of 1
$R_{\rm J}$. Most of the rest of the candidates have false positive
probabilities below 10%. We do note that as discussed in Section 4, 4754460
(for which we calculate a 5% false positive probability) does show a partial
deep eclipse right at the end of the Kepler data that indicates that it is
most likely an EB. Apart from this, the expected number of false positives
among the candidates with $R<R_{\rm J}$, according to these calculations, is
about one.

In light of this result, we demonstrate the sensitivity of our measured
occurrence rates on contamination by computing a second constraint on
$\Gamma_{V}$ for each volume with one candidate removed. In this case,
Equation (16) would be replaced by

$\displaystyle\tilde{\Gamma}_{V}$ $\displaystyle\equiv$
$\displaystyle\frac{C(P_{\mathrm{min}},\,P_{\mathrm{max}};\,R_{\mathrm{min}},\,R_{\mathrm{max}})-1}{Z(P_{\mathrm{min}},\,P_{\mathrm{max}};\,R_{\mathrm{min}},\,R_{\mathrm{max}})}\quad.$
(21)

These updated rates are listed in Table 6. In each case, the results are
consistent within the uncertainties but the difference can be used to get a
qualitative sense of the systematic uncertainty introduced by the false
positive population.

We note that in the above procedure we have not corrected our predictions for
the fact that our search has explicitly excluded KOIs that host known Kepler
EBs—if any of these excluded systems show fewer than three eclipses and do not
present both primary and secondary eclipses, then they should also should
perhaps count towards the number of EBs we should have expected to find in
this survey. However, as the Kepler EB catalog does not provide information on
whether both primary and secondary eclipses are detected, we neglect this
correction. We note that this is a conservative decision, in the sense that
accounting for the effect of excluded EBs on our predictions would only
further decrease the FPP of the planetary signals, as they would be even less
likely to be caused by EBs.

Figure 8: Predicted eclipse shape distributions for the two false positive
scenerios and exoplanet transits (grayscale heatmap). In this figure, the
relative normalization of the maps are arbitrary but the absolute
normalization is discussed in Section 7. The green points show the shape
parameters of the long-period exoplanet candidates from Table 2.

##  Section 8 Comparison with the literature

The population of long-period planets has been previously studied using radial
velocity, microlensing, and direct imaging surveys. These methods all measure
the occurrence rate as a function of planet mass instead of radius. Using the
transit method, however, we do not directly measure the mass of the planet.
Therefore, to compare our results with the literature, we must rely on a
mass–radius (M–R) relationship constructed using exoplanets with both mass and
radius measurements (for example Weiss & Marcy, 2014; Wolfgang et al., 2016;
Chen & Kipping, 2016) to predict the expected masses of the transiting
planets.

Table 7 lists gives constraints on the predicted masses of the exoplanet
candidates using the probabilistic M–R relationship from Chen & Kipping (2016)
and taking uncertainties in the planet radius and statistical uncertainties in
the M–R parameters. We compare the predictions with the predictions from
Wolfgang et al. (2016) and find similar values with smaller uncertainties and
choose to use the Chen & Kipping (2016) relationship because it is more
conservative in the relevant range of parameter space.

A detailed discussion of the systematic effects introduced by the use of a M–R
relationship is beyond the scope of this paper but it is worth noting that all
published relationships are based on exoplanets much closer to their host star
than any of the candidates discussed here. This effect would cause the masses
of these cool planets to be systematically underestimated.

Using the same M–R relationship, we also compute the completeness of our
transit search as a function of planet mass and orbital semi-major axis. This
function is plotted in Figure 9 with the same color scheme as Figure 6. These
injections and the predicted masses and measured semi-major axes of the
candidates can then be used to estimate the occurrence rate in mass–semi-major
axis units using the method from Section 6. One small change to Equation (18)
is necessary to account for the fact that the injections were not made
uniformly in $\ln M$ and $\ln a$. We numerically estimate the prior
distribution in mass and semi-major axis from which the injections were drawn
$\tilde{p}(\ln a,\,\ln M)$ and Equation (18) becomes

$\displaystyle
Z(a_{\mathrm{min}},\,a_{\mathrm{max}};\,M_{\mathrm{min}},\,M_{\mathrm{max}})$
$\displaystyle\approx$
$\displaystyle\frac{K\,V}{J}\,\sum_{j=1}^{J}\frac{Q_{k}(w^{(j)})}{\tilde{p}(\ln
a^{(j)},\,\ln M^{(j)})}\quad.$ (22)

Using this result, we find that the mean occurrence rate density in the range
$0.01\,M_{\mathrm{J}}\leq M<20\,M_{\mathrm{J}}$ and $1.5{\,\mathrm{au}}\leq
a<9{\,\mathrm{au}}$ is

$\displaystyle\frac{\,\mathrm{d}^{2}N}{\,\mathrm{d}\ln M\,\,\mathrm{d}\ln a}$
$\displaystyle=$ $\displaystyle 0.068\pm 0.019$ (23)

where $N$ is the expected number of planets per G/K dwarf. This result and the
equivalent result as a function of planet mass and orbital period are listed
in Table 8.

The uncertainty in Equation (23) and Table 8 does not take into account the
uncertainties in the mass estimates or any systematic noise in the mass–radius
relationship. Therefore, these specific results should be taken with the
appropriate grain of salt but predictions in these parameter spaces ease
comparison with occurrence rates measured computed using different methods.

Clanton & Gaudi (2016) studied the occurrence rate of long-period giant
planets orbiting M dwarfs by combining results from radial velocity,
microlensing, and direct imaging surveys. In the period and mass range
$10^{3}-10^{4}{\,\mathrm{days}}$ and $10-10^{4}\,M_{\oplus}$, they find a mean
occurrence rate density of

$\displaystyle\frac{\,\mathrm{d}^{2}N}{\,\mathrm{d}\ln M\,\,\mathrm{d}\ln P}$
$\displaystyle=$ $\displaystyle 0.023$ (24)

per M dwarf with large uncertainty. This result is slightly lower than our
estimated rate for a similar range of masses and periods but around G/K
dwarfs. This difference is consistent with previous observational and
theoretical results that cooler stars host fewer long-period giant planets
(for example Laughlin et al., 2004; Cumming et al., 2008; Clanton & Gaudi,
2016).

Recently, Bryan et al. (2016) studied the frequency of long-period giant
planets in systems with inner hot Jupiters based on long-baseline radial
velocity monitoring of these systems (Knutson et al., 2014). In this sample,
the computed occurrence rate of long-period giant planets was found to be

$\displaystyle\frac{\,\mathrm{d}^{2}N}{\,\mathrm{d}\ln M\,\,\mathrm{d}\ln a}$
$\displaystyle=$ $\displaystyle 0.125\pm 0.012$ (25)

in the range $1-20\,M_{\mathrm{J}}$ and $5-20{\,\mathrm{au}}$. This result is
about a factor of two larger than our estimate (Equation 23) once again
suggesting that cold Jupiters might preferentially occur in systems with inner
planets – or that the presence of cold Jupiters encourages the formation of
hot Jupiters.

A recent review of the occurrence rate estimates based on direct imaging
surveys (Bowler, 2016) reports the upper limit on the occurrence rate of giant
planets orbiting F/G/K dwarfs as $<6.8\%$ in the range $5-13\,M_{\mathrm{J}}$
and $10-100{\,\mathrm{au}}$. Converted to a rate density, this gives

$\displaystyle\frac{\,\mathrm{d}^{2}N}{\,\mathrm{d}\ln M\,\,\mathrm{d}\ln a}$
$\displaystyle<$ $\displaystyle 0.03\quad.$ (26)

This value is lower than the value computed using our sample in Table 8 but
this is consistent with the fact that direct imaging is sensitive to the
potentially less common large planets at wider separations than detections
with the transit method.

As a final comparison, we repeated the analysis of Burke et al. (2015) and fit
a double power-law occurrence rate to the short-period Kepler planet
candidates777The analysis was adapted from publicly available code that was
demonstrated to reproduce the same results as Burke et al. (2015) by Foreman-
Mackey (2015). and extrapolated to the center of the two bins where we
computed the occurrence rate. At a period of 7 years and a radius of
$0.2\,R_{\mathrm{J}}$, the extrapolated occurrence rate density is $0.73\pm
0.28$ and at a radius of $0.6\,R_{\mathrm{J}}$, the extrapolated rate density
is $0.15\pm 0.05$. These extrapolated values are qualitatively consistent with
the rates listed in Table 6 but we note that extrapolations and their
statistical uncertainties should not be taken too seriously.

Table 7: The predicted masses and radial velocity semi-amplitudes for the
candidates from Table 2

kic id | Kp | radius | mass | period | $t_{0}$ | semi-major axis | semi-amplitude | $K/P$
---|---|---|---|---|---|---|---|---
|  | $R_{\mathrm{J}}$ | $M_{\mathrm{J}}$ | year | KBJD | au | m s-1 | m s-1 yr-1
3218908 | $14.6$ | $0.514_{-0.093}^{+0.092}$ | $0.079_{-0.038}^{+0.074}$ | $7.0_{-3.4}^{+9.5}$ | $766.6722_{-0.0114}^{+0.0096}$ | $3.4_{-1.2}^{+2.7}$ | $1.41_{-0.77}^{+1.57}$ | $0.20_{-0.15}^{+0.49}$
3239945 | $14.0$ | $0.876_{-0.039}^{+0.039}$ | $6.5_{-6.3}^{+38.1}$ | $2.9328721_{-0.0000026}^{+0.0000026}$ | $420.28714_{-0.00068}^{+0.00069}$ | $1.864_{-0.023}^{+0.025}$ | $157.5_{-153.0}^{+900.5}$ | $53.7_{-52.2}^{+307.1}$
4754460 | $14.9$ | $0.67_{-0.15}^{+0.16}$ | $0.140_{-0.078}^{+9.289}$ | $5.9_{-3.0}^{+11.8}$ | $826.8369_{-0.0046}^{+0.0046}$ | $3.1_{-1.2}^{+3.4}$ | $2.5_{-1.5}^{+185.3}$ | $0.46_{-0.39}^{+35.87}$
6551440 | $13.6$ | $0.282_{-0.083}^{+0.093}$ | $0.028_{-0.015}^{+0.033}$ | $4.0_{-1.2}^{+4.2}$ | $1039.0589_{-0.0037}^{+0.0037}$ | $2.50_{-0.52}^{+1.55}$ | $0.57_{-0.33}^{+0.77}$ | $0.133_{-0.094}^{+0.230}$
8410697 | $13.4$ | $0.698_{-0.078}^{+0.107}$ | $0.157_{-0.081}^{+11.382}$ | $2.8688097_{-0.0000054}^{+0.0000053}$ | $542.1231_{-0.0013}^{+0.0013}$ | $1.925_{-0.038}^{+0.054}$ | $3.6_{-1.9}^{+262.6}$ | $1.26_{-0.65}^{+91.55}$
8426957 | $13.6$ | $1.04_{-0.25}^{+0.30}$ | $3.8_{-3.6}^{+47.5}$ | $54.2_{-36.1}^{+88.4}$ | $784.677_{-0.013}^{+0.013}$ | $14.1_{-7.4}^{+13.0}$ | $33.8_{-32.4}^{+389.5}$ | $0.73_{-0.71}^{+10.01}$
8505215 | $13.0$ | $0.277_{-0.017}^{+0.017}$ | $0.028_{-0.012}^{+0.022}$ | $9.1_{-3.4}^{+9.5}$ | $140.0492_{-0.0018}^{+0.0017}$ | $4.0_{-1.1}^{+2.5}$ | $0.45_{-0.21}^{+0.39}$ | $0.048_{-0.032}^{+0.070}$
8738735 | $13.9$ | $0.355_{-0.044}^{+0.045}$ | $0.042_{-0.019}^{+0.037}$ | $9.9_{-5.0}^{+14.9}$ | $697.8538_{-0.0049}^{+0.0059}$ | $4.8_{-1.8}^{+4.0}$ | $0.54_{-0.27}^{+0.54}$ | $0.053_{-0.040}^{+0.127}$
8800954 | $13.4$ | $0.386_{-0.025}^{+0.025}$ | $0.049_{-0.022}^{+0.039}$ | $1.9279957_{-0.0000091}^{+0.0000092}$ | $492.7652_{-0.0024}^{+0.0024}$ | $1.420_{-0.028}^{+0.026}$ | $1.39_{-0.62}^{+1.11}$ | $0.72_{-0.32}^{+0.58}$
9306307 | $14.0$ | $1.22_{-0.36}^{+0.49}$ | $4.6_{-4.3}^{+97.9}$ | $4.3_{-1.1}^{+3.3}$ | $1191.35648_{-0.00018}^{+0.00018}$ | $2.39_{-0.45}^{+1.09}$ | $106.0_{-99.4}^{+2186.4}$ | $22.4_{-21.1}^{+488.4}$
10187159 | $14.4$ | $0.43_{-0.13}^{+0.21}$ | $0.061_{-0.036}^{+0.096}$ | $4.9_{-1.8}^{+7.6}$ | $604.1102_{-0.0031}^{+0.0023}$ | $2.70_{-0.69}^{+2.35}$ | $1.25_{-0.82}^{+2.42}$ | $0.26_{-0.21}^{+0.69}$
10287723 | $13.4$ | $0.266_{-0.024}^{+0.027}$ | $0.026_{-0.012}^{+0.022}$ | $4.9_{-1.3}^{+4.3}$ | $393.5976_{-0.0029}^{+0.0031}$ | $2.58_{-0.46}^{+1.35}$ | $0.62_{-0.28}^{+0.52}$ | $0.117_{-0.068}^{+0.129}$
10321319 | $11.9$ | $0.163_{-0.037}^{+0.046}$ | $0.0120_{-0.0058}^{+0.0115}$ | $5.5_{-2.1}^{+8.1}$ | $554.3562_{-0.0063}^{+0.0064}$ | $2.93_{-0.81}^{+2.44}$ | $0.23_{-0.12}^{+0.25}$ | $0.041_{-0.031}^{+0.069}$
10602068 | $14.9$ | $2.00_{-0.35}^{+0.66}$ | $162.4_{-161.6}^{+73.9}$ | $3.16_{-0.83}^{+2.65}$ | $830.80892_{-0.00015}^{+0.00015}$ | $2.11_{-0.40}^{+1.09}$ | $3100.4_{-3083.7}^{+1797.2}$ | $889.4_{-885.4}^{+890.8}$
10842718 | $14.6$ | $0.74_{-0.16}^{+0.16}$ | $0.19_{-0.12}^{+22.15}$ | $12.7_{-6.6}^{+20.2}$ | $226.2344_{-0.0047}^{+0.0047}$ | $5.3_{-2.1}^{+4.7}$ | $2.6_{-1.7}^{+324.0}$ | $0.25_{-0.22}^{+30.41}$
11709124 | $14.5$ | $0.83_{-0.11}^{+0.12}$ | $0.93_{-0.82}^{+35.44}$ | $4.3_{-1.3}^{+4.7}$ | $657.2674_{-0.0016}^{+0.0018}$ | $2.54_{-0.56}^{+1.62}$ | $18.1_{-16.0}^{+668.5}$ | $3.5_{-3.1}^{+146.2}$

Note. — The masses and radial velocity amplitudes are estimated based on the
measured radius using a probabilistic mass–radius relation (Chen & Kipping,
2016).

Figure 9: The same as Figure 6 converted into the planet mass and semi-major
axis plane. Since the completeness function depends on the planet’s radius and
not its mass, a probabilistic mass–radius relationship (Chen & Kipping, 2016)
was used to convert radius to mass.

Table 8: The occurrence rate computed in mass units

volume | rate densityaaThe rate density is given by Equation (16) and the value in parentheses is computed assuming one candidate is a false positive (Equation 21). | integrated rate
---|---|---
$2{\,\mathrm{yr}}<P<25{\,\mathrm{yr}};\,0.01\,M_{\mathrm{J}}\leq M<20\,M_{\mathrm{J}}$ | $0.046\pm 0.013$ | $0.882\pm 0.245$
$1.5{\,\mathrm{au}}<a<9{\,\mathrm{au}};\,0.01\,M_{\mathrm{J}}\leq M<20\,M_{\mathrm{J}}$ | $0.068\pm 0.019$ | $0.925\pm 0.257$

aafootnotetext: The rate density is measured in natural logarithmic units; see
Equation (23).

Note. — These values are computed assuming that the occurrence rate is flat in
the logarithmic parameters.

##  Section 9 Prospects for follow-up

A real concern about the detection of exoplanets from a single transit is that
follow-up and confirmation is difficult. Since the period of the orbit is
poorly constrained and transits are sparse, any prediction of a subsequent
transit time will be too uncertain to schedule targeted photometric follow-up
(Beichman et al., 2016; Dalba & Muirhead, 2016). Instead, follow-up using
radial velocity and astrometry are more promising. For both radial velocity
and astrometry, there is information about the orbiting planets in
measurements made at all times – not just during transit. This allows
observations to be scheduled without a well constrained orbital period.
Furthermore, follow-up of any of these candidates using radial velocity or
astrometry would provide a measurement of the density of a planet that would
be valuable for the study of planetary compositions.

Table 7 lists the posterior predictions for the semi-amplitude $K$ of the
radial velocity signal produced by each candidate using the mass predictions
from the previous section. Since the orbital periods are long, we also include
a simple prediction for the slope of the radial velocity trend induced by this
planet by taking the ratio of the semi-amplitude and the orbital period. Any
radial velocity follow-up of the candidates presented here would be an
ambitious undertaking because the stars are relatively faint and, in most
cases, the radial velocity trends are small. Some candidates should, however,
be within the reach of current state-of-the-art facilities.

In principle the Gaia Mission will be very sensitive to the astrometric wobble
produced by a long-period exoplanet (Perryman et al., 2014). To leading order,
the astrometric signal strength is proportional to the semi-major axis of the
stellar (primary) reflex motion in angular units. That is, detectability is
related to the angle $\alpha$ given by

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{a}{D}\,\frac{M_{\mathrm{p}}}{M_{\mathrm{s}}}\quad,$ (27)

where $a$ is the semi-major axis, $D$ is the distance from the observer to the
primary, and $M_{\mathrm{p}}/M_{\mathrm{s}}$ is the planet-to-star mass ratio.

The single-visit precision of Gaia will vary with magnitude but is expected to
be on the order of 40 $\mu$as at these magnitudes. In detail the confidence
with which an exoplanet can be detected or measured in the final Gaia data
depends on this precision, the number of crossings of the star through the
Gaia field-of-view, and details of how the projected orbit is sampled by the
time history of the focal-plane crossings. However, it is not expected that
Gaia can detect or precisely measure exoplanet-induced astrometric wobbles
that are much smaller in amplitude than the single-visit precision (Perryman
et al., 2014).

The primary stars in the Kepler Field are typically at distances of $\sim 500$
pc, and typical mass ratios are in the $10^{-4}$ range. We therefore expect
astrometric amplitudes in the 0.3 to 3 $\mu$as range. These planets will not
be detectable or measurable in the Gaia data under any circumstances, but it
may be possible to identify which candidates are false-alarms caused by
eclipsing binaries. However, similar planets around closer stars _will_ be
detectable with Gaia. This means that there will be a comparable exoplanet
occurrence rate measurement from the Gaia data. It also means that many of the
discoveries of the K2 and TESS missions could be followed up and precisely
measured by the Gaia Mission.

##  Section 10 Summary

We have developed a fully automated method to search for the transits of long-
period exoplanets with only one or two observable transits in the Kepler
archival light curves. This method uses probabilistic model comparison to veto
non-transit signals. Applying this method to the brightest 39,036 G/K dwarfs
in the Kepler target list, we discover 16 systems with likely astrophysical
transits and eclipses. We fit the light curve for each candidate with a
physical generative model and informative priors on eccentricity and stellar
density to estimate the planet’s orbital period. The constraint on the period
is also informed by the simplifying assumption that no other transit could
occur during the baseline of Kepler observations of the target. Simulations of
the false positive population – lone primary or secondary eclipses of binary
systems or background eclipsing binaries – suggest that 13 of these candidates
have high probability of being planetary in nature.

We measure the empirical detection efficiency function of our search procedure
by injecting simulated transit signals into the target light curves and
measuring the recovery rate. By combining the measured detection efficiency
and the catalog of exoplanet candidates, we estimate the integrated occurrence
rate of exoplanets with orbital periods in the range $2-25{\,\mathrm{years}}$
and radii in the range $0.1-1\,R_{\mathrm{J}}$ to be $2.00\pm 0.72$ planets
per G/K dwarf. This result is qualitatively consistent with estimates of the
occurrence rate of long-period giant planets based on data from radial
velocity and direct imaging surveys. The occurrence rate measured here – for
Sun-like hosts – is higher than microlensing results for generally lower mass
stars (Gaudi, 2012; Clanton & Gaudi, 2014, 2016) but this discrepancy is
consistent with predictions from the core-accretion model (Laughlin et al.,
2004).

Using a probabilistic mass–radius relationship, we predict the masses of our
candidates and report predictions for the radial velocity semi-amplitudes.
Unfortunately, since the target stars are faint and the amplitudes are small,
these targets are unlikely to be accessible with even the current state-of-
the-art high-precision instruments. We also discuss the potential for
astrometric follow-up using the forthcoming data from the Gaia Mission with
similarly discouraging results.

Any detailed analysis of individual systems detected with only a single
transit requires follow-up observations to convincingly rule out false
positive scenarios and to better characterize the stellar host parameters
(with, for example, parallax measurements from Gaia). The conclusions of this
work – and all other occurrence rate results based on Kepler data – are
conditioned on the assumption that the stellar characterization of the target
sample is systematic and un-biased. The main population-level results should
be fairly insensitive systematic issues with the sample but a rigorous
analysis of these effects will be required to come to more detailed
conclusions about this population of long-period transiting planets.

Our method of transit discovery is especially relevant for future photometric
surveys like K2, TESS, and PLATO where the survey baseline is shorter than
Kepler. The transits of planets with orbital periods longer than the
observation baselines will be plentiful in these forthcoming data sets and
this method can, in principle, be trivially generalized to discover these
planets, prioritize follow-up, and study their population.

All of the code used in this project is available from
https://github.com/dfm/peerless under the MIT open-source software license.
This code (plus some dependencies) can be run to re-generate all of the
figures and results in this paper; this version of the paper was generated
with git commit a485057 (2016-09-30). The parameter estimation results
represented as MCMC samplings and the injection results are available for
download from Zenodo at http://dx.doi.org/10.5281/zenodo.58273.

It is a pleasure to thank Jeff Coughlin, So Hattori, Heather Knutson, Phil
Muirhead, Darin Ragozzine, Hans-Walter Rix, Dun Wang, and Angie Wolfgang for
helpful discussions and contributions. We thank the anonymous referee for
comments that improved the presentation and clarity of this manuscript.

T.D.M. was supported by the National Aeronautics and Space Administration,
under the Kepler participating scientist program (grant NNX14AE11G), and is
grateful for the hospitality of both the Institute for Advanced Study and
Carnegie Observatories that helped support this work. D.W.H. was partially
supported by the National Science Foundation (grant IIS-1124794), the National
Aeronautics and Space Administration (grant NNX12AI50G), and the Moore–Sloan
Data Science Environment at NYU. E.A. acknowledges support from NASA grants
NNX13AF20G, NNX13AF62G, and NASA Astrobiology Institutes Virtual Planetary
Laboratory, supported by NASA under cooperative agreement NNH05ZDA001C.

This research made use of the NASA Astrophysics Data System and the NASA
Exoplanet Archive. The Exoplanet Archive is operated by the California
Institute of Technology, under contract with NASA under the Exoplanet
Exploration Program.

This paper includes data collected by the Kepler mission. Funding for the
Kepler mission is provided by the NASA Science Mission directorate. We are
grateful to the entire Kepler team, past and present. Their tireless efforts
were all essential to the tremendous success of the mission and the successes
of K2, present and future.

These data were obtained from the Mikulski Archive for Space Telescopes
(MAST). STScI is operated by the Association of Universities for Research in
Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST is provided
by the NASA Office of Space Science via grant NNX13AC07G and by other grants
and contracts.

Computing resources were provided by High Performance Computing at New York
University.

Kepler .

##  Appendix A Details of the light curve models

In Section 2.2, the five light curve models were listed. In this section, we
give the mathematical details of each model and list the parameters that are
fit. Each model – except the transit model – can be easily differentiated with
respect to its parameters. As discussed in the following section, this feature
is crucial for efficient and robust likelihood maximization.

  * •

The box model is given by

$\displaystyle m_{\mathrm{box}}(t)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}a,&\mbox{if $t\leq
t_{\mathrm{min}}$}\\\ b,&\mbox{if $t_{\mathrm{min}}<t\leq
t_{\mathrm{max}}$}\\\ c,&\mbox{if $t_{\mathrm{max}}<t$}\end{array}\right.$
(A4)

where $a$, $b$, and $c$ are free parameters, and $t_{\mathrm{min}}$ and
$t_{\mathrm{max}}$ are fixed. In practice, we include two different box models
where $t_{\mathrm{min}}$ and $t_{\mathrm{max}}$ are set using different
heuristics. The first box has the bounds set to match the ingress and egress
of the best fit transit. The second box is chosen based on the largest change
points in the light curve.

  * •

The step model is given by

$\displaystyle m_{\mathrm{step}}(t)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}m_{1}+h_{1}\,\exp\left([t-t_{0}]/w_{1}\right),&\mbox{if
$t<t_{0}$}\\\ m_{2}+h_{2}\,\exp\left([t_{0}-t]/w_{2}\right),&\mbox{if
$t_{0}\leq t$}\end{array}\right.$ (A7)

where all of the parameters – including $t_{0}$ – are included in the fit. To
ensure that the widths $w_{1}$ and $w_{2}$ remain positive, we fit for $\log
w_{1}$ and $\log w_{2}$.

  * •

For the outlier model, we iterate through all cadences $t_{n}$ within 0.3 days
of the candidate transit time and evaluate the model as

$\displaystyle m_{\mathrm{outlier}}(t)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}f(t_{0}),&\mbox{if $t=t_{0}$}\\\
\mathrm{median}[f(t\neq t_{0})],&\mbox{if $t\neq t_{0}$}\end{array}\right.$
(A10)

where $f(t)$ is the observed time series. With this model, no non-linear
optimization is required and the final value of $t_{0}$ is the one with the
maximum likelihood in this grid search.

  * •

The variability model only has one parameter, the flux $m_{0}$ and
$m_{\mathrm{variability}}(t)=m_{0}$ at all times. The variability is captured
by the Gaussian Process residual model.

  * •

Finally, the transit model is an exposure time integrated, limb darkened light
curve (Mandel & Agol, 2002; Kipping, 2010) parameterized by the radius ratio
between the planet and star, the transit duration, the transit time, the
impact parameter, and two quadratic limb darkening coefficients (Kipping,
2013a). Analytically computing the gradient of a simple transit model is
possible (Pál, 2008) but it becomes substantially more tedious as the model
becomes more realistic. Therefore, we instead use a compile-time automatic
differentiation library888More specifically, we use the Jet object from the
BSD-licensed Ceres Solver http://ceres-solver.org (Agarwal et al., 2016) to
efficiently compute first derivatives of the full transit model with respect
to the orbital and physical parameters to machine precision.

##  Appendix B Gaussian process regression

Gaussian Processes (GPs) are a class of non-parametric, stochastic models that
have been demonstrated to be good effective models for the variability in
Kepler light curves. A simple GP model can be used to capture residual non-
transit variability in light curves. In this paper, we use a GP model for two
steps: light curve–level transit shape vetting and parameter estimation. A
full discussion of GPs is beyond the scope of this paper, so we will only
summarize the most relevant points here and direct an interested reader to
Rasmussen & Williams (2006) for more details.

A GP model is specified by the following likelihood function

$\displaystyle\mathcal{L}=\ln
p({\boldsymbol{y}}\,|\,{{\boldsymbol{\theta}}},\,{{\boldsymbol{\alpha}}})$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\,{\boldsymbol{r}}({{\boldsymbol{\theta}}})^{\mathrm{T}}\,K({{\boldsymbol{\alpha}}})^{-1}\,{\boldsymbol{r}}({{\boldsymbol{\theta}}})-\frac{1}{2}\log\det
K({{\boldsymbol{\alpha}}})-\frac{N}{2}\log{2\,\pi}$ (B1)

where $\boldsymbol{y}$ is a list of measurements in a scalar time series – in
this case, fluxes – measured at the times $\boldsymbol{t}$, and

$\displaystyle{\boldsymbol{r}}({{\boldsymbol{\theta}}})$ $\displaystyle=$
$\displaystyle{\boldsymbol{y}}-m({\boldsymbol{t}};\,{{\boldsymbol{\theta}}})$
(B2)

is the vector of residuals away from the mean model
$m({\boldsymbol{t}};\,{{\boldsymbol{\theta}}})$. For the purposes of this
paper, we model the covariance matrix $K({{\boldsymbol{\alpha}}})$ using the
Matérn-3/2 kernel. Under this model, the elements of
$K({{\boldsymbol{\alpha}}})$ are given by

$\displaystyle\left[K({{\boldsymbol{\alpha}}})\right]_{ij}$ $\displaystyle=$
$\displaystyle\sigma_{i}^{2}\,\delta_{ij}+\alpha^{2}\left[1+\frac{|t_{i}-t_{j}|}{\sqrt{3}\,\tau}\right]\exp\left(-\frac{|t_{i}-t_{j}|}{\sqrt{3}\,\tau}\right)$
(B3)

where $\sigma_{i}$ is the reported uncertainty on the $i$-th measurement in
the time series and $\delta_{ij}$ is the Kronecker delta.

This covariance function (Equation B3) is specified by an amplitude $\alpha$
and a time scale $\tau$ and we will simultaneously fit for these
hyperparameters ${{\boldsymbol{\alpha}}}=(\alpha,\,\tau)$ and the parameters
of the mean model ${\boldsymbol{\theta}}$. To efficiently find the parameter
set that maximizes Equation (B1) using a non-linear optimization routine999We
use the L-BFGS-B method as implemented in SciPy
http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html.,
it is useful to be able to compute the gradient of Equation (B1) with respect
to the parameters ${\boldsymbol{\theta}}$ and ${\boldsymbol{\alpha}}$. These
gradients are given by

$\displaystyle\frac{\,\mathrm{d}\ln
p({\boldsymbol{y}}\,|\,{{\boldsymbol{\theta}}},\,{{\boldsymbol{\alpha}}})}{\,\mathrm{d}{{\boldsymbol{\theta}}}}$
$\displaystyle=$
$\displaystyle\frac{\,\mathrm{d}m({\boldsymbol{t}};\,{{\boldsymbol{\theta}}})}{\,\mathrm{d}{{\boldsymbol{\theta}}}}^{\mathrm{T}}\,K({{\boldsymbol{\alpha}}})^{-1}\,{\boldsymbol{r}}({{\boldsymbol{\theta}}})$
(B4)

and

$\displaystyle\frac{\,\mathrm{d}\ln
p({\boldsymbol{y}}\,|\,{{\boldsymbol{\theta}}},\,{{\boldsymbol{\alpha}}})}{\,\mathrm{d}{{\boldsymbol{\alpha}}}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\,\mathrm{Tr}\left(\left[{\boldsymbol{\phi}}\,{\boldsymbol{\phi}}^{\mathrm{T}}-K({{\boldsymbol{\alpha}}})^{-1}\right]\,\frac{\,\mathrm{d}K({{\boldsymbol{\alpha}}})}{\,\mathrm{d}{{\boldsymbol{\alpha}}}}\right)$
(B5)

where

$\displaystyle{\boldsymbol{\phi}}$ $\displaystyle=$ $\displaystyle
K({{\boldsymbol{\alpha}}})^{-1}\,{\boldsymbol{r}}({{\boldsymbol{\theta}}})\quad.$
(B6)

## References

  * Agarwal et al. (2016) Agarwal, S., Mierle, K., & et al. 2016, Ceres Solver, http://ceres-solver.org, ,
  * Ambikasaran et al. (2016) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D., & O Neil, M. 2016, IEEE
  * Ballard & Johnson (2016) Ballard, S., & Johnson, J. A. 2016, ApJ, 816, 66
  * Beichman et al. (2016) Beichman, C., Livingston, J., Werner, M., et al. 2016, ApJ, 822, 39
  * Borucki et al. (2011) Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19
  * Bowler (2016) Bowler, B. P. 2016, ArXiv e-prints, arXiv:1605.02731
  * Bryan et al. (2016) Bryan, M. L., Knutson, H. A., Howard, A. W., et al. 2016, ApJ, 821, 89
  * Bryson et al. (2013) Bryson, S. T., Jenkins, J. M., Gilliland, R. L., et al. 2013, PASP, 125, 889
  * Burke & McCullough (2014) Burke, C. J., & McCullough, P. R. 2014, ApJ, 792, 79
  * Burke et al. (2014) Burke, C. J., Bryson, S. T., Mullally, F., et al. 2014, ApJS, 210, 19
  * Burke et al. (2015) Burke, C. J., Christiansen, J. L., Mullally, F., et al. 2015, ApJ, 809, 8
  * Cassan et al. (2012) Cassan, A., Kubas, D., Beaulieu, J.-P., et al. 2012, Nature, 481, 167
  * Chen & Kipping (2016) Chen, J., & Kipping, D. M. 2016, ArXiv e-prints, arXiv:1603.08614
  * Christiansen et al. (2013) Christiansen, J. L., Clarke, B. D., Burke, C. J., et al. 2013, ApJS, 207, 35
  * Christiansen et al. (2015) —. 2015, ApJ, 810, 95
  * Clanton & Gaudi (2014) Clanton, C., & Gaudi, B. S. 2014, ApJ, 791, 91
  * Clanton & Gaudi (2016) —. 2016, ApJ, 819, 125
  * Coughlin et al. (2016) Coughlin, J. L., Mullally, F., Thompson, S. E., et al. 2016, ApJS, 224, 12
  * Cumming et al. (2008) Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531
  * Dalba & Muirhead (2016) Dalba, P. A., & Muirhead, P. S. 2016, ApJ, 826, L7
  * Dalba et al. (2015) Dalba, P. A., Muirhead, P. S., Fortney, J. J., et al. 2015, ApJ, 814, 154
  * Dawson & Johnson (2012) Dawson, R. I., & Johnson, J. A. 2012, ApJ, 756, 122
  * Dressing & Charbonneau (2015) Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45
  * Fang & Margot (2012) Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92
  * Fischer et al. (2012) Fischer, D. A., Schwamb, M. E., Schawinski, K., et al. 2012, MNRAS, 419, 2900
  * Foreman-Mackey (2015) Foreman-Mackey, D. 2015, An experiment in open science: exoplanet population inference, http://dfm.io/posts/exopop/, Zenodo, doi:10.5281/zenodo.57218
  * Foreman-Mackey (2016) —. 2016, The Journal of Open Source Software, 24, doi:10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64
  * Foreman-Mackey et al. (2015) Foreman-Mackey, D., Montet, B. T., Hogg, D. W., et al. 2015, ApJ, 806, 215
  * Foreman-Mackey & Morton (2016) Foreman-Mackey, D., & Morton, T. 2016, dfm/transit: v0.3.0, , , doi:10.5281/zenodo.159478
  * Gaidos & Mann (2012) Gaidos, E., & Mann, A. W. 2012, ApJ, 762, 41
  * Gaudi (2012) Gaudi, B. S. 2012, ARA&A, 50, 411
  * Girardi et al. (2005) Girardi, L., Groenewegen, M. A. T., Hatziminaoglou, E., & da Costa, L. 2005, A&A, 436, 895
  * Gould et al. (2004) Gould, A., Gaudi, B. S., & Han, C. 2004, ArXiv Astrophysics e-prints, astro-ph/0405217
  * Gould et al. (2010) Gould, A., Dong, S., Gaudi, B. S., et al. 2010, ApJ, 720, 1073
  * Howard et al. (2012) Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398
  * Huber et al. (2014) Huber, D., Silva Aguirre, V., Matthews, J. M., et al. 2014, ApJS, 211, 2
  * Hunter et al. (2007) Hunter, J. D., et al. 2007, Computing in science and engineering, 9, 90
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, ,
  * Kipping (2010) Kipping, D. M. 2010, MNRAS, 408, 1758
  * Kipping (2013a) —. 2013a, MNRAS, 435, 2152
  * Kipping (2013b) —. 2013b, MNRAS, 434, L51
  * Kipping et al. (2014a) Kipping, D. M., Bastien, F. A., Stassun, K. G., et al. 2014a, ApJ, 785, L32
  * Kipping et al. (2012) Kipping, D. M., Dunn, W. R., Jasinski, J. M., & Manthri, V. P. 2012, MNRAS, 421, 1166
  * Kipping et al. (2014b) Kipping, D. M., Torres, G., Buchhave, L. A., et al. 2014b, ApJ, 795, 25
  * Kipping et al. (2016) Kipping, D. M., Torres, G., Henze, C., et al. 2016, ApJ, 820, 112
  * Kirk et al. (2016) Kirk, B., Conroy, K., Prša, A., et al. 2016, AJ, 151, 68
  * Knutson et al. (2014) Knutson, H. A., Fulton, B. J., Montet, B. T., et al. 2014, ApJ, 785, 126
  * Kovács et al. (2002) Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161
  * Laughlin et al. (2004) Laughlin, G., Bodenheimer, P., & Adams, F. C. 2004, ApJ, 612, L73
  * Mandel & Agol (2002) Mandel, K., & Agol, E. 2002, ApJ, 580, L171
  * Moriarty & Ballard (2015) Moriarty, J., & Ballard, S. 2015, ArXiv e-prints, arXiv:1512.03445
  * Morton & Foreman-Mackey (2016) Morton, T., & Foreman-Mackey, D. 2016, timothydmorton/exosyspop: v0.1-alpha, , , doi:10.5281/zenodo.159487
  * Morton (2012) Morton, T. D. 2012, ApJ, 761, 6
  * Morton (2015a) —. 2015a, isochrones: Stellar model grid package, Astrophysics Source Code Library, , , ascl:1503.010
  * Morton (2015b) —. 2015b, VESPA: False positive probabilities calculator, Astrophysics Source Code Library, , , ascl:1503.011
  * Morton et al. (2016) Morton, T. D., Bryson, S. T., Coughlin, J. L., et al. 2016, ApJ, 822, 86
  * Morton & Johnson (2011) Morton, T. D., & Johnson, J. A. 2011, ApJ, 738, 170
  * Mulders et al. (2015) Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 814, 130
  * Mullally et al. (2016) Mullally, F., Coughlin, J. L., Thompson, S. E., et al. 2016, PASP, 128, 074502
  * Osborn et al. (2016) Osborn, H. P., Armstrong, D. J., Brown, D. J. A., et al. 2016, MNRAS, 457, 2273
  * Pál (2008) Pál, A. 2008, MNRAS, 390, 281
  * Perryman et al. (2014) Perryman, M., Hartman, J., Bakos, G. Á., & Lindegren, L. 2014, ApJ, 797, 14
  * Petigura et al. (2013) Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academy of Science, 110, 19273
  * Rasmussen & Williams (2006) Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian processes for machine learning (MIT Press)
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Rogers (2015) Rogers, L. A. 2015, ApJ, 801, 41
  * Rowe et al. (2015) Rowe, J. F., Coughlin, J. L., Antoci, V., et al. 2015, ApJS, 217, 16
  * Shvartzvald et al. (2016) Shvartzvald, Y., Maoz, D., Udalski, A., et al. 2016, MNRAS, 457, 4089
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000
  * Stumpe et al. (2012) Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985
  * Tremaine & Dong (2012) Tremaine, S., & Dong, S. 2012, AJ, 143, 94
  * Uehara et al. (2016) Uehara, S., Kawahara, H., Masuda, K., Yamada, S., & Aizawa, M. 2016, ApJ, 822, 2
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Wang et al. (2013) Wang, J., Fischer, D. A., Barclay, T., et al. 2013, ApJ, 776, 10
  * Wang et al. (2015) —. 2015, ApJ, 815, 127
  * Weiss & Marcy (2014) Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, L6
  * Winn (2010) Winn, J. N. 2010, ArXiv e-prints, arXiv:1001.2010
  * Wolfgang et al. (2016) Wolfgang, A., Rogers, L. A., & Ford, E. B. 2016, ApJ, 825, 19
  * Yee & Gaudi (2008) Yee, J. C., & Gaudi, B. S. 2008, ApJ, 688, 616