Revisiting the Long-Period Transiting Planets from Kepler

TL;DR

This study revisits long-period transiting planets from Kepler data, identifying 15 candidates, and finds that cold Jupiters are common around Sun-like stars, especially in systems with inner planets, with occurrence rates consistent with radial velocity surveys.

Contribution

It applies an improved pipeline to a larger sample with updated stellar data, confirming the prevalence of long-period giant planets and their association with systems containing inner planets.

Findings

15 planet candidates identified, including 2 new ones.

Planet occurrence rate for 2-10 year periods is about 0.7 per star.

Large outer planets are mostly found in systems with small inner planets.

Abstract

Currently, we have only limited means to probe the presence of planets at large orbital separations. Foreman-Mackey et al. searched for long-period transiting planets in the Kepler light curves using an automated pipeline. Here, we apply their pipeline, with minor modifications, to a larger sample and use updated stellar parameters from Gaia DR2. The latter boosts the stellar radii for most of the planet candidates found by FM16, invalidating a number of them as false positives. We identify 15 candidates, including two new ones. All have sizes from 0.3 to 1 $R_{\rm J}$, and all but two have periods from 2 to 10 yr. We report two main findings based on this sample. First, the planet occurrence rate for the above size and period ranges is $0.70^{+0.40}_{-0.20}$ planets per Sun-like star, with the frequency of cold Jupiters agreeing with that from radial velocity surveys. Planet occurrence rises with decreasing planet size, roughly describable as $dN/d\log R \propto R^{\alpha}$ with $\alpha = -1.6^{+1.0}_{-0.9}$, i.e., Neptune-sized planets are some four times more common than Jupiter-sized ones. Second, five out of our 15 candidates orbit stars with known transiting planets at shorter periods, including one with five inner planets. We interpret this high incidence rate to mean: (1) almost all our candidates should be genuine; (2) across a large orbital range (from $\sim 0.05$ to a few astronomical units), mutual inclinations in these systems are at most a few degrees; and (3) large outer planets exist almost exclusively in systems with small inner planets.

Tables

Table 1. Table 1 System Parameters of the Long-period Transiting Exoplanet Candidates

KIC ID	$R_{\ast }$	$K_p$	$t_0$	Period	Radius	$b$	S/N	FAP	KOI/KepleraaThe KOI number and target Kepler number, as available. KIC 8505215 (also known as KOI-99) and KIC 11342550 (also known as KOI-1421) are misidentified by the Kepler pipeline as shorter period candidates. Visual inspection shows that these putative candidates only transit once and instead have much longer periods. Here, their periods are determined by the transit cord.	$n_{\mathrm{inner}}$
	($R_{\odot }$)		(KBJD)	(days)	($R_J$)
3218908c	${1.22}_{-0.06}^{+0.05}$	$14.6$	${766.6769}_{-0.0094}^{+0.0078}$	${1246.35}_{-83.09}^{+150.22}$	${0.783}_{-0.044}^{+0.039}$	${0.00}_{-0.39}^{+0.39}$	$90.6$	$0.04$	1108/770	3
3239945d	${0.74}_{-0.03}^{+0.03}$	$14.0$	${420.2867}_{-0.0004}^{+0.0004}$	${1071.23}_{-0.01}^{+0.01}$	${0.904}_{-0.038}^{+0.038}$	${0.14}_{-0.17}^{+0.10}$	$358.6$	$0.27$	490/167	3
4754460	${1.17}_{-0.09}^{+0.09}$	$14.9$	${826.8364}_{-0.0044}^{+0.0045}$	${1772.14}_{-326.60}^{+174.75}$	${0.662}_{-0.060}^{+0.058}$	${0.88}_{-0.10}^{+0.03}$	$30.6$	$0.05$
6551440	${1.14}_{-0.08}^{+0.08}$	$13.6$	${1039.0583}_{-0.0054}^{+0.0051}$	${1018.14}_{-97.08}^{+194.53}$	${0.374}_{-0.031}^{+0.032}$	${0.76}_{-0.07}^{+0.08}$	$25.0$	$0.03$
8410697b,d	${1.05}_{-0.07}^{+0.08}$	$13.4$	${542.1321}_{-0.0018}^{+0.0017}$	${1047.83}_{-0.01}^{+0.01}$	${0.727}_{-0.050}^{+0.057}$	${0.09}_{-0.29}^{+0.28}$	$73.2$	$0.05$
8505215c	${0.8}_{-0.03}^{+0.04}$	$13.0$	${140.0488}_{-0.0019}^{+0.0018}$	${2608.45}_{-114.47}^{+195.89}$	${0.310}_{-0.013}^{+0.017}$	${0.04}_{-0.35}^{+0.32}$	$104.7$	$0.04$	99/none
8800954d	${0.83}_{-0.03}^{+0.04}$	$13.4$	${492.7665}_{-0.0010}^{+0.0009}$	${704.20}_{-0.01}^{+0.01}$	${0.411}_{-0.016}^{+0.021}$	${0.18}_{-0.26}^{+0.22}$	$121.2$	$0.05$	1274 /421
9306307	${1.38}_{-0.11}^{+0.11}$	$14.0$	${1191.3562}_{-0.0001}^{+0.0001}$	${730.83}_{-0.24}^{+0.52}$	${3.296}_{-0.266}^{+0.264}$	${0.64}_{-0.02}^{+0.01}$	$31.8$	$1.00$
10187159c	${0.78}_{-0.05}^{+0.06}$	$14.4$	${604.1067}_{-0.0025}^{+0.0024}$	${1006.63}_{-47.67}^{+95.27}$	${0.577}_{-0.039}^{+0.048}$	${0.02}_{-0.38}^{+0.35}$	$29.2$	$0.09$	1870/989	1
10602068	${0.84}_{-0.06}^{+0.07}$	$14.9$	${830.8089}_{-0.0002}^{+0.0002}$	${769.19}_{-15.98}^{+15.45}$	${2.368}_{-0.171}^{+0.199}$	${0.60}_{-0.01}^{+0.01}$	$693.6$	$1.00$
10842718b	${0.86}_{-0.06}^{+0.07}$	$14.6$	${226.2347}_{-0.0041}^{+0.0040}$	${8375.64}_{-267.18}^{+473.31}$	${0.580}_{-0.041}^{+0.048}$	${0.01}_{-0.29}^{+0.28}$	$38.2$	$0.10$
11709124c	${1.07}_{-0.05}^{+0.05}$	$14.5$	${657.2675}_{-0.0018}^{+0.0020}$	${1183.93}_{-88.71}^{+148.27}$	${0.926}_{-0.048}^{+0.050}$	${0.21}_{-0.63}^{+0.32}$	$155.7$	$0.06$	435/154	5
6186417	${1.34}_{-0.13}^{+0.15}$	$15.4$	${958.7534}_{-0.0059}^{+0.0062}$	${936.06}_{-107.32}^{+203.22}$	${0.733}_{-0.081}^{+0.093}$	${0.05}_{-0.57}^{+0.57}$	$32.6$	$0.17$
6234593	${0.98}_{-0.08}^{+0.09}$	$15.7$	${1147.9908}_{-0.0003}^{+0.0003}$	${732.02}_{-1.13}^{+2.50}$	${2.298}_{-0.188}^{+0.212}$	${0.42}_{-0.01}^{+0.01}$	$293.1$	$1.00$
7906827d	${1.17}_{-0.11}^{+0.11}$	$15.7$	${772.1934}_{-0.0019}^{+0.0020}$	${737.11}_{-0.01}^{+0.01}$	${0.955}_{-0.099}^{+0.093}$	${0.03}_{-0.29}^{+0.27}$	$106.8$	$0.02$
7947784	${1.35}_{-0.12}^{+0.12}$	$15.5$	${905.2546}_{-0.0006}^{+0.0006}$	${741.58}_{-8.02}^{+15.24}$	${2.713}_{-0.245}^{+0.243}$	${0.75}_{-0.02}^{+0.02}$	$266.1$	$0.98$
9704149b	${0.91}_{-0.06}^{+0.07}$	$15.1$	${419.7203}_{-0.0036}^{+0.0035}$	${1245.41}_{-119.54}^{+312.09}$	${0.464}_{-0.033}^{+0.039}$	${0.03}_{-0.56}^{+0.60}$	$35.8$	$0.32$
10525077b,d	${1.15}_{-0.09}^{+0.1}$	$15.4$	${335.2493}_{-0.0065}^{+0.0080}$	${854.09}_{-0.01}^{+0.01}$	${0.541}_{-0.047}^{+0.053}$	${0.47}_{-0.98}^{+0.31}$	$40.0$	$0.54$	5800/459	1
11342550c	${1.04}_{-0.06}^{+0.06}$	$15.3$	${524.2814}_{-0.0019}^{+0.0020}$	${1632.13}_{-111.35}^{+151.10}$	${0.889}_{-0.057}^{+0.058}$	${0.08}_{-0.52}^{+0.42}$	$169.5$	$0.03$	1421/none

Table 2. Table 2 Parameter Distributions for Simulated Transits

Name	Parameter	Distribution
Period	$\log P$	$𝒰(\log 2\mathrm{yr},\log 25\mathrm{yr})$
Radius Ratio	$\log R_p/R_{\star }$	$𝒰(\log 0.02,\log 0.3)$
Impact Parameter	$b$	$𝒰(0,1+R_p/R_{\star })$
Eccentricity	$e$	$\beta (0.867,3.03)$aaKipping (2013)
	$\omega$	$𝒰(0,2\pi )$
Limb Darkening	$q_1$	$𝒰(0,1)$
	$q_2$	$𝒰(0,1)$

Table 3. Table 3 Occurrence Rates of Outer Planets

$R_p\left(R_J\right)$	$K_p$ Range	Occurrence RateaaThe occurrence rate is expressed as the number of planets per Sun-like star for the given radius bin, and integrated over $2$ to $10$ yr. The reported values are taken from the peak of the posterior distribution, and the uncertainties are taken as the values at $\exp \left(-0.5\right)$ from the peaks.	$n_{\mathrm{pl}}$
$0.30-0.45$	bright	${0.52}_{-0.29}^{+0.46}$	2
	mid
	dim
	all	${0.50}_{-0.28}^{+0.44}$	$2$
$0.45-0.67$	bright	${0.07}_{-0.05}^{+0.10}$	1
	mid	${0.31}_{-0.17}^{+0.28}$	2
	dim	${0.33}_{-0.23}^{+0.45}$	1
	all	${0.17}_{-0.07}^{+0.10}$	$4$
$0.67-1.00$	bright	${0.20}_{-0.08}^{+0.12}$	4
	mid	${0.06}_{-0.04}^{+0.09}$	1
	dim	${0.13}_{-0.07}^{+0.12}$	2
	all	${0.14}_{-0.04}^{+0.06}$	$7$
$0.30-1.00$	totalbbIncluded in Wang et al. (2015).	${0.70}_{-0.20}^{+0.40}$	$13$

Equations

$$p_{\mathrm{geo}}=\frac{R_{\ast }+R_{\mathrm{p}}}{a}\frac{1+e\sin \omega }{1-e^2}$$

$$p_{\mathrm{geo}}={\left(\frac{4{\pi }^2}{G{M}_{\ast }}\right)}^{1/3}\left(\frac{1+e\sin \omega }{1-e^2}\right)\left(R_{\ast }+R_{\mathrm{p}}\right)P^{-2/3}$$

$$p_{\mathrm{win}}=\left\{\begin{array}[]{ll}1-\left(1-f_{\rm duty}\right)^{T/P}&\quad P\leq T\\ Tf_{\rm duty}/P&\quad P>T\end{array}\right.$$

$$p_{\det }=p_{\mathrm{rec}}{p}_{\mathrm{geo}}{p}_{\mathrm{win}}.$$

$$\frac{d^{2}N}{d\log Pd\log R_{\mathrm{p}}}=A{\left(\frac{R_{\mathrm{p}}}{R_{\mathrm{J}}}\right)}^{\alpha }{\left(\frac{P}{\mathrm{year}}\right)}^{\beta }.$$

$${\bar{N}}_{\mathrm{lp}}=\sum _i{N}_{\star }^{(i)}\int {\text{displaylimits}}_{2\mathrm{yrs}}^{10\mathrm{yrs}}\int {\text{displaylimits}}_{R_{p,\min }}^{R_{p,\max }}A{\left(\frac{R_{\mathrm{p}}}{R_J}\right)}^{\alpha }{p}_{\det }d\log Pd\log R_{\mathrm{p}},$$

$$P(N_{\mathrm{lp}}|n_{\mathrm{lp}})=\frac{{\bar{N}}_{\mathrm{lp}}^{N_{\mathrm{lp}}}\exp (-{\bar{N}}_{\mathrm{lp}})}{N_{\mathrm{lp}}!}.$$

$$P(n_{\mathrm{lp}}|N_{\mathrm{lp}})=\frac{P(N_{\mathrm{lp}}|n_{\mathrm{lp}})P(n_{\mathrm{lp}})}{P(N_{\mathrm{lp}})}\propto P(N_{\mathrm{lp}}|n_{\mathrm{lp}}),$$

$$\mathcal{L}=A\log (5)e^{-{\bar{N}}_{\mathrm{lp}}}\prod _i^{N_{\mathrm{lp}}}{\left(\frac{R_{\mathrm{p},i}}{R_{\mathrm{J}}}\right)}^{\alpha }{p}_{\det ,i},$$

$$\frac{N_{\mathrm{lp}}^{\mathrm{in}}}{N_{\mathrm{lp}}}=P(\mathrm{SE}|\mathrm{Outer})p_{\mathrm{geo}}^{\mathrm{in}}{p}_{\mathrm{rec}}^{\mathrm{in}}.$$

$$p_{\mathrm{geo}}\approx \frac{R_{\ast }/a_{\mathrm{in}}}{\sin i_0}.$$

$$P(\mathrm{SE}|\mathrm{Outer})=89\%\left(\frac{N_{\mathrm{lp}}^{\mathrm{in}}/N_{\mathrm{lp}}}{5/13}\right)\left(\frac{0.03}{R_{\ast }/a_{\mathrm{in}}}\right)\left(\frac{i_0}{4^{\circ }}\right).$$

$$\xi \equiv \frac{T_{\mathrm{dur},\mathrm{in}}/P_{\mathrm{in}}^{1/3}}{T_{\mathrm{dur},\mathrm{out}}/P_{\mathrm{out}}^{1/3}}.$$

$$\frac{|R_{\mathrm{p},\mathrm{left}}-R_{\mathrm{p},\mathrm{right}}|}{R_{\mathrm{p},\mathrm{original}}}\le 0.005.$$

$$S/N\propto {\left(\frac{R_{\mathrm{p}}}{R_{\ast }}\right)}^2\sqrt{\Delta T\times F}.$$

$${\left(\frac{dN}{d\log R_{\mathrm{p}}}\right)}_{\mathrm{obs}}$$

Full text

I. Introduction

Currently, we have limited means to probe planet populations at orbital separations on the order of an astronomical unit or above. These planets are difficult to identify in transit surveys like the Kepler mission as they only transit once or twice within the observational baseline. Furthermore, unless they have Jovian masses, their radial velocity (RV) signals are hard to detect. Though planet microlensing studies have yielded an interesting sample, including planets with Neptune masses, such a sample is small and is mostly composed of planets orbiting M-dwarfs (Gould et al., 2010). The elusive nature of these long-period planets is a roadblock on our path toward a complete census of planet populations and a successful theory of planet formation. It has also recently become apparent that long-period planets may be correlated with short-period ones (Bryan et al., 2019; Zhu & Wu, 2018), and may influence the dynamical evolution of the latter. Given this, it seems relevant and worthwhile to expand our knowledge of such planets. In this work, we pursue this task by searching for transiting long-period planets in stellar light curves obtained by the Kepler mission.

The Kepler mission (Borucki et al., 2010) has been responsible for the discovery of thousands of transiting exoplanets and planetary candidates (e.g., Thompson et al., 2018). Relative to other transit surveys such as the Transiting Exoplanet Survey Satellite (TESS) (Ricker et al., 2014) or even the upcoming PLATO mission (Rauer et al., 2014), Kepler’s long observational baseline makes it uniquely capable of probing planets at a wide range of orbital separations. However, as the transit probability decreases rapidly with separation ($\sim 10^{-3}$ at 5 au), many recent exoplanet population studies have restricted themselves to the distribution and occurrence rate of relatively short-period planets (e.g., Fressin et al., 2013; Petigura et al., 2013; Dressing & Charbonneau, 2015).

Long-period transit events are overlooked by standard search procedures, which require three or more transits to be observed. Yee & Gaudi (2008) first introduced the idea of searching for single-transit events in the Kepler data, and a handful of studies have since sought to increase the sample of known long-period planets through visual inspection of individual light curves (e.g., Wang et al., 2015; Uehara et al., 2016). However this method has its own shortcomings; it is difficult to determine the detection efficiency of such a search procedure, which is critical for studying the underlying planet population.

Foreman-Mackey et al. (2016) (hereafter FM16) were the first to perform an automated search for long-period transiting planets. FM16 sought long-period, high signal-to-noise ratio (S/N) transit signals in the light curves of 39,036 bright main sequence GK stars. They identified 16 long-period planet candidates (two of which they label as likely false positives). By injecting artificial transits into the light curves, they obtained the search completeness and detection efficiency for their pipeline, allowing them to infer the occurrence rate of planets in the outer region. FM16 reported a rate of $0.42\pm 0.16$ planets per star within a radius range of $0.4~{}R_{\rm J}<R_{\rm p}<1~{}R_{\rm J}$ and a period range of 2 yr $<P<~{}$25 yr. Their work constitutes significant progress toward understanding the population of outer planets.

However, a major change has occurred since FM16 that warrants a new study. FM16 adopted stellar parameters from the Kepler Input Catalog (KIC) (Brown et al., 2011; Huber et al., 2014), but the recent second data release of the Gaia mission (Gaia Collaboration et al., 2016, 2018) has updated the stellar radii for a significant fraction of the stars in FM16’s target catalog (e.g., Berger et al., 2018). In particular, many of the stars identified as GK dwarfs in the Kepler catalog have larger radii than initially thought, meaning that the sizes of many of FM16’s planet candidates may have been underestimated. A number of candidates are now too large to be compatible with sub-stellar objects, and many small planets are now Jovian in size. This upgrade necessarily requires FM16’s planet occurrence rate to be adjusted.

FM16 adopted the philosophy of open-source software and have made their entire pipeline publicly available. This affords us great ease in our work. Here, we repeat their exercise but with a number of modifications. First, we apply the search procedure of FM16 to more than 61,000 Kepler light curves to identify long-period planet candidates, using the revised stellar properties of Berger et al. (2018) to determine their transit characteristics. We also identified a missing factor of $\pi^{1/3}$ in the transit probability calculation in FM16’s code. This alone reduces FM16’s occurrence rate by a factor of $1.46$. We adopt a procedure different from that in the FM16 code to convert transit depth to planet radius, after a discrepancy in our earlier draft was identified (K. Masuda 2019, private communication). We further employ a procedure to flatten stellar variability to increase the transit S/N. We also provide an appended explanation on why so few double-transit systems are found relative to single-transit ones, a puzzle posed but not well-explored by FM16.

After performing injection and recovery tests to compute the completeness of the search procedure, we use our new sample of planet candidates to provide a new estimate for the planet occurrence rate at long periods. We compare this to earlier results from RV and microlensing studies. Moreover, we carry out a number of detailed analyses to clarify the significance of our results. One intriguing finding of FM16 is that a large fraction of transiting outer planets also have inner transiting companions, including one with five planets. Given the rarity of transiting systems observed by Kepler, this is not expected unless the inner and outer planetary systems are correlated in occurrence rate and are inclined similarly. We investigate this behavior in detail. Our sample also affords us an opportunity to study the size distribution of outer planets, which provides important input for theories of planet formation.

II. Searching for Candidates

II.1. Target Selection

Our candidate search method makes use of the open-source code peerless111https://github.com/dfm/peerless from FM16, closely following their procedure with a few modifications. We refer the reader to FM16 for a detailed description of their original search and vetting procedure, and we discuss the changes we implement in Appendix A. Below, we describe the sample of light curves to which this modified procedure is applied.

Our initial stellar sample of quiet Sun-like stars is selected so as to match that of FM16, with an additional 50,817 stars dimmer than their magnitude limit of 15 (up to $K_{p}=17$ mag). Our initial selection criteria are:

• •

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

• •

$R_{*}\leq 1.15~{}R_{\odot}$,

• •

data span $\geq$ 2 yr,

• •

duty cycle $\geq$ 0.6,

• •

$K_{p}\leq 17$ mag, and

• •

CDPP${}_{7.5\rm hr}\leq$ 1000 ppm.

However, we further restrict our sample based on the updated stellar radii provided by Berger et al. (2018), who combine precise parallax measurements from Gaia DR2 (Gaia Collaboration et al., 2016, 2018) with the stellar parameters of the Kepler DR25 catalog to produce revised stellar radii for a large fraction of stars in the Kepler field. Berger et al. (2018) claim that their parallax-derived radii have an uncertainty of $\sim$8%, 4–5 times more precise than those provided by the latest Kepler catalog of stellar properties. For stars with poorly constrained parallax/distance measurements or missing photometry, no updated radius is reported. From their catalog, we find revised radii for $\sim$95% of our targets, and we discard the remaining targets for which no revised radius is available. For a comparison of all stellar radii reported by Kepler and Berger et al. (2018) in our target sample, see Figure 1.

For a considerable fraction of our target sample, the new radii push the stars to sizes inconsistent with a Sun-like sample. Keeping these stars in our sample would make it difficult to constrain the occurrence rate of long-period transiting planets without taking the differences in stellar type into account; for instance, many of the assumed GK dwarfs in our sample should now be considered giants based on their revised radii. The properties of such stars may have an effect on planet occurrence, and we therefore choose to limit our sample to targets with $0.7~{}R_{\odot}\leq R_{*}\leq 1.5~{}R_{\odot}$ while keeping all other parameters the same. This reduces our sample by $\sim$ 28%, resulting in a final selection of 61,418 Sun-like stars. While this is by no means the most straightforward method of cutting our stellar sample, the general conclusions we draw should have little dependence on the selection method, given the large number of stars we consider. A histogram displaying the magnitude distribution of our sample is shown in Figure 2.

Following our target selection, we apply our modified search and vetting procedure to the pre-search data conditioning (PDC) light curves. A complete description of this process can be found in FM16. This yields $19$ transit candidates, out of $61,418$ targets. To constrain the transit parameters of each candidate, we adopt an MCMC fit (described below), rather than the output from peerless. 222This revision was prompted by Kento Masuda, who notified us that the planet radii in our earlier draft, using peerless, did not match those expected based on the transit depths and stellar radii. Surprisingly, results in FM16 did not share this problem. The cause for this discrepancy is unclear. Finally, we use peerless to determine the false-positive probability for each event.

II.2. MCMC Fit to the Transit Light Curves

The physical parameters of these 19 candidates are determined by performing an MCMC fit to the light curves. We use the batman package (Kreidberg, 2015) to generate the transit light curve for each candidate in isolation, assuming a quadratic limb darkening law and circular orbit, and use the emcee package (Foreman-Mackey et al., 2013) to constrain the transit parameters.

We first detrend the light curve, using a third-order spline, after masking the transit based on its duration, which is estimated in the peerless search output. The flux errors are determined from the standard deviation of the out-of-transit light curve. The fitting parameters are the time of mid-transit ($t_{0}$), logarithmic period ($\log(P)$), planet-star radius ratio ($R_{\rm p}/R_{*}$), impact parameter ($b$), and the quadratic limb darkening coefficients ($u_{1}$ and $u_{2}$). The likelihood for the MCMC sampling $\mathcal{L}$ is calculated as $\mathcal{L}\propto\exp(-\chi^{2}/2)$, where $\chi^{2}$ is the sum of the standard chi-squared for each transit.

For candidates with two transits, we fit each transit separately, then set the period as the difference between their respective $t_{0}$ values. We then phase-fold the transits on this period and fit the resulting phase-folded transit. For candidates with only one transit, we first use the isocrones package (Morton, 2015) to determine the mean stellar density and its uncertainty from the Gaia DR2 parallax (Gaia Collaboration et al., 2018), corrected for zero-point offset with 0.03 mas, and the KIC broadband photometry. We then use these values as a prior in our transit fits to constrain the planet’s orbital period (see Seager & Mallén-Ornelas, 2003).333This fails for eclipsing binaries where both components contribute to the gravitational mass. This explains the periods of the largest candidates in our catalog, which we are not concerned with as they are likely eclipsing binaries.We also place a constraint on the period such that no additional transits can occur within the Kepler baseline. As stated earlier, we assume a circular orbit when fitting each transit. Introducing a non-zero eccentricity only increases the uncertainty on the inferred orbital period while having little influence on the median (for an illustration, see Section 3.3 of Villanueva et al., 2019).

II.3. Results of the Candidate Search

The results of our MCMC fits to the candidate transits are given in Table 1. Twelve of these candidates are common to FM16, and we provide new constraints on their transit parameters based on the revised stellar radii of Berger et al. (2018). We note that, by restricting our target list based on the revised stellar radii, we lose four of the candidates identified in FM16 (KIC 8426957, KIC 8738735, KIC 10287723, and KIC 10321319) as their stellar parameters are no longer consistent with those of Sun-like stars. The remaining seven transit light curves not common to FM16’s sample are shown in Figure 3.

We also consider the overlap between our sample and that of both Wang et al. (2015) and Uehara et al. (2016). Our list of target light curves contains ten candidates identified by Wang et al. (2015), four of which we recover (see Table 1). The remainder either do not meet our S/N threshold (KIC 5536555, KIC 9413313, KIC 9662267, KIC 11716643, and KIC 12454613) or are rejected for close proximity to a data gap (KIC 6191521). Seven of the candidates found in Uehara et al. (2016) are also in our target sample; of these, we recover five candidates. The other two are rejected for having a grazing impact parameter (KIC 3230491) and for close proximity to a data gap (KIC 6191521, also found in Wang et al., 2015). We note that these exclusions are automatically included in our completeness calculation and so do not affect our determination of the occurrence rate described below. The radii and orbital periods of all planet candidates are plotted in Figure 4.

Of the seven candidates displayed in Figure 3, two have two observed transits within the baseline of the Kepler observations (KIC 7906827 and KIC 10525077), the latter of which was previously discovered by Wang et al. (2015) (Kepler-459b). Two additional single-transit candidates, KIC 9704149 and KIC 11342550, were reported in Wang et al. (2015) and Uehara et al. (2016) (KOI-1421), respectively. Thus, four of the long-period transit candidates in our sample have never been published before. However, based on our false-alarm probability (FAP) estimates (see Table 1), it is very likely that two of these candidates are not transiting planets but are instead eclipsing binary stars. In Section II.5, we further discuss potential astrophysical false positives among our candidates.

Interestingly, five of the total 19 transit candidates have known inner companions. This means that over a quarter of our long-period candidates are found in systems with at least one shorter-period planet. Considering that only 1745 of the targets in our sample of 61,418 stars contain previously discovered planets and/or planet candidates, this abundance of long-period planets in multiple-planet systems is intriguing. Could these be special systems for which transits are more detectable? We find this to be unlikely. Photometric uncertainties in the host-star light curves are similar to those without inner systems. The S/N ratio of these candidates are also comparable to the others. In Section III, we discuss the implications of this abundance of inner planets in terms of the mutual inclination of bodies in multiple-planet systems and the architecture of such systems.

II.4. Search Completeness and Detection Efficiency

To constrain the transit detection efficiency for our target sample, we first characterize the completeness of the search procedure. We follow the commonly used method of injecting artificial transits with known parameters into the Kepler light curves, and then measure the recovery rate of the procedure (e.g., Petigura et al., 2013; Christiansen et al., 2015; Dressing & Charbonneau, 2015; Foreman-Mackey et al., 2016, and many others). To achieve this, we first randomly choose a star from the target sample and load its PDC light curve and stellar properties. We then sample planet properties from the distributions given in Table 2, and the transit is calculated and multiplied into the light curve. The search procedure is run on this transit-injected light curve, and if at least one transit within one transit duration passes all steps of the vetting process, the candidate is considered recovered. Note that, depending on the injected period and time of mid-transit, a candidate can display up to two transits in a single light curve. A more detailed analysis of the recovery rate for double- versus single transit events is described in Appendix B.

The completeness or the fraction of recovered injections as a function of radius and orbital period for our stellar sample is shown in Figure 5. We report this completeness in three separate magnitude bins, each of which contains roughly one-third of our full target sample ($6\leq K_{p}<14.8$; $14.8\leq K_{p}<15.4$; $15.4\leq K_{p}\leq 17$). The magnitude dependence of the completeness is most obvious at small fractional planet radii, where the percentage of recovered injections is essentially zero for the dimmest stars. Note that the radius is reported in terms of radius ratio between the simulated planet and its host star, $R_{\rm p}/R_{*}$. We deem this more illuminating than reporting the radius in terms of $R_{\rm J}$, as our stellar sample spans a range of radii; a $1~{}R_{\rm J}$ planet around a $0.75~{}R_{\odot}$ star would be easier to recover, for instance, compared to the same planet around a $1.5~{}R_{\odot}$ star if all other factors are kept the same. The decrease in the completeness toward the largest $R_{\rm p}/R_{*}$ values is due in part to our choice of injected planetary properties, chosen to maximize the recovery of signals with smaller radii (FM16). For instance, signals with $R_{\rm p}/R_{*}$ values above $\sim 0.1$ are more likely to be rejected during the vetting process of the search pipeline for having a grazing impact parameter. This parameter is influenced by the radius both in its injected distribution ($b\leq 1+R_{\rm p}/R_{*}$) (Table 2) and acceptance criteria ($b\leq 1-R_{\rm p}/R_{*}$) (FM16).

This completeness provides an estimate of the probability of detecting a transit with some set of planetary parameters, given that the planet transits its host star within the observing period. We denote this probability as $p_{\mathrm{rec}}$. The total detection efficiency is then determined by combining this with the geometric transit probability and the window function – in other words, the probability that a planet will transit along our line of sight given its physical parameters, and the probability that it transits within the baseline of our observations. The geometric transit probability for a given system is expressed as

[TABLE]

(Winn, 2010) where $R_{*}$ is the stellar radius, $R_{\rm p}$ is the planet radius, $a$ is the semi-major axis, $e$ is the orbital eccentricity, and $\omega$ is the argument of periastron. This equation can be rewritten using Kepler’s third law such that

[TABLE]

where $M_{*}$ is the stellar mass and $P$ is the orbital period of the planet.

Like Burke & McCullough (2014) and FM16, we use the binomial probability of observing a minimum of one transit to approximate the window function

[TABLE]

where $f_{\rm duty}$ is the duty cycle and $T$ is the baseline of observations for the star being considered. The total detection efficiency is then

[TABLE]

Because our initial completeness is calculated for three separate ranges in magnitude, this detection efficiency is also magnitude-dependent. We take this into account when computing the planet occurrence rate in the following section.

II.5. False Positives

Before we proceed to analyze our results, we discuss the possibility of false positives within our catalog of 19 planet candidates.

The general increase in stellar radii following the second data release of Gaia (see Figure 1) plays an important role when identifying false positives in our sample. The resulting jump in the size of our planet candidates necessarily changes their physical association; while some candidates simply move from Neptune-sized into the Jovian range, others exceed our expectations for maximum planet radii and thus are unlikely to be planets. Some of the candidates identified by FM16 belong to the latter class.

Among our catalog of 19 candidates (Table 1), the four largest ones have radii larger than $2~{}R_{\rm J}$. The peerless package assigns them a low planet probability (FAP $\geq 0.9$). They are more likely eclipsing binaries blended with foreground stars or with physically associated stellar companions (Foreman-Mackey et al., 2016; Morton et al., 2016). We reject these four from the outset. Two candidates, with respective periods of $704$ days and $8375$ days, fall outside the period range ($2-10$ yr) we are considering. In the following, we consider the remaining 13 candidates.

Five of the remaining $13$ have known inner transiting companions. As the inner and the outer companions should be mutually inclined and do not necessarily transit simultaneously, this surprising fact itself suggests that most of these 13 candidates should be real planets (we discuss this in more detail below).

Four of our candidates are larger than Saturn. Santerne et al. (2012) reported that such large bodies have higher false-positive rates than do smaller planets, if their periods are within a couple hundred days. Currently there is no study on whether such a trend also extends to colder planets. Two of these four giants have inner transiting companions (see Figures 4 and 6), boosting their credibility as genuine planets (Lissauer et al., 2012), but are the other two likely real? To make an educated guess, we review known sources of astrophysical confusion and their prevalence:

• •

Brown Dwarfs. Because their sizes are comparable to that of giant planets, brown dwarfs may initially appear to be an important source of confusion. However, brown dwarfs are exceedingly uncommon companions for solar-type stars, particularly at periods $<2000$ days, with less than $1\%$ of companions being brown dwarfs in mass, far below the rate of planetary companions (Grether & Lineweaver, 2006). This phenomenon is commonly referred to as the brown dwarf desert (Murphy et al., 2018; Grether & Lineweaver, 2006, and references therein).

• •

M-dwarfs. With masses above $0.1~{}M_{\odot}$, M-dwarfs can have radii as low as that of Jupiter. By analyzing Figures 11 and 16 of Raghavan et al. (2010), we find a binary fraction of just $0.6\%$ for companions with size below $0.15~{}R_{\odot}=1.5R_{\rm J}$ and with period between 2 and 10 yr.

So together, brown dwarf or low-mass M-dwarf companions may masquerade as Jovian planets, but with a total occurrence rate well below $2\%$. This is much lower than the occurrence rate of cold Jupiters obtained by RV studies ($\sim 10\%$). We therefore argue that at most one, if any, of the large candidates could be a false-positive.

We have little information regarding the false positive rate for the smaller candidates in our sample. In the following analysis, we assume that they are all genuine planets. This is supported by their low FAP values, and the fact that three out of nine have inner systems.

II.6. Possible Flux Dilution

We have been prompted by the referee to consider the issue of flux dilution, which has been shown in the past to interfere with the Kepler planet sample. Many of our planet candidates are around faint stars, and multiple stars of comparable magnitudes may cohabit the same Kepler CCD pixel ($\sim 4\arcsec$ across). These stars can be either physically associated or unassociated with the host star, and their light can dilute the signal of a deeper transit into that of a shallower, planet-like, transit. Cohabiting stars are not rare events – Furlan et al. (2017) found that $\sim 30$ % of Kepler planet host stars have at least one neighbor lying within $4\arcsec$.

Querying ExoFOP444https://exofop.ipac.caltech.edu/ revealed that eight of our $15$ candidates have AO follow-up imaging (Law et al., 2014; Hirsch et al., 2017; Furlan et al., 2017; Ziegler et al., 2018), and among these, three have stellar companions with distances from $1\arcsec$ to $4\arcsec$. However, among the $13$ candidates on which we base our analysis, only KIC 3239945 has a known companion resolved by AO, with a $K$-magnitude difference of $3.6$ between the two stars and a separation of $2.2\arcsec$. Fortunately, a thorough vetting of this system, which also hosts three inner planets, has been performed by Kipping et al. (2016) and the planetary nature of the outer planet (Kepler-167e) is firmly established.

We also query Gaia DR2 for companions to the 13 stars in our final candidate sample, using their corresponding Gaia DR2 IDs (Table 1 of Berger et al., 2018). Gaia DR2 has an effective angular resolution of $\sim 0.4\arcsec$ (Gaia Collaboration et al., 2018), and can completely resolve any stellar companion beyond $1\arcsec$ (Ziegler et al., 2018). We find companions within $4\arcsec$ of three of our targets: KIC 3239945 (also listed in ExoFOP, see above), KIC 4754460, and KIC 8505215. In all cases, the point sources have Gaia $G$ magnitude differences greater than 5. The large magnitude differences make it unlikely that our candidates are stellar transits diluted by light from a companion. If the transit signal comes from the brighter source, its depth will be little affected by the companion light; if instead the transit signal arises from the dimmer object, the transit would have to be intrinsically deep ($47\%$ depth for KIC 4754460 and $19\%$ for KIC 8505215) and V-shaped, accompanied by secondary transits. None of these features are compatible with observations.

What about star pairs that cannot be resolved by AO and Gaia? This is increasingly rare for those that are not physically associated, but more likely for physically bound ones. Given our sample, the only possible scenario for false-positives is that in which a solar-type star is orbited by two M-dwarfs that transit each other. This would, however, create a very different transit shape, which should be identified by the peerless vetting procedure and FAP analysis. As such, we consider this unlikely.

III. Discussion

Here, we discuss the implications of our work on three issues: the occurrence rate of long-period planets, the size distribution of such planets, and the inner-outer planet correlation.

III.1. Occurrence Rate for Outer Planets

Having justified that the majority of our planet candidates are likely true planets, we proceed to constrain the planet occurrence rate. We perform this calculation, using the detection efficiency as determined in Section II.4, for planets from 0.3 to 1 $R_{\rm J}$.

As alluded to earlier, the occurrence estimate by FM16 has to be revised using the updated stellar radii from Berger et al. (2018). In fact, many of their candidates are now shifted to sizes larger than Jupiter. Furthermore, while the equation describing the geometric transit probability given in FM16 is correct, its implementation in peerless misses a factor of $\pi^{1/3}$. This decreases their reported occurrence rate by a factor of 1.46. We report a new occurrence rate, taking into account the revised radii, this minor calculation error, and the seven additional transit candidates found in our sample.

To estimate this occurrence rate, we assume that the uncertainties on the planet properties, such as radius and orbital period, are negligible. This is a permissible simplification because the mean occurrence rate is calculated in bins considerably larger than the uncertainties. Additionally, like FM16 we assume that no other transits occur in a data gap for those candidates with a single-transit event, i.e.,the orbital period we determine is assumed correct. Each candidate is also assumed to orbit a star accurately described by the Kepler catalog of Huber et al. (2014), with updated radii from Berger et al. (2018), rather than a companion or background star. We also exclude the likely false positives contained in our sample when computing the occurrence rate, and make the simplifying assumption that all of our planet candidates have zero eccentricity; the latter has a negligible effect on the outcome of our results. We further assume that the distribution of planets in the parameter space of interest follows the relation

[TABLE]

In the following, we set $\beta=0$. This assumed period insensitivity is justified both because we are in a relatively narrow range ($2-10$ yr), and because RV studies have not found any strong period dependence for giant planets in this range (Cumming et al., 2008; Bryan et al., 2016). The total expected number of detected planets, for $N^{(i)}_{*}$ stars in magnitude bin $i$, is

[TABLE]

where the detection efficiency $p_{\rm det}$ is given by equation (4).

We divide the radius range covered by our candidates ($0.3-1.0~{}R_{\rm J}$) into three equal logarithmic bins. Within each narrow radius bin, we can adopt $\alpha=0$, but we let $A$ vary across the bins.

The probability of getting the actual number of planets we detect, $N_{\rm lp}$, given its expectation value $\bar{N}_{\rm lp}$, is

[TABLE]

According to Bayes’ theorem, the posterior distribution of $n_{\rm lp}$ is given by

[TABLE]

because both $P(n_{\rm lp})$ and $P(N_{\rm lp})$ are essentially constant. The posterior probability distributions for the three radius bins are presented in Figure 7. The separation by magnitude bins is a precaution observed because the detection completeness is dependent on the stellar magnitude. In practice, we find such a precaution unnecessary – it makes little difference to the results whether we separate the computation by magnitude bin and combine the results to determine the total occurrence rate (black curves in Figure 7), or use the search completeness computed for all magnitudes (the bottom right panel of Figure 5) to determine the total rate.

The derived occurrence rates for different radius bins and separate magnitude bins are presented in Table 3. The center values are taken to be the peak of the posterior distribution, and the $1-\sigma$ upper and lower bounds are taken to be where the probabilities have fallen by a factor $\exp{(-0.5)}$ from the peak.555Although the standard procedure is to report the 68% range centered on the median as the 1-$\sigma$ values, we would like to be consistent with our later approach of using maximum likelihood to constrain the values of $A$ and $\alpha$. Our results show that planet occurrence rises with decreasing planet sizes, from $0.14$ for the largest size bin (Jupiter-sized) to $0.50$ for the lowest bin (Neptune-sized). There are more small planets than large giants – the number of detections are smaller for the smaller planets, but their detection efficiency drops even faster, demanding a larger underlying occurrence rate. This result is not sensitive to our choice of bin size – the trend persists when we increase the number of bins to four or five.

To account for this size dependency in our calculation of the total occurrence rate, we need to introduce a nonzero $\alpha$ into our analysis. This is detailed in Section III.2. The final result is $0.70^{+0.4}_{-0.2}$ planets per Sun-like star, with sizes between 0.3 and 1 $R_{\rm J}$, and periods between 2 and 10 yr.

Here, we attempt to compare our results to those from other studies. Long-term RV surveys have monitored a large number of stars and have obtained occurrence rates for cold giant planets. Mayor et al. (2011) concluded a planet occurrence rate of $0.139\pm 0.017$ for planet mass $m\sin i>50~{}M_{\oplus}$ and period $P<$ 10 yr. A more convenient comparison is against that of Cumming et al. (2008), where they expressed the planet occurrence rate as a function of planet mass and period. Evaluating their parameterized form for the mass range $M_{p}\geq M_{\rm Saturn}$ and the period range $2-10$ yr, we obtain an occurrence rate of $0.05-0.1$ planets per FGK star. To compare this against our results, we assume that planets with masses $M_{p}\geq M_{\rm Saturn}$ are also larger than Saturn. For the size range $R_{\rm Saturn}-R_{\rm J}$, we find an occurrence rate of $0.07_{-0.03}^{+0.04}$ planets per Sun-like star, in good agreement with Cumming et al. (2008).

III.2. Size Distribution of Outer Planets

A unique outcome of our study is the size distribution for outer planets, something hitherto unexplored. In the following, we investigate the observed size distribution and discuss its implications.

The likelihood of producing the observed radius distribution from the underlying distribution (Equation 5) is given by Poisson statistics (e.g., Gould et al., 2010; Suzuki et al., 2016):

[TABLE]

where $R_{{\rm p},i}$ is the radius of the $i^{\rm th}$ planet and $p_{{\rm det},i}$ its detection efficiency (a function of radius and period).

We evaluate the likelihood values over a range of $\alpha$ and $A$ and present the results in Figure 8. Both parameters are constrained reasonably well. In particular, the power-law slope $\alpha$ is determined to be $\alpha=-1.6_{-0.9}^{+1.0}$, meaning the occurrence rate rises toward smaller planets. Figure 9 shows the best-fit size distribution for our sample, as well as the 1-$\sigma$ allowed range. The occurrence rates derived from the three radius bins separately (assuming $\alpha=0$ in each bin) are also overplotted. These two different approaches obtain consistent results.

Using this derived $\alpha$ value, we can obtain an overall occurrence rate of $0.70^{+0.40}_{-0.20}$, integrated over the size range and period range we defined previously. This is consistent with, though somewhat lower than, what one would obtain by a simple summation of the three radius bins in Table 3.

How does our result compare to those from microlensing studies, which probe the mass function (but not radius distribution) of outer planets? Planet microlensing is sensitive to planets with a planet-star mass ratio down to Neptune (at $5\times 10^{-5}$ that of the Sun), orbiting at a distance of a few AUs. Most of the planet hosts in these studies are M-dwarfs (median mass $0.6~{}M_{\odot}$), due to the latter’s prevalence in the Galaxy. Suzuki et al. (2016) analyzed 30 planetary microlensing events and concluded that the planet occurrence is $0.43$ planets per star, with mass ratio $\geq 3\times 10^{-5}$ and the projected planet-star separation within $[0.45,2.22]$ Einstein radii (see also Gould et al., 2010; Cassan et al., 2012). These separations are similar in physical dimensions to those of our planets, but span a larger dynamic range (a factor of $2.35$ wider). For a crude comparison, we scale down their occurrence rate accordingly to obtain $0.18$. This is some four times lower than our value ($0.70^{+0.40}_{-0.20}$), suggesting that M-dwarfs may have fewer outer planets than do solar-type stars. On the other hand, the two occurrence rates are marginally consistent with each other, given the large error bars from both surveys and the uncertainties involved in mass-size conversion. More data is required.

The size distribution we obtain in this study is qualitatively similar to the mass distribution inferred by Suzuki et al. (2016), both indicating that Neptune-massed planets are a few times more abundant than Jovian-massed ones.

What does our size distribution imply for theories of planet structure and formation? First, it is interesting to notice the absence of candidates between $1$ and $2R_{\rm J}$ in size, despite the high detection completeness for bodies at these sizes. This is what one would expect from calculations of planetary structures: independent of mass and composition, all cold planets should have sizes below $1R_{\rm J}$ (e.g., Fortney et al., 2007).

Second, while the rising size function toward smaller planets is not surprising, the fact that it is continuous across the size range is unexpected. Planets with sizes between Neptune and Jupiter (our middle size bin, $0.45-0.67R_{\rm J}$) are not present in our solar system. Lacking appropriate analogs, one may call them super-Neptunes or sub-Saturns. We find that they are at least as common as Jovian-sized ones, if not more so.

In contrast, population synthesis models based on core-accretion theory (see review by Benz et al., 2014, and references therein), predict a deficit of such planets. If these planets have core-masses of, say, $25~{}M_{\oplus}$, their sizes would require hydrogen envelopes with masses from $5$ to $30~{}M_{\oplus}$ (see Figure 8 of Fortney et al., 2007). These latter values straddle the range over which runaway gas accretion is predicted to occur, and models predict that the envelopes will grow rapidly (both in mass and size), transforming the planets from Neptunes to Jupiters. As such, the abundance of these intermediate planets is puzzling.

Interestingly, a similar contradiction is noticed by Suzuki et al. (2018). They argued that the smooth planet mass function observed by microlensing surveys is incompatible with, and challenges, our current theory of core-accretion. Our results here provide another piece of evidence in this debate.

III.3. Correlation between Inner and Outer Systems

Among our list of planet candidates, five out of 13 have inner transiting companions (Figure 6). This high incidence is useful both for ascertaining the planetary nature of our candidates (Section II.5) and for studying the correlation of (a) orbital orientation and (b) occurrence between the inner and outer planetary systems. We discuss the latter point here.

The fraction of detected long-period planets with inner transiting companions depends only on the detectability of the inner planets, when the outer one is found to be transiting. Mathematically, this is

[TABLE]

Here, $P({\rm SE}|{\rm Outer})$ quantifies the correlation between the inner and outer systems (i.e.,the probability of the inner system being present when the star is known to have an outer planet), and $p_{\rm geo}^{\rm in}$ and $p_{\rm rec}^{\rm in}$ are the geometric transit probability and the recovery rate for the inner planets, respectively. In this case, we are only interested in those planets detectable by Kepler, for which reliable statistical results are available. We therefore limit ourselves to inner ones with orbital periods shorter than 400 days and radii $\gtrsim 1~{}R_{\oplus}$ (i.e., super-Earths). The planet search pipeline developed by the Kepler team (Bryson et al., 2013) has a recovery rate of nearly 100% for such planets,  so we simply adopt $p_{\rm rec}^{\rm in}=1$.

Given the presence of a known outer planet, the geometric transit probability for an inner planet will depend on the inclination dispersion of the system $i_{0}$ such that

[TABLE]

where $a_{\rm in}$ is the orbital separation of the inner one. Below, we adopt $a_{\rm in}\sim 30R_{*}$, as is representative for the inner planets here. With this and $p_{\rm rec}^{\rm in}=1$, we may now express equation (10) as

[TABLE]

Before discussing the implications of the above results, we first consider evidences for the value of $i_{0}$.

• •

In our solar system, the four terrestrial planets are inclined by $\sim 2-6^{\circ}$ from the invariable plane, i.e., the plane that is mostly determined by the orbits of Jupiter and Saturn.

• •

From analyzing Kepler statistics for the close-in Kepler planets, Zhu et al. (2018) concluded that $i_{0}\sim 0.8^{\circ}(k/5)^{\beta}$, where $k$ is the intrinsic multiplicity of the (inner) system and the power-law index $\beta\in[-4,-2]$. They also found that typically $k\sim 3$ within $400$ days. This yields an average inclination of $i_{0}\sim 4^{\circ}$.

• •

A direct constraint on the mutual inclination $i_{0}$ can be derived from the ratio of transit durations $T_{\rm dur}$ of planet pairs weighted by their respective orbital velocities ($\propto P^{-1/3}$; Steffen et al. 2010)

[TABLE]

Taking the nominal periods for the outer planets,666We use the orbital periods derived from our transit fits assuming zero eccentricity (Table 1). This is justified, given that Fabrycky et al. (2014) have shown that the parameter $\xi$ is insensitive to orbital eccentricities. we plot the resulting distribution of $\xi$ for the $5$ systems in Figure 10, where the $\pm 1\sigma$ uncertainties of individual $\xi$ measurements are also shown. We also compute the $\xi$ distributions for models with three different mutual inclinations ($1^{\circ}$, $3^{\circ}$, and $5^{\circ}$), following the method of Fabrycky et al. (2014) (see also Zhu et al. 2018). As Figure 10 shows, models with larger mutual inclinations better resemble the data. Even though we cannot distinguish between the models with $i_{0}=3^{\circ}$ and $5^{\circ}$, this test is robust in excluding the small mutual inclination case ($i_{0}=1^{\circ}$). Therefore, it is not unreasonable to use $i_{0}=4^{\circ}$ for the mutual inclination between the inner and outer planetary systems.

We conclude from the above that $i_{0}=4^{\circ}$ may be a reasonable value to describe the inner-outer mutual inclinations777This is surprising from the perspective of planet eccentricities from radial velocity surveys. Typical cold Jupiters are found to have an eccentricity dispersion $e_{0}\sim 0.3$. If we assume that $i_{0}\sim e_{0}/2$, as is often found in dynamical experiments, this would imply a large $i_{0}\sim 10^{\circ}$ among cold Jupiters.. Any larger value of $i_{0}$ would further strengthen our argument below.

Given the above argument, equation (12) suggests that there must be a strong correlation between small inner planets and outer planets. Almost all systems with outer planets should also have small inner planets.888 Whether the reverse is true, i.e., that all small inner planets have large outer planets, has not yet been answered. Interestingly, such a conclusion was also obtained by Zhu & Wu (2018) through analysis of RV data. They found that almost all cold Jupiters should have inner Kepler systems. Such a conclusion, now obtained independently from two studies (and extended to planets smaller than Jupiters) suggests that the formation of cold planets is conditional upon the presence of inner planets.

Turning the above argument around, we may also argue that the mutual inclinations between inner and outer planets, spanning a period range from a few days to a few years, must be low. Even assuming $P({\rm SE}|{\rm CJ})=100\%$, the mutual inclination cannot be much greater than $4^{\circ}$. These are very flat dynamical systems.

Anomalously, one of the cold giant candidates has five inner transiting planets. At face value, this is statistically unlikely given that systems with five inner transiting planets only account for $0.5\%$ (3/589) of all Kepler systems. In other words, one should expect to have 0.025 of such systems among a sample of five inner systems. This is, however, naturally explained by the above cited result of Zhu et al. (2018): because the inclination dispersion of the inner system correlates with its multiplicity, adopting $i_{0}=0.8^{\circ}$ for $k=5$ boosts the possibility of detecting a single system containing five transiting planets to 0.6.

IV. Conclusions

By employing the automated transit search pipeline of FM16 to an extended sample of target stars, we detect 19 planet candidates, including 12 in common with FM16. Updated stellar radii from Gaia DR2 allow us to exclude a number of them as false positives, and we are left with 15 likely candidates. We further restrict ourselves to $13$ of those for our analysis. Among these, five have known inner transiting companions.

The FM16 pipeline also yields the detection efficiency and search completeness. Using these, we are able to report a total occurrence rate of $0.70_{-0.20}^{+0.40}$ per Sun-like star for outer planets within a radius range of $0.3-1R_{\rm J}$ and a period range of $2-10$ yr. We also find that the radius distribution of cold planets can be expressed as a single power-law of the form $dN/d\log{R_{\rm p}}\propto R_{\rm p}^{\alpha}$, where $\alpha=-1.6^{+1.0}_{-0.9}$.

Our occurrence rate for Jovian planets is compatible with those found by long-term RV monitoring, and our finding that Neptunes are more common than Jupiters is an independent confirmation of the results from microlensing studies. While microlensing has discovered a nearly continuous mass function going from Jupiter- to Neptune-mass planets, we find a continuous size distribution between the two ends, challenging current theories of core accretion.

Finally, we investigate the seemingly surprising fact that five of our 13 long-period planets orbit stars with known transiting inner companions. This brings us to two conclusions. First, there is a strong correlation between the presence of cold planets, including both Jovian and smaller planets, with the inner system. The formation of cold planets appears to be conditioned upon the presence of small inner planets. This independently confirms the conclusions drawn by Zhu & Wu (2018) and Bryan et al. (2019) using RV data. It also provides a valuable clue for planet formation, one that unfortunately remains undeciphered. Second, the mutual inclinations between the inner and outer parts of these planetary systems must be quite small. We find that, across two decades in period, the average inclination should remain below $4^{\circ}$.

Despite our small sample size, our ability to constrain the prevalence of outer planets and elucidate the relationship between inner and outer planets is a testament to the quality and versatility of the Kepler data. Looking forward, it is important to note that current and upcoming transit surveys (e.g., TESS and PLATO) will have considerably shorter observational baselines. However, the number of planets with periods longer than these baselines will still undoubtedly be numerous (Villanueva et al., 2019). Such detections will enable further statistical analysis of the population of long-period planets, better illuminating the architecture of and relationship between planets occupying the inner and outer reaches of their systems.

The authors would like to thank Dan Foreman-Mackey and his coauthors to FM16 for the peerless code, which is well-documented and entirely open-source. The authors also thank Kento Masuda for informing us of an earlier mistake in our draft, and the anonymous referee for improving the clarity of our manuscript. This work has made use of data collected by the Kepler mission and retrieved 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. Funding for the Kepler mission is provided by the NASA Science Mission directorate. This research has also made use of the NASA Exoplanet Archive and the Exoplanet Follow-up Observation Program website, which are operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. M.K.H. and Y.W. additionally acknowledge funding from the Natural Sciences and Engineering Research Council (NSERC) of Canada. W.Z. was supported by the Beatrice and Vincent Tremaine Fellowship at CITA.

Appendix A Modification to the Search and Vetting Procedure of FM16

We use peerless to complete a fully automated search for long-period planetary transits in our target light curves. In this appendix, we provide a description of the modifications we make to the search and vetting procedure of FM16.

To determine whether a candidate signal is the result of a planetary transit or simply a systematic effect, FM16 fit five separate models to the potential transit. We elect to discard their box model, which was designed to account for signals that are neither convincingly transit-like nor well-described by the other models. We find that this model has little physical motivation behind it and in practice tends to be too stringent, rejecting planetary transits that would otherwise be accepted at all other levels and/or by eye.

We make two supplemental changes to the peerless code to eliminate false positive signals caused in part by the running windowed median (RWM) correction. The RWM correction is applied with a two-day half-width to all PDC light curves (Smith et al., 2012; Stumpe et al., 2012) to reduce the effects of stellar variability – but it does not work in all cases. A subset of light curves display large-amplitude stellar variability on timescales of a few to tens of days, which are not excluded from our target sample because each star’s CDPP is low on the 7.5 hr timescale we consider for our selection cuts. This stellar variability is not fully removed by the RWM with a two-day half-width, and in some cases, the variability is actually made to look more transit-like after the correction (see Figure 11). As such, these signals are not rejected by the variability model described in FM16. To eliminate such signals, rather than use the RWM, we fit a second-order polynomial to the PDC light curves, using a running window with an identical two-day half-width. This polynomial correction considerably decreases the noise of such variable light curves, and thus increases the S/N of planetary transits found in those light curves. However, transits are not well-modeled by a polynomial (the change in flux at ingress and egress is far too abrupt). We therefore revert to the median correction at any point in the light curve where the polynomial residuals are more than $3\sigma$ above their median value.

Unfortunately, this running windowed polynomial correction is very computationally expensive. We therefore choose to only apply it to the light curves of real signals that have already passed all other steps in the vetting process. If any vetted signals are removed by the polynomial correction to the original PDC light curves, we discard the signals as artifacts of the RWM correction applied to stellar variability.

The third modification we make to peerless is to introduce a symmetry criterion to the vetting process for real signals. For each signal, we mirror the data set about the time of mid-transit to produce a ‘left’ transit and a ‘right’ transit. We fit these using the same transit model as FM16, and compare the fitted radii such that

[TABLE]

If the difference between the ‘left’ and ‘right’ radii (scaled by the fitted radius of the original signal) is larger than 0.5%, we consider the signal asymmetric and reject it. This limit is chosen to reflect similar methods of asymmetry-based transit rejection, such as that of Turner et al. (2016) and references therein.

This symmetry criterion is particularly important for removing signals resulting from sudden drops followed by gradual rises in pixel sensitivity, as well as signals that are produced at the union of two sections of the PDC light curve. In these cases, the step function model (or any other model) is insufficient to describe the change in flux well enough to exceed the Bayesian information criterion of the transit model (see FM16 for details). This symmetry restriction therefore serves as a final automated check of the validity of any identified candidates. We also confirm by eye that none of the signals rejected in this step are convincingly transit-like in appearance.

We further correct the factor of $\pi^{1/3}$ missing from peerless in the calculation of transit probability, which resulted in an underestimation of the occurrence rate by 1.46 in FM16.

Additionally, though it is not a modification we chose to adopt in this paper, we deem it worth addressing the fact that we experimented with lowering the S/N threshold for detecting transits. Decreasing this threshold from 25 to 20 introduced a plethora of signals that were not convincingly transit-like by eye, yet could not be removed by the vetting process. We therefore chose to retain the S/N threshold used by FM16 so as to avoid these signals and preserve the automated nature of the search and vetting procedure.

Appendix B Double- versus Single-transit Recovery

An important caveat of the injection and recovery results presented in Section II.4 is that, depending on the injected orbital period and time of mid-transit, a set of simulated transit parameters can produce up to two transits in a single light curve, but only one transit needs to pass all steps of the vetting process to consider the injection recovered. This increases the chance of a double-transit injection being recovered, as it essentially provides twice as many opportunities for the candidate to be accepted. As one might expect, however, a candidate whose orbital period permits it to transit twice within the baseline of observations could masquerade as a candidate with a single-transiting event (i.e., having a much longer period) if the second transit is rejected by the vetting process, takes place during a gap in the observations, or simply is not picked up by the search procedure for any number of reasons. The results of our artificial injections allow us to investigate the prevalence of this phenomenon, assuming a logarithmically flat distribution in orbital period. We may then compare the ratio of double- to single-transiting candidates expected to be detected under this assumption to the ratio observed in our candidate catalog.

First, we note that a candidate can only produce two transits if its orbital period is shorter than the baseline of observations (less than $\sim 4$ yr for Kepler) and if the time of mid-transit is appropriate. Recall also that we only consider orbital periods greater than two years, as this value will be just out of range of transit searches requiring three or more transits within the Kepler baseline. For instance, consider a transiting planet with an orbital period of three years. Its transit could appear twice in a Kepler light curve only if its time of mid-transit occurred in the first quarter of the observations. For simplicity, we restrict our consideration to candidates with orbital periods between 2 and 3.8 yr (coinciding with the first column of each plot in Figure 5), but allow the full range of values for all other transit parameters and stellar magnitudes. We also neglect the effect of decreasing transit probability across this relatively small period range for our simple analysis.

With a logarithmically flat period distribution and our chosen period range, the number of injected planets capable of transiting twice in a Kepler light curve (we will call these double injections) is in a nearly 3:2 ratio with the number of single injections: $\frac{D}{S}|_{\rm inj}=1.56$. Following the application of our search and vetting procedure to these artificial signals, however, the ratio of accepted double- versus single-transit events is very different: $\frac{D}{S}|_{\rm acc}=0.54$. This is because about one-half of all accepted double injections are mistaken for single-transit events. While these injected planets are still considered recovered, their periods are likely to be more crudely estimated (based on their transit duration) compared to double injections where both transits are accepted and a precise period prediction can be made. A breakdown of each step in our double vs. single transit analysis is shown in Figure 12.

How, then, does this compare to the observed ratio of double- versus single-transit events among our planet candidates? Of the 15 candidates we deem genuine with a period range of $2-3.8$ yr and a fractional planet radius of $0.03\leq R_{\rm p}/R_{*}\leq 0.3$, we find four double-transit candidates and six single-transit candidates, such that $\frac{D}{S}|_{\rm obs}=0.67$. This ratio is not totally removed from our injection and recovery results assuming a logarithmically flat period distribution, but it is difficult to confidently compare the two given such a small sample size. Minimally adjusting our chosen range in period or fractional radius substantially changes this $\frac{D}{S}|_{\rm obs}$ ratio while having a much smaller effect on $\frac{D}{S}|_{\rm acc}$ from our injection results.

Therefore, we find it difficult to determine whether a logarithmically flat distribution in period is representative of the population of outer planets in nature. While this assumed distribution is acceptable for the purposes of this work, we recommend further analysis following the inevitable expansion of the long-period planet catalog in the coming years.

Appendix C An Alternative Derivation of the Size Distribution of Outer Planets

The intrinsic size distribution is indicative of the physical compositions of outer planets. The occurrence estimate of Section III.1 already points to the higher incidence of smaller planets. Here, we illustrate our resulting size distribution, using a slightly different argument, independent from the search completeness calculated by peerless.

Adopting the same power law as in equation (5), and taking $\beta=0$, we can constrain the value of $\alpha$ using the observed sample, correcting for the fact that smaller planets are more difficult to detect. Planet candidates in our sample are detected when their transit signal surpasses an S/N threshold of 25. While the signal strength is determined by the transit depth, noise in Kepler light curves is a combination of photon (shot) noise, measurement error (e.g., pointing error and instrument noise), and stellar variability (Koch et al., 2010). The Kepler team encapsulates the total noise of each star into quarterly transit durations of 3, 6 and 12 hr, known as the Combined Differential Photometric Precision (CDPP) (Christiansen et al., 2012). For dwarf stars (nongiants) dimmer than 12th magnitude, Table 2 of Christiansen et al. (2012) shows that the noise is dominated by photon shot noise and the CDPP value scales with stellar magnitude roughly as $CDPP\propto\sqrt{{F}}$, where $F$ is the stellar flux. This fails for brighter stars for which the CDPP is dominated by intrinsic stellar variability (and for which the CDPP is constant). Due to their prevalence, the former population is of relevance for our search. Therefore, the single-transit S/N goes as

[TABLE]

Here, we ignore the dependence on transit duration $\Delta T$ and on stellar radius $R_{*}$. In reality, both affect the S/N to some degree. The minimum planet size at which a transit is detected thus goes as $F^{-1/4}$, meaning smaller planets can only be detected around brighter stars999Indeed, the smallest planet in our sample (KIC 8505215) is detected around the brightest star in our list, at $K_{p}=13.0$. It is also among the smallest stars, making detection easier.. We denote the minimum flux at which a planet of size $R_{\rm p}$ is detected as $F_{\rm min}(R_{\rm p})$. The stellar sample we use has a flux distribution of the form $dN/d\log F\propto F^{\gamma}$ (Figure 2), with $\gamma\approx-1.4$, so the observed size distribution is related to the intrinsic one (equation 5) as

[TABLE]

In other words, because we can only detect small planets around stars that are brighter, a single detection of one such small planet carries a lot of weight. From Figure 13, we find that the observed size distribution scales as $R_{\rm p}^{n}$ where, very roughly, $n\in[1.5,3]$. Therefore, $\alpha=n+4\gamma\in[-4.1,-2.6]$. Using a median value of $\alpha\sim-3.5$ would imply that the occurrence rate rises by a factor of $\sim 16$ between $1~{}R_{\rm J}$ and $0.45~{}R_{\rm J}$. In other words, cold Neptunes are much more prevalent than cold Jupiters; a conclusion similar to that of Suzuki et al. (2016).

Figure 13 also shows the detailed size distribution calculated using detection completeness. Our size distribution above captures the essence of the rising behavior, but is steeper, likely because our simple derivation has ignored effects of stellar variability and variance in stellar radius.