A Closer Look at Exoplanet Occurrence Rates: Considering the Multiplicity of Stars without Detected Planets

Arjun B. Savel; Courtney D. Dressing; Lea A. Hirsch; David R. Ciardi; Jordan P. C. Fleming; Steven A. Giacalone; Andrew W. Mayo; Jessie L. Christiansen

doi:10.3847/1538-3881/abc47d

1. Introduction

For four years, the Kepler spacecraft (Batalha et al. 2010; Borucki et al. 2010; Bryson et al. 2010b; Haas et al. 2010; Jenkins et al. 2010; Koch et al. 2010) pointed at a single patch of sky between the constellations of Cygnus and Lyra, searching for transiting planets orbiting a variety of stars. The survey was groundbreaking—thanks to the Kepler mission, the scientific community has on hand precise photometric data for more than 150,000 stars. These data have thus far yielded 2392 confirmed planets, in addition to 2368 planet candidates that have yet to be confirmed.⁶

The size and uniformity of the Kepler sample were designed to advance the mission objective of constraining the frequency of Earth-like planets around Sun-like stars (Borucki 2016); thus far, data from the Kepler mission have led to the discovery of more than half of all confirmed exoplanets, profoundly informing estimates of exoplanet occurrence rates (e.g., Youdin 2011; Howard et al. 2012; Fressin et al. 2013; Burke et al. 2015; Dressing & Charbonneau 2015; Petigura et al. 2018). Clarifying the properties of the stars observed by Kepler, then, has a direct impact on our understanding of the exoplanet population.

One of the challenges of establishing planet occurrence rates from Kepler data is related to a physical aspect of the spacecraft itself: the size of its pixels. Though significantly smaller than the 21'' × 21'' pixels used by the subsequent NASA TESS mission, Kepler's 398 × 398 pixels were still large enough that individual target stars were not always the only stars within a given pixel or selected stellar aperture. Accordingly, if multiple stars happened to lie within that region, they would contribute to the total recorded flux. While the Kepler pipeline compensated for the flux of known stars near Kepler targets (Bryson et al. 2010a), close-in companions (which are unresolved in Kepler data products) can jeopardize the legitimacy of a possible planet detection.

First, if a target star is part of an eclipsing binary system, the companion star could partially obscure the target star during the course of an orbit, producing a shallow dip in the total system's measured flux that can be mistakenly attributed to a transiting planet. The companion star is gravitationally bound to the target star of interest in this scenario.The second case involves a companion that is near the target star on-sky solely as a result of a line-of-sight projection, with no physical association between the stars. Here, the companion star could itself be eclipsed by yet another unresolved star—known in general as a foreground or background eclipsing binary configuration, with the designation contingent on whether the target star or the eclipsing binary system is closer to the observer. Even in the situation in which an exoplanet truly transits the target star, the "flux dilution" caused by an additional star in the aperture can make the measured transit depth artificially shallow, thereby shrinking the exoplanet's calculated radius (Ciardi et al. 2015). Given the ubiquity of binary stars (Raghavan et al. 2010; Winters et al. 2019) and the size of the Kepler field, constraining stellar multiplicity for every Kepler target star is not possible in a single effort—there is simply too much sky to cover. Furthermore, the faintness of these targets, with half of them being fainter than 14.6 mag in the Kepler bandpass (Huber et al. 2014), poses an additional observational challenge.

A number of studies have sought to address the issue of unresolved binaries in the Kepler field by performing high-resolution, ground-based imaging follow-up of Kepler planet candidate host stars (e.g., Adams et al. 2012; Ziegler et al. 2016; Furlan et al. 2017; Hirsch et al. 2017). In practice, adaptive optics (AO) imaging has proven to be an effective follow-up method, probing regions of parameter space unexplored by spectroscopic methods. In addition, imaging provides useful information for detected companions, such as magnitude differences, position angle, and angular separation (Teske et al. 2015). AO also is generally able to reveal companions at intermediate angular separations of 1'' < Δθ < 4'', where Δθ refers to the angular separation between the primary and companion stars. This is a region complementary to that of speckle imaging, a technique that involves taking a number of short exposures of a target and reconstructing a single image, often by interferometric methods. The resulting data are thereby guarded against noise introduced by turbulence in Earth's atmosphere. Speckle imaging is generally sensitive to nearby stars with angular separations between the diffraction limit⁷ and 12 (Matson et al. 2019). At larger angular separations, sources can fall out of a typical speckle camera's field of view, and the necessary atmospheric isoplanicity breaks down—that is, different regions of the field of view will be distorted to varying degrees, thus reducing the quality of the final, reconstructed image (Horch et al. 2014). The speckle imaging method has thus been useful, for instance, in making statistically robust claims about stars less than an arcsecond from Kepler targets (Horch et al. 2014) and comparing the binary fraction of transit survey target stars to field stars (Matson et al. 2018). AO complements this by allowing one to image stars at greater angular separations that are also important for assessing flux dilution.

The aforementioned imaging studies, however, have not completely characterized the full population of Kepler target stars; most notably, they prioritized stars that are planet candidate hosts. Studying the companion rate of stars that are not known to be planet hosts is necessary to determine accurate occurrence rates, given that the presence of nearby stars affects the depth of transit searches. For instance, in simulating the effects of single-star assumptions on exoplanet occurrence, Bouma et al. (2018) found that calculations based on transiting planet catalogs operating under said assumptions would overpredict exoplanet occurrence, especially for small planets. The full effect of the assumption of target star singularity on estimates of planet occurrence rates has not yet been sufficiently explored because stars without detected planets have not yet been probed with the same follow-up intensity as planet hosts.

The first step of an occurrence rate calculation generally involves the selection of a stellar sample; this step itself frequently includes removing those stars that are suspected to be binaries (e.g., Howard et al. 2012; Fressin et al. 2013; Hsu et al. 2019). Binaries were not culled in the original KIC (Kepler Input Catalog) stellar selection, nor in the selection of Kepler Objects of Interest; this was due in part to the imaging camera used to construct the KIC having a median FWHM of 25, resulting in the blending of close binaries (Brown et al. 2011). Subsequently derived stellar and photometric properties all assumed single target stars. Accordingly, the KIC does not introduce biases related to avoidance of stellar binaries, but the assumed properties of some target stars were distorted by the presence of unresolved binaries.

Some recent estimates (e.g., Petigura et al. 2018) have incorporated the results of imaging analyses such as Furlan et al. (2017) to identify and remove stellar binaries. Biases in the sample selection of these imaging papers, then, can directly influence the stellar and planetary samples considered in occurrence rate studies.

Since the publication of Kepler-era imaging papers, the scientific landscape has been altered by the introduction of data from the Gaia spacecraft (Gaia Collaboration et al. 2016). In Gaia Data Release 2 (DR2), the mission provided high-precision parallaxes and proper motions for more than 1.3 billion sources (Gaia Collaboration et al. 2018). This wealth of data can be used as an independent probe to complement ground-based studies. For example, Berger et al. (2018) revised the radii of 177,911 Kepler stars using Gaia DR2 data in conjunction with the Kepler DR25 data release; using these stellar radius updates, they also revised 2123 confirmed and 1922 candidate exoplanets. Their process of identifying stellar binaries incorporated both ground-based data from Ziegler et al. (2018) and Gaia temperature–radius data (i.e., flagging as binaries those stars that had a larger radius for a given temperature than an equivalent-temperature, main-sequence star).

Throughout this paper, we use several terms to refer to our observed stars. The phrase "target stars" refers to stars that were monitored by Kepler as part of its planet search and subsequently imaged by our team. The primary of a system is the star that is brightest in our observations (obtained in Ks and J). The term companion refers to any fainter star near the primary (i.e., Δθ < 4''); further designation as bound denotes likely physical association.

In this work, we aim to identify and determine the bound status of stars with low angular separation from Kepler target stars. These target stars are not known planet candidate hosts, but the presence of nearby stars would have affected the completeness of Kepler's search for transiting planets and biased the corresponding planet occurrence rate estimates. In Section 2, we present our target sample. We then discuss our data reduction and observations in Section 3. Next, we analyze the physical association of our detected companions to our target stars in Section 4. Finally, Section 5 demonstrates our use of publicly available completeness and reliability tools from the Kepler DR25 data release to explore the broader impacts of this study on estimates of exoplanet occurrence.

2. Target Sample

The broader goal of our work is to refine estimates of the frequency of Earth-like planets by more accurately determining which Kepler target stars would have permitted the detection of potentially habitable planets. We constructed a list of potential targets by using the stellar properties and light-curve noise levels (the Combined Differential Photometric Precision = CDPP) reported in the Kepler Stellar table⁸ to identify stars for which an Earth-sized planet receiving the same insolation flux as Earth would have been detected with a multiple event statistic (MES, the Kepler spacecraft proxy for signal-to-noise ratio, S/N; Burke & Catanzarite 2017) greater than 7.1 (Jenkins et al. 2010). Throughout the paper, we will refer to this group of stars as "amenable Kepler targets." We then selected stars from this list to observe at Lick Observatory (this paper) and Palomar Observatory (J. Christiansen et al. 2020, in preparation). To complement our high-resolution AO imaging campaign, we are also obtaining speckle images of the sample (K. Hardegree-Ullman et al. 2020, in preparation).

In this paper, we discuss the Lick sample, which contains 71 Kepler targets that are not known to host planet candidates. Figure 1 displays the magnitudes, effective temperatures, surface gravities, and metallicities of our targets. Figure 2 displays the positions of our targets in Kepler's field of view. We selected stars that were not known to have associated ground-based AO imaging follow-up data. They were additionally selected on the basis of being visible at a reasonable air mass (<2) from Lick Observatory during our 2018 July observing run. Furthermore, we imposed a faint-end magnitude constraint; any source fainter than roughly 14 Kp mag would have been difficult to observe from Lick, even when employing the laser guide star (LGS). This brightness constraint would preclude us from observing 63% of stars in the Kepler DR25 data release (Mathur et al. 2017). Additional observing constraints included the necessity of pointing more than 20° from the Moon to avoid contamination and keeping the telescope within 45° of zenith when observing in LGS mode.

**Figure 2.** The on-sky distribution of the Kepler DR25 stars (Mathur et al. 2017, white dots) and this paper's target sample (red stars). We highlight star K09776236 (green square), visible in this figure within a circular region of higher-than-average density. This region is in fact the open cluster NGC 6811—notable given the utility of precise ages in exoplanet detection.
Download figure:
Standard image High-resolution image

The Kepler magnitude of our sample ranges from 12.2 to 14.0, with a median magnitude of 13.3. The stellar mass of our sample ranges from 0.1M_⊙ to 0.96M_⊙, with a median mass of 0.67M_⊙. The stellar effective temperature of our sample ranges from 3261 K to 6147 K, with a median temperature of 4870 K.

We note that the metallicity distribution of our sample is discrepant from the metallicity distribution of all amenable Kepler stars. While a full investigation of target selection biases will be undertaken in forthcoming work examining a larger population of stars (J. Christiansen et al. 2020, in preparation), we highlight two considerations: first, our sample of stars may have less reliable metallicity estimates on average than all amenable stars on average, as planet hosts would be subjected to more detailed follow-up characterization. Second, Furlan & Howell (2020) showed that fitting unresolved binary stars as single stars results in generally underestimated metallicities, which could further alter our sample's metallicity distribution.

The stellar mass range of this sample can be motivated by a signal-to-noise ratio (S/N) argument. Transit S/N, to first order, scales as

$\begin{eqnarray*}&&{\rm{S}}/{\rm{N}}\propto \displaystyle \frac{\delta }{\sigma }\sqrt{{N}_{\mathrm{transit}}},\end{eqnarray*}$

with δ being the transit depth (which in turn scales with ${R}_{{\rm{p}}}/{R}_{\star }$ ), σ being the data noise on the timescale of the planetary transit, and N_transit being the number of observed transits. For a planet of a given radius (in our case 1 R_⊕), a search for detectable (i.e., high-S/N) targets favors smaller stars. Stars with small radii that fall within our magnitude ranges are likely to be main-sequence stars with masses less than the Sun's. Applying a cut in MES instead would yield a similar stellar sample, given that the MES is proportional to the above transit S/N expression (Burke & Catanzarite 2017). The explicit effects of our sample biases on our conclusions here are not addressed here, as they will be explored further in an upcoming paper.

The magnitude panel of Figure 1 (top left) displays fewer possible target stars than the other panels because the Huber et al. (2014) catalog contained a number (79) of stars with NaN Kepler magnitudes reported. On average, stars in our target sample tend to be brighter than the full sample of possible target stars observed by Kepler. The majority of our sample stars have T_eff < 6000 K and log(g) > 4.0. The on-sky distribution of our target sample is shown in Figure 2.

3. Adaptive Optics Observations

We conducted our observations at Lick Observatory using the ShARCS camera (Kupke et al. 2012; Gavel et al. 2014) on the Shane 3 m telescope; our targets are presented in Table 1. We observed in both the LGS and natural guide star (NGS) AO modes. Our observing bands were K_S (centered at 2.150 μm with an FWHM of 0.320 μm) and J (centered at 1.238 μm with an FWHM of 0.271). Under ideal conditions for point sources, ShARCS has predicted limiting magnitudes in K_S and J of 19.24 and 22.28 for a detection of S/N = 5 in LGS mode, respectively (Srinath et al. 2014). For our purposes, "ideal conditions" include low atmospheric turbulence, minimal cloud cover, ideal laser performance, and a lack of overresponse of the deformable mirror. For targets with detected companions in K_S, we observed the same target in the J band both to verify that the candidate companions were not chromatic speckles and to measure the colors of bona fide stars to inform subsequent investigations of stellar properties. Integration times varied between 15 and 200 s per frame in K_S and between 20 and 180 s per frame in J. We performed four-point dithers with 4'' spacing over all of our sources to improve image quality. We obtained our initial set of observations of all 71 target stars on 2018 July 21–25. For the four stars that were determined to have a companion only in later analysis, we reobserved the star in both filters; these subsequent observations were taken on 2019 September 5.

Table 1. Constraints on Undetected Companions

Object	Kp mag	Obs. Date	Exp. Time (s)	5σ Detection Threshold ^a
				01	02	05	10	20	40
K02576087	13.7560	2018-07-25	100	0.2	0.4	1.5	3.3	4.9	5.4
K03442846	13.0500	2018-07-25	100	0.1	0.3	1.3	3.2	5.0	6.2
K03640967	13.2730	2018-07-24	100	0.2	0.5	2.1	4.5	6.2	6.8
K03648376	13.9080	2018-07-25	100	0.1	0.2	1.0	2.7	5.2	6.1
K04175216	13.8240	2018-07-23	150	0.4	0.7	2.1	4.2	4.5	6.7
K04458370	12.3860	2018-07-21	70; 50	0.2	0.3	1.4	3.5	5.0	5.6
K04841888	12.3340	2018-07-21	35; 9; 2.91	0.8	1.7	4.3	6.2	7.7	8.4
K05016839	13.1940	2018-07-25	80	0.0	0.0	0.5	1.6	3.3	4.3
K05183442	12.6110	2018-07-21	140	0.3	0.6	1.8	4.0	6.0	6.9
K05184292	13.6470	2018-07-22	150	0.4	0.7	2.2	4.2	5.9	6.4
K05952338	12.7630	2018-07-21	110	0.3	0.6	1.9	3.8	5.9	6.8
K06029168	13.6410	2018-07-22	150; 180	0.1	0.2	1.1	2.8	4.4	4.9
K06292948	13.8360	2018-07-23	150	0.4	0.7	2.4	4.7	5.7	6.2
K06359882	13.3620	2018-07-24	100	0.2	0.4	1.1	2.4	4.9	5.6
K06468660	12.2620	2018-07-24	30	0.4	0.7	2.4	4.4	6.1	6.9
K06511203	13.2500	2018-07-24	140	0.1	0.3	0.9	2.4	4.9	5.7
K06766663	13.8460	2018-07-23	150	0.3	0.7	1.9	4.0	6.3	7.1
K06851483	13.1410	2018-07-24	140	0.2	0.4	1.5	3.4	5.8	6.8
K06947459	13.9880	2018-07-25	100	0.2	0.5	1.6	3.5	4.9	5.4
K07201740	12.2720	2018-07-24	10; 20; 140	0.4	0.7	2.3	4.9	6.2	6.8
K07294931	13.1660	2018-07-24	60	0.4	0.7	2.8	4.9	5.9	6.4
K07304605	13.1230	2018-07-24	60	0.4	0.7	2.5	4.5	6.2	6.9
K07625082	12.2610	2018-07-23	35; 80	0.2	0.5	2.3	4.6	6.2	6.8
K07677767	12.4410	2018-07-22	40; 80	0.1	0.2	1.1	3.1	5.3	7.1
K07692454	12.6610	2018-07-23	2.91; 7.28	0.7	1.4	3.6	5.6	7.5	8.1
K07878293	12.6850	2018-07-21	140	−0.1	−0.2	0.2	1.2	2.4	2.9
K07902000	12.7600	2018-07-23	7.28; 80	0.3	0.5	2.1	4.3	6.6	7.8
K08015478	13.0920	2018-07-24	140	0.2	0.5	1.6	3.6	5.7	6.3
K08040551	12.4440	2018-07-23	2.91; 7.28	0.7	1.4	3.5	5.5	7.0	7.6
K08086729	13.6660	2018-07-22	150; 80	0.4	0.8	2.6	4.9	6.5	7.1
K08108098	13.9520	2018-07-23	80; 150	0.4	0.7	2.1	4.3	5.2	5.7
K08218516	13.9840	2018-07-23	150	0.3	0.7	2.1	4.4	6.3	7.1
K08284195	12.8300	2018-07-21	35	0.4	0.9	2.7	4.1	5.7	6.5
K08301912	12.5610	2018-07-22	40	0.8	1.6	3.9	5.7	6.9	7.5
K08352531	13.6530	2018-07-22	150	0.3	0.7	2.2	4.4	6.2	6.9
K08429578	13.9330	2018-07-25	100	0.2	0.5	1.8	3.8	5.5	6.1
K08479329	13.4180	2018-07-24	140	0.3	0.6	2.2	4.3	5.3	5.7
K08940022	13.8300	2018-07-25	4; 150; 100	0.2	0.4	1.5	3.4	5.6	6.4
K08952590	13.7920	2018-07-25	100	0.2	0.3	1.4	3.2	4.4	4.9
K08959520	13.9250	2018-07-23	7.28; 150	0.4	0.9	2.4	4.6	6.4	6.9
K09019339	12.2370	2018-07-21	30; 45	0.7	1.4	3.5	5.5	6.9	7.4
K09031901	13.9930	2018-07-23	80	0.5	0.9	2.6	4.9	6.6	7.3
K09080730	13.1050	2018-07-24	140	0.3	0.6	1.9	4.4	6.5	7.4
K09084569	12.6420	2018-07-21	140	0.2	0.5	1.7	3.8	6.1	7.1
K09203794	13.7640	2018-07-23	150; 35	0.5	0.9	2.6	4.8	6.6	7.6
K09225647	13.9190	2018-07-23	150; 80	0.5	1.0	2.7	5.1	6.7	7.2
K09266177	13.9390	2018-07-23	35; 150	0.4	0.7	2.2	4.5	5.8	6.3
K09274201	13.7340	2018-07-22	80; 150	0.3	0.6	2.1	4.2	5.4	5.9
K09407276	12.6830	2018-07-24	60	0.3	0.5	2.2	4.5	5.8	6.3
K09459198	12.8820	2018-07-21	60	0.2	0.5	1.8	3.8	5.6	6.4
K09466939	13.6910	2018-07-25	100	0.1	0.2	1.1	2.8	4.4	4.9
K09693806	13.1670	2018-07-25	15; 150	0.0	0.0	0.5	1.7	3.9	5.3
K09760867	13.5400	2018-07-22	150	0.4	0.7	2.4	4.6	6.3	6.9
K09766587	13.8180	2018-07-23	150	0.4	1.0	2.5	4.7	6.5	7.1
K09776236	12.5810	2018-07-24	60; 30	0.6	1.2	3.5	5.5	7.0	7.6
K09787380	13.3970	2018-07-22	150	0.4	0.7	2.2	4.4	6.2	6.9
K09825693	13.7030	2018-07-22	150	0.4	0.9	2.5	4.7	5.9	6.3
K09883594	13.6620	2018-07-25	100; 180	0.1	0.2	0.9	1.9	4.5	5.7
K10063075	12.9820	2018-07-21	35	0.4	0.8	2.5	4.6	5.8	6.4
K10066774	13.1080	2018-07-25	100	0.2	0.3	1.3	3.1	5.1	5.7
K10258023	13.6770	2018-07-22	80; 150	0.2	0.4	1.5	3.2	4.3	4.8
K10412127	13.8800	2018-07-23	150	0.3	0.6	1.9	3.9	5.2	5.7
K10423465	12.4240	2018-07-24; 2018-07-25	30; 80	0.0	0.0	0.5	1.0	1.3	1.7
K10526805	13.5510	2018-07-25	100	0.2	0.4	1.7	3.7	5.0	5.4
K11080894	12.9010	2018-07-21	120	0.3	0.5	1.9	3.8	5.6	6.3
K11086854	13.4290	2018-07-25	80	0.1	0.1	0.7	1.9	3.3	4.0
K11294748	12.5200	2018-07-22	120; 30	0.8	1.7	4.0	5.9	7.1	7.6
K11395188	12.6230	2018-07-22	40	0.7	1.5	3.7	5.6	6.9	7.4
K11550055	12.7960	2018-07-22	40	0.4	0.7	2.3	4.9	6.2	6.8
K11666439	12.9030	2018-07-25	80	0.1	0.1	0.7	2.2	4.1	5.2
K12504315	13.3970	2018-07-25	100	0.3	0.7	2.3	4.5	6.0	6.6

Note.

^aThe threshold at each radius bin has an error of ±0.2 mag.

Download table as: ASCIITypeset images: 1 2

3.1. Data Reduction

We reduced our data using SImMER,⁹ an open-source, Python-based pipeline, first developed and used in Hirsch et al. (2019). Our pipeline employs standard practices for dark subtraction and flat fielding of our science images. We opted for a dither pattern to remove noise from the sky in our observing program. Throughout the data reduction process, the images must be aligned for stacking; for science images, accordingly, the pipeline identifies the center of the brightest star in the image to allow for subsequent alignment. In the image registration mode chosen for this observation program, science images are rotated around points selected within a specified search radius. The initial image is subtracted from each rotated image, and the residuals are summed for rotations of 90°, 180°, and 270°. The center of rotation corresponding to the lowest total residuals is taken to be the center of the image. This pipeline, with algorithms based on Morzinski et al. (2015), will be described in detail in a forthcoming paper (A. Savel et al. 2020, in preparation).

After our science data are passed through our pipeline, we restrict our analysis to an 8'' by 8'' region centered on the primary star. This region is chosen given the image scale on a single Kepler spacecraft pixel (398 × 398),¹⁰ the potential for on-sky drift of stars since the Kepler observations due to stellar proper motion, and the plate scale of ShARCS (0033 pixel⁻¹). Additionally, we discard individual exposures that were taken under particularly detrimental seeing conditions or poor AO performance. To detect sources on the image, we run an automated script that isolates pixel regions with counts greater than 5σ above the background. In line with existing literature (e.g., Marois et al. 2006; Nielsen et al. 2008; Janson et al. 2011; Wang et al. 2015a), we determine this threshold by building concentric annuli out from the central source, calculating the mean and standard deviation of the flux within each annulus, and defining our threshold as 5σ above the mean. Characteristic errors for our detection threshold are computed by calculating 4σ above the mean flux and 6σ above the mean flux. This source-detection algorithm can be triggered by hot pixels, cosmic-ray hits, and other artifacts. Moreover, certain "ghost" images are also present because of secondary reflections in the mirror; these are identified on the basis of their recurring position angle and angular separation from the target, and they (and other artifacts) are ignored in subsequent data processing.

3.2. Detecting Potential Stellar Companions

By running our script on the final, reduced images of each target, we generate a list of primary stars with potential companions. For all sources, including those that do not have a detected companion, we establish constraints on the existence of undetected companions using the aforementioned 5σ detection threshold. For sources that are determined to have potential companions, we calculate Δm (the magnitude difference between the primary and companion) by performing PSF photometry with photutils (Bradley et al. 2019) in each filter.

Unless otherwise noted, we perform PSF photometry on our images. For cases in which PSF fitting is not feasible because of a source's faintness or suboptimal seeing and there is sufficient angular separation between the primary and the companion, we perform aperture photometry on our data. We do so with the aperture photometry capabilities provided by the photutils Python package. The sky count per pixel value is then multiplied by the area of the central source aperture and subtracted from the summed circular aperture counts. Finally, these counts are converted to instrumental magnitudes, which are used to determine Δm between the primary and companions. To account for different magnitude estimates resulting from different aperture sizes (we test aperture radii of 1–32 pixels), we calculate and apply aperture corrections derived from growth curves (Howell 1989; Stetson 1990).

For the majority of our target stars (58 of 71), our observations did not reveal any stellar companions, but we identified 14 companions with Δθ < 4'' for 13 stars (Table 2, Figure 3). Within our sample, 3 (4.2%) stars have companions within 1'', 6 (8.5%) have companions within 2'', and 14 (21.1%) have companions within 4''. Gaia is able to resolve eight of these companions. Our observations of one star in our sample (K05184292) revealed a companion just beyond our 4'' separation cutoff (411 ± 0046). Although Gaia reported a separation of only 382 for the star and companion, we exclude K05184292 from our sample of stars with nearby companions. This observed companion rate is roughly consistent with those of planet host imaging studies (e.g., Horch et al. 2014; Furlan et al. 2017; Hirsch et al. 2017); see Figure 4.

**Figure 3.** Images of the subset of our target sample for which we detected companions. Each box is 4'' × 4'', with ICRS East pointing upward. The color scaling for K03442846, K09766587, K08479329, and the bottom two rows (with the bottom two rows being reobserved stars) is square root in both filters. All other images are scaled linearly. Filters are noted in the bottom left of each image, and the bottom right indicates our physical association designation. Blue boxes are drawn around images that are of the same star; a dotted box is drawn around the target star K08108098, which has two companions.
Download figure:
Standard image High-resolution image

**Figure 4.** A comparison of our data with large surveys of Kepler planet hosts in separation/Δm space. From Hirsch et al. (2017), Δm is represented in Kp mag; from Furlan et al. (2017), in U mag; from Ziegler et al. (2017), in a long-pass filter cutting on at 600 nm; and from Adams et al. (2012) and this work, in K_S mag. While our study differs in its focus on stars not known to host planets, our detections are broadly consistent with those made by these other surveys.
Download figure:
Standard image High-resolution image

Table 2. Detected Companions

Object	Kp mag	Crowding ^a	ΔKp	ΔK_S	ΔJ	Δθ ('')	PA (^o)	Planet-detection Properties (I_p = I_⊕) ^b
								P (days)	t_dur (hr)	Rms CDPP (ppm) ^c	R_p,min (R_⊕) ^d
K03442846	13.050	0.975	1.971	2.957	2.567	2.511 ± 0.046	359.8 ± 0.7	138.5	7.4	57.422	0.87
K03640967	13.273	0.927	0.563	0.501	0.516	0.288 ± 0.046	117.7 ± 9.7	90.9	5.6	66.127	0.72
K03648376	13.908	0.959	0.371	0.371	0.101	0.344 ± 0.046	249.4 ± 5.5	24.4	2.5	122.275	0.39
K04175216	13.824	0.974	2.700	2.226	2.335	1.980 ± 0.046	276.8 ± 1.4	81.2	5.1	100.015	0.81
K06468660	12.262	0.961	4.869	4.964	5.647	2.829 ± 0.046	290.9 ± 0.9	211.8	9.5	31.81	0.87
K08108098	13.952	0.951	1.785	2.918	2.848	3.685 ± 0.046	94.3 ± 1.6	309.1	11.4	50.18	1.01
⋯	⋯	⋯	4.495	3.676	2.458	1.671 ± 0.046	120.8 ± 0.7	⋯	⋯	⋯	⋯
K08284195	12.830	0.989	3.838	2.377	2.803	1.246 ± 0.046	74.0 ± 2.0	307.7	12.6	24.195	0.93
K08301912	12.561	0.986	5.587	6.606	6.140	1.977 ± 0.046	337.9 ± 1.1	270.7	10.6	28.238	0.85
K08479329	13.418	0.979	4.495	3.676	3.965	2.944 ± 0.046	185.2 ± 0.6	228.1	8.7	51.534	0.99
K09766587	13.818	0.964	3.654	3.770	4.251	2.516 ± 0.046	145.7 ± 0.9	95.1	5.8	78.034	0.75
K09883594	13.662	0.983	1.038	0.522	0.645	0.807 ± 0.046	26.6 ± 2.7	97.7	5.7	80.248	0.77
K11294748	12.520	0.991	8.164	5.665	6.399	3.421 ± 0.046	207.0 ± 0.6	241.8	10.3	31.169	0.93
K11550055	12.561	0.988	4.367	4.540	4.188	2.112 ± 0.046	325.7 ± 1.0	133.0	6.3	36.672	0.5

Notes.

^aThe Target Aperture Definition crowding metric (Bryson et al. 2010a) is a per-target metric that indicates the fraction of the flux in the optimal Kepler aperture that can be attributed to the target, as opposed to nearby stars.^bThese are planet properties associated with hypothetical planets; they do not correspond to any known planet or planet candidate.^cThese CDPP values are associated with the hypothetical planet's transit duration (t_dur); however, they do not correspond to the exact transit duration value, as CDPP is provided on a t_dur grid.^dThis minimum value is computed on the basis of detectability; it does not take into account any priors based on planet-formation theories.

Download table as: ASCII Typeset image

To express our measurements with respect to the Kepler bandpass, we utilize the color conversion equations detailed in Howell et al. (2012) to convert our measurements in K_S and J to ΔKp and a predicted Kp mag for companions.

The Kepler pipeline, as previously stated, corrects for the contaminating flux of stars known to be near Kepler targets. For reference, Table 2 also includes the corresponding "crowding metric" for each Kepler star near which we detected a companion. This value represents the Kepler pipeline's estimate of the fraction of flux from the chosen aperture that can be attributed to the target star—values close to 1 indicate low contamination from the background, while values close to 0 indicate high contamination from the background. Notably, all of our stars with detected companions had a crowding metric of greater than 0.92.

For each of our target stars for which we detect a stellar companion, we also compute properties of hypothetical temperate planets that might orbit that target star using the KeplerPORTS code.¹¹ We first compute the period at which a planet would receive the Earth's insolation flux, as per the stellar T_eff, radius, and log(g) reported in the Kepler DR25 catalog (Thompson et al. 2018). This period value allows us to then calculate the duration of a hypothetical temperate planet's transit duration, assuming zero impact parameter (i.e., that the planet transits across the center of its host star). With the corresponding host star, CDPP values taken from NASA Exoplanet Archive's Kepler Stellar table, we can solve for the transit depth of the minimum detectable planet using the following relation from Christiansen et al. (2012):

$\begin{eqnarray*}&&{\rm{S}}/{\rm{N}}=\sqrt{\displaystyle \frac{{t}_{\mathrm{obs}}\,\ast \,{f}_{o}}{P}}\displaystyle \frac{\delta }{{\mathrm{CDPP}}_{\mathrm{eff}}},\end{eqnarray*}$

where t_obs is the total time span over which the target was observed, f_o is the fraction of that time during which the target was actually observed, and CDPP_eff is scaled from the gridded value provided from the Kepler Stellar table. We set the S/N to the minimum detectable 7.1σ to arrive at the minimum detectable transit depth, δ_min.

This leads us to the minimum detectable temperate planet radius for a given star:

$\begin{eqnarray*}&&{R}_{{\rm{p}},\min .}=\sqrt{{\delta }_{\min }}{R}_{\star }.\end{eqnarray*}$

These calculations are performed prior to including any effects of our detected stellar companions; the associated minimum planet sized are determined via the CDPP of the target star's light curve. As such, their radii have not been corrected to account for the flux dilution of nearby stellar companions.

3.3. Detection Limits

In general, our ShARCS observations are sensitive to companions with magnitude differences as large as ΔK_S = 0.3 at 01, ΔK_S = 3.8 at 1'', and ΔK_S = 6.1 at 4''. Table 3 details limiting magnitudes at radii of 01, 02, 05, 10, 20, and 4''; Figure 5 shows detection limits for all nights of observation, and Table 1 presents detection limits on a per-target basis.

**Figure 5.** Contrast curves for all observations, with our detected companions overplotted. The average contrast curve for our 2018 July observing run is shown in black; the gray curves are individual contrast curves. Blue stars represent a bound companion and red squares an unbound companion, while gray circles represent a companion of ambiguous bound status. The green diamond represents the star that Gaia detected but we did not in our AO observations, K09459198. Note that its magnitude difference is expressed in the Gaia G band, and the K_S-band magnitude difference may be closer to the corresponding computed contrast curve.
Download figure:
Standard image High-resolution image

Table 3. Properties of Detected Companions

Object	Gaia Properties					Isochrone Properties of Companion
	Prim. ID	d_prim. (pc)	RUWE_prim	Section ID	d_sec. (pc)	T_eff (K)	R (R_⊙)	M (M_⊙)	σ_diff ^a
K03442846	2052650671329776640	${194.1}_{-0.6}^{+0.6}$	1.020359	2052650671329776768	${4612.7}_{-990.6}^{+1532.4}$	N/A	N/A	N/A	2.64
K03640967	2052895961203523456	${557.5}_{-8.2}^{+8.4}$	1.037085	N/A	N/A	4085	0.32	0.46	0.64
K03648376	2052718737971131648	N/A	N/A	N/A	N/A	3356	0.16	4.97	0.42
K04175216	2076199186744463488	${245.4}_{-1.01}^{+1.019}$	1.0207567	2076199186733958656	${229.8}_{-6.2}^{+6.5}$	3632	0.16	0.25	0.34
K06468660	2075399322092699136	${187.6}_{-1.0}^{+1.0}$	0.947195	2075399326399523200	${1237.7}_{-601.3}^{+3827.7}$	N/A	N/A	N/A	1.53
K08108098	2078110859504453760	${746.6}_{-20.6}^{+21.7}$	2.039871	N/A	N/A	3686	0.18	0.30	3.04
⋯	⋯	⋯	⋯	2078110859504454144	${2380}_{-423}^{+637}$	N/A	N/A	N/A	2.34
K08284195	2105947332818606336	${359.4}_{-2.9}^{+2.9}$	1.122675	N/A	N/A	3655	0.29	0.47	0.20
K08301912	2126348255678300928	${370.5}_{-4.0}^{+4.1}$	0.985490	N/A	N/A	N/A	N/A	N/A	5.20
K08479329	2106721766960549376	${444.3}_{-3.7}^{+3.7}$	1.005669	2106721766957199872	${893.7}_{-150.3}^{+223.2}$	3470	0.13	0.21	0.63
K09766587	2127716151220299392	${209.6}_{-0.7}^{+0.7}$	1.050296	2127716151220582528	${4400}_{-1030}^{+1510}$	N/A	N/A	N/A	1.07
K09883594	2130770216564087552	${397.450695}_{-27.5}^{+31.8}$	11.709280	2130770216559401344	${304.1}_{-8.8}^{+9.4}$	4189	0.36	0.50	0.19
K11294748	2129558962769802112	${247.4}_{-1.3}^{+1.3}$	1.025728	2129558967064981760	${2143.6}_{-796.0}^{1424.7}$	N/A	N/A	N/A	4.34
K11550055	2132074684030519168	${310.9}_{-9.3}^{+9.9}$	4.393098	2132074684030883072	${2849.7}_{-596.6}^{+913.9}$	N/A	N/A	N/A	1.70

Note.

^aThis refers to the number of σ difference in J − K_S between this study's observed companions and their associated primary's recovered model isochrone.

Download table as: ASCII Typeset image

We verified our estimated detection limits by investigating whether Gaia data contains stellar companions to Kepler targets with 1'' < Δθ < 4'' other than the ones that we observed. To do so, for each target star, we performed a cone search with a radius of 20'' around KIC coordinates, allowing for drift due to proper motion. We subsequently calculated the angular separation between all sources in that cone, flagging any pairs with separations of less than 4''. For these pairs, we aimed to identify our target stars by determining whether the brighter source's Gaia G-band measurement was within 2 mag of the KIC magnitude measurement, given the measurement of magnitudes in different wave bands. With one exception, we determined from this analysis that our pipeline caught all companions with 1'' < Δθ < 4'' with respect to our Kepler targets—within Gaia's depth and resolution limits.

Of the 71 Kepler target stars we observed, only one had a companion out to 4'' that our observations missed but was detected by Gaia. This star, a companion to K09459198, is at a separation of 2.60'' from the primary and at a position angle of 196°, with a Δm of 4.64 in Gaia G band. The respective Bailer-Jones distance estimates ( ${149.9}_{-0.4}^{+0.4}$ pc and ${1610}_{-351}^{+583}$ pc) are sufficiently discrepant so as to support the notion of the pair being a chance, unbound alignment. As discussed in Section 4.2, the proper motion comparison method demonstrated that the pair is unbound, as well. It is possible that our source-detection algorithm missed this companion due to its smeared nature, which can be attributed to poor seeing or subpar AO performance (Figure 6).

**Figure 6.** K09459198, the only target star with a companion detected by Gaia but missed in our AO observations. A circle (solid white) with radius determined by the Gaia-determined angular separation from host to companion is plotted over the image that we observed with ShARCS. A second circle (dotted, blue) is overplotted around a candidate for this companion in our image. The image position angle (1905 ± 06) is discrepant from the Gaia-determined position angle (196779 ± 0002); this inconsistency (significant at >3σ) could be attributed either to proper motion drift between the Gaia observation epoch and our observation run or to the faintness of the companion in the image.
Download figure:
Standard image High-resolution image

4. Determining Physical Association

The physical association of stars is important for investigating the role of stellar multiplicity on exoplanet systems: companion stars of the same magnitude and angular separation will dilute exoplanet transits similarly, but only a physically bound companion will affect the formation and evolution of an exoplanet. Additionally, stellar multiplicity has been shown to have a suppressive effect on the formation of exoplanets (Wang et al. 2014a, 2014b, 2015a, 2015b; Kraus et al. 2016). For the purposes of this paper, explicit ties to planet formation are out of scope, as we primarily focus on the impacts of dilution on light-curve observations. The canonical method of ascertaining physical associations for surveys such as this, as described in Section 4.1, is by using relative photometry to determine if a discovered companion lies along the same isochrone as the primary star (e.g., Everett et al. 2015; Teske et al. 2015; Hirsch et al. 2017). In this paper, using data from Gaia DR2 is our preferred method for determining physical association when an AO-detected companion is present in Gaia data, as described in Sections 4.2 and 4.3. However, using both methods links our study to prior ones, establishing the utility of Gaia DR2 in surveys such as these by demonstrating that both provide commensurate results. Moreover, Gaia DR2 cannot recover all the companions that we detect (see Section 4.3), so using the isochrone recovery method allows us to compare our results to past research on subarcsecond companions (e.g., Horch et al. 2014). We are able to perform the isochrone method on all detected companions, whereas we are able to perform Gaia-related analysis for 8 out of 14 stellar companions.

4.1. Isochrone Analysis

Stars that are born together retain similar properties: namely, age and metallicity. With this in mind, one can assess the likelihood that two stars are bound by noting whether or not they lie along the same isochrone, as model isochrones can be parameterized by age and metallicity (assuming other abundances and parameters are held constant).

Hence, in order to determine whether companions are bound to our target stars, we place each target star on an appropriate isochrone (choosing for convenience to use the same isochrones as Huber et al. 2014) and determine whether, based on our observations, the corresponding companion lies along the same isochrone. Following Teske et al. (2015), Everett et al. (2015), Hirsch et al. (2017), and others, we do not derive and fit our own isochrones independently. Rather, we pull the model-dependent values from Huber et al. (2014) for T_eff, [Fe/H], log(g), and the corresponding uncertainties to determine the isochrone model used in the analysis that produced the Huber et al. (2014) catalog. Our choice of the Huber et al. (2014) catalog was motivated by our desire to further cement the linkage between our work and previous, similar studies. Our choice to use [Fe/H] in our procedure as opposed to other metallicity estimates or abundances is thus based on its usage in the initial catalog, which was in turn motivated by narrowband photometry being able to provide some sensitivity to [Fe/H] (Huber et al. 2014). This isochrone analysis procedure is illustrated in Figure 7.

**Figure 7.** A color–magnitude diagram showing a companion determined to be bound by the isochrone method. In each panel, the red cross denotes the absolute magnitude and color of the primary as determined by the recovered catalog isochrone; the blue cross represents the companion, with its position determined by photometric measurements. In the left panel, the primary star is K08284195, and the detected companion is shown to be potentially bound to the primary. Conversely, in the right panel, the primary star is K03442846, and the detected companion is shown to be likely unbound to the primary. In each diagram, the primary and secondary are both plotted with their respective 1σ error bars.
Download figure:
Standard image High-resolution image

Specifically, we generate a grid of Dartmouth isochrones (Dotter et al. 2008) spanning metallicities of −2.5 to 0.5 dex (step size of 0.02 dex) and ages of 1–15 Gyr (step size of 0.5 Gyr) with a mass range of 0.1–4M_⊙ (step size of 0.02 M_⊙). Assuming independence of parameters (which is not strictly true), the catalog values of a given target star and their associated errors are used to define a multivariate distribution. The corresponding probability distribution function is subsequently evaluated at each equivalent evolutionary point (EEP) of each isochrone in T_eff/[Fe/H]/log(g) space. The EEP at which the probability distribution function is at its greatest, then, represents the primary as modeled by Huber et al. (2014). Conveniently, each EEP corresponds to a specific absolute magnitude in a specified wave band. Using this absolute magnitude, we place the primary star, its 1σ error bars in K_S and J − K_S, and its corresponding isochrone in color–magnitude space.

We then use the magnitude differences between the primary and companion derived from our ShARCS photometry to calculate the absolute magnitudes of associated companions in K_S and J, assuming the same distance as the primary. If the companion's 1σ error bars in color–magnitude space intersect the isochrone (allowing an isochrone width of ±1σ in [Fe/H]), we conclude that the companion is likely bound to the primary. In this case, we then identify the EEP on our primary's isochrone that best reproduces the position of the companion in color–magnitude space, allowing for the extrapolation of companion stellar parameters—namely the companion's radius, which is subsequently used to determine exoplanet radius corrections. Otherwise, we assume that the pair represents an unbound, chance alignment.

For those companions that isochrone analysis indicates are bound to their primaries, we determine stellar effective temperature temperature by their placement on said the color–magnitude diagram. This, in conjunction with the notion that the bound companions being fainter than their nongiant primaries, implies that these companions would themselves not be giants, indicates that three (75%) of our bound companions could be affected by model errors. Specifically, the derived temperatures, radii, and masses of these stars could be inaccurate because of inaccuracies in the adopted stellar models, especially with respect to cool dwarfs (e.g., Boyajian et al. 2012; Zhou et al. 2013; Newton et al. 2015). These shortcomings would affect nine primary stars (13%) of our sample, per the classification schema of Pecaut & Mamajek (2013). Our justification for continuing to use these models is that we simply seek to recover the isochrones used to populate our sample's catalog stellar parameters, as opposed to independently determining best-fit parameters. In doing so, we remain internally consistent with regard to errors incurred in the catalogs.

4.2. Proper Motion Comparison

Comparing the proper motions of stars is an established practice for determining physical association (e.g., Luyten 1988; Lépine & Bongiorno 2007; Deacon et al. 2015; Janes 2017; Godoy-Rivera & Chanamé 2018). The overarching goal of these studies has been to identify pairs of stars that occupy similar regions in phase space with respect to the background, as these stars are likely to be physically associated. Simply put, stars that move together were likely born together.

Here, we opt to implement a modified approach to the method described in Lépine & Bongiorno (2007) to identify whether AO-detected companions present in Gaia data are physically associated. As in the aforementioned studies, Lépine & Bongiorno (2007) sought wide binaries in proper motion catalogs by identifying "common proper motion pairs" (CPM pairs). The catalog of choice for Lépine & Bongiorno (2007) was the LSPM-north proper motion catalog (Lépine & Shara 2005). In particular, Lépine & Bongiorno (2007) aimed to find LSPM-north CPM pairs in which at least one star was represented in the Hipparcos catalog (Perryman et al. 1997)—with the goal of constraining the projected physical separation between stellar components with a Hipparcos parallax measurement. To do so, Lépine & Bongiorno (2007) constructed all possible pairs between two groups: those stars that are represented in both LSPM-north and Hipparcos, and the entire LSPM-north catalog. For each of these pairs, Lépine & Bongiorno (2007) calculated their angular separation (Δθ) and proper motion difference (Δμ) and constructed a distribution with respect to a variable dependent on both terms, ΔX = ΔX(Δθ, Δμ). Lépine & Bongiorno (2007) then performed this process once more—this time, creating a simulated population introducing an angular offset in every calculation of the angular separation for each pair. This offset disrupts the distribution of "genuine" CPM pairs—as physically associated pairs are now much farther on the sky than they were before—without perturbing the overall statistics of chance alignments. The distributions of ΔX for both the initial and simulated populations are then compared against one another. The threshold for binarity is taken to be the value of ΔX at which ΔX_sim > ΔX_initial—that is, the value of ΔX greater than which the chance alignment pairs are more common than genuine CPM pairs. Every initial pair's calculated ΔX can thus be analyzed: if its value of ΔX is less than the threshold, then it is more likely to be a bound, genuine CPM pair than a chance alignment. Otherwise, it is more likely to be a chance alignment than a bound, genuine CPM pair.

The analysis performed by Lépine & Bongiorno (2007) was tailored for stars with large proper motions, with the rationale being that they are often close to the observer, meaning that intrinsically faint companions to these stars would in turn often be resolved in imaging surveys. Additionally, Lépine & Bongiorno (2007) note that, at the time, this restriction would allow for greater representation of stars with accurate parallax measurements in the sample. For our approach, we need not limit ourselves to stars with large proper motions, as we are already restricting ourselves to stars that we have resolved with our AO data—our magnitude limitations already implicitly restrict the maximum distance (and thus the minimum proper motion) of our stars. Moreover, the uniformity and precision of Gaia DR2 data ensures high-quality parallax measurements even for sources with submilliarcsecond/year proper motion in the Kepler field.

Notably, Lépine & Bongiorno (2007) performed their analysis for every pair between stars represented in both LSPM-north and Hipparcos and the stars in LSPM-north catalog. For the Gaia DR2 data, this is not feasible at once; at best, ignoring duplicates, the problem scales as Θ(n² − n), thereby implying billions of pairs across the full Kepler field. As such, we employ an adaptation of the method described in Lépine & Bongiorno (2007) to select smaller, representative subsamples.

Usage of Gaia DR2 data in phase-space characterizations necessitates other qualifiers. First, the Gaia team notes that DR2 data alone cannot resolve sources with Δθ < 04–05 (Arenou et al. 2018). Indeed, comparisons to Robo-AO binary yields revealed a DR2 recovery rate (that is, the fraction of Robo-AO-observed binaries that are also represented in Gaia DR2) of 22.4% within Δθ < 1'', as opposed to a recovery rate of 93% for Δθ > 2'' (Ziegler et al. 2018). Consequently, our detected companions at subarcsecond separations are in general unlikely to be represented in Gaia DR2 data.

Epoch propagation with Gaia data for our sample is not possible for two reasons. First, propagating forward KIC-blended binaries to Gaia DR2 source locations is out of scope for this paper; it is one of the concerns that is being factored into the next Gaia data release. Second, we cannot propagate forward from Gaia DR2 sources to our observation times, as Gaia's native ADQL epoch propagation command requires radial velocity components for all sources in the selected region of sky, which are not available.

Our proper motion comparison method is as follows:

1.
We take a slice of Gaia sources in the sky with a size commensurate with the Kepler field, centered on R.A. = 19^h22^m40^s and decl. = +44°30'00''. We fit a piecewise power-law function to the number density of all the stars in this region as a function of proper motion; the break is empirically determined to be about 8 mas yr⁻¹, though the fits of the functions on either side of the break are robust to changes on the order of 0.5 mas yr⁻¹. Our empirical fit to values greater than the break has the form N(μ) = μ^−3.79 for proper motion μ (Figure 8). This step is performed only once.
2.
We next perform a smaller cone search around each primary. The radius of this cone search varies on a case-by-case basis between 03 and 07. As clarified below, the search radius is chosen in the interest of computational tractability.
3.
Each source in this second cone search is then paired with every other source, producing n² − n initial pairs for n sources in the cone search. For each pair, Δθ is first calculated; given that the sources that we are interested in are close to one another on the sky, we make use of the approximation ${\rm{\Delta }}\theta \simeq \sqrt{{({\alpha }_{1}-{\alpha }_{2})}^{2}{\cos }^{2}{({\delta }_{1})+({\delta }_{1}-{\delta }_{2})}^{2}}$ for R.A. and decl. (α_i, δ_i). If Δθ > 4'', the pair is rejected, as we are uninterested in the rate of chance alignments beyond this radius (see Godoy-Rivera & Chanamé 2018 for a similar treatment on Kepler CCDs). Otherwise, Δμ and μ_average are then calculated.
4.
The previous step is then repeated, this time after introducing an angular offset Θ in the calculation of Δθ for every pair within the second cone search. In practice, this amounts to subtracting Θ from each R.A. term, such that ${\rm{\Delta }}\theta \simeq \sqrt{{({\alpha }_{1}-[{\alpha }_{2}-{\rm{\Theta }}])}^{2}{\cos }^{2}{({\delta }_{1})+({\delta }_{1}-{\delta }_{2})}^{2}}$ , per Lépine & Bongiorno (2007). Doing so effectively "shuffles" sources within the second cone search, preserving the overall statistics of unbound companions while disrupting the statistics of bound companions. Accordingly, this population (hereafter denoted with a "sim" subscript to relevant variables) can be compared with the unmodified Gaia population (hereafter denoted with a "Gaia" subscript to relevant variables). We follow Godoy-Rivera & Chanamé (2018) in setting Θ = 10''.
5.
A variable ΔX is defined to model the probability that a pair represents a genuine binary. Functionally, it assumes the form ${\rm{\Delta }}X=\sqrt{{\mu }^{-3.79}{\rm{\Delta }}\theta {\rm{\Delta }}\mu }$ , per Lépine & Bongiorno (2007). Hereafter, ΔX_sim = ΔX_sim(Δθ_sim, Δμ) for the simulated ("shuffled") population and ΔX_Gaia = ΔX_Gaia(Δθ_Gaia, Δμ).
6.
We next calculate ΔX_sim and ΔX_Gaia for each possible pair in the cone search and overplot the binned respective number densities, N(ΔX_sim) and N(ΔX_Gaia) (see the left panel of Figure 9).
7.
We compare N(ΔX_Gaia) to N(ΔX_sim). The value of ΔX at which the latter exceeds the former is taken to be our threshold for physical association, ΔX_thresh.
8.
Our remaining steps with respect to Gaia involve placing our AO targets in the context of these greater patches of sky. For each target, we take a cut around each primary with a radius of 20'', locating it and its putative companion by matching the Δm from our photometry. If multiple sources meet this criterion, we attempt more stringent matches based on position angle and Δθ. These comparisons, however, are less desirable than the Δm comparison, as they are more liable to change on the timescale of a few years.
9.
Finally, we calculate ΔX for our primary and companion and determine where the pair lies with respect to ΔX_thresh (Figure 9). If ΔX for our pair is less than ΔX_thresh, the pair is likely to be bound; otherwise, the pair is likely to be unbound.
10.
We repeat steps 2–8 for each primary and companion from our observations.

**Figure 8.** Gaia DR2 proper motion of all sources in the Kepler field, with draws from our MCMC fit to these data overplotted. We find a break in the power-law distribution of the number density of Gaia sources binned by proper motion, with there being very few sources at the high and low proper motion tails. The functional form of this distribution is taken into account when constraining the physical association of stars with our proper motion method (Section 4.2).
Download figure:
Standard image High-resolution image

**Figure 9.** An example of the proper motion comparison method for the target K11550055 and its companion, which we observed to be 4.37 K_S magnitudes fainter than the primary and have Δθ = 206. Our analysis indicates that the target and the companion that we detected near it are not likely to be physically associated. This is determined based on the similarity of the target and companion's proper motions and positions, quantified by the variable ΔX (defined as $\sqrt{{\mu }^{-3.79}{\rm{\Delta }}\theta {\rm{\Delta }}\mu }$ ). The variable N refers to the number of pairs lying within a specific bin of ΔX. The distributions of ΔX are calculated for both the population near the companions and target on the sky and a simulated, or "shuffled," version of that population (left). The value at which N_sim > N_Gaia is ΔX = 1.5 × 10⁴, so we adopt this as our threshold for physical association (right). Because ΔX for K11550055 and its companion is greater than this threshold, we determine that the pair is unlikely to be physically associated. It should be emphasized that this boundary is target dependent, as it varies based on the surrounding on-sky population.
Download figure:
Standard image High-resolution image

**Figure 9.** An example of the proper motion comparison method for the target K11550055 and its companion, which we observed to be 4.37 K_S magnitudes fainter than the primary and have Δθ = 206. Our analysis indicates that the target and the companion that we detected near it are not likely to be physically associated. This is determined based on the similarity of the target and companion's proper motions and positions, quantified by the variable ΔX (defined as $\sqrt{{\mu }^{-3.79}{\rm{\Delta }}\theta {\rm{\Delta }}\mu }$ ). The variable N refers to the number of pairs lying within a specific bin of ΔX. The distributions of ΔX are calculated for both the population near the companions and target on the sky and a simulated, or "shuffled," version of that population (left). The value at which N_sim > N_Gaia is ΔX = 1.5 × 10⁴, so we adopt this as our threshold for physical association (right). Because ΔX for K11550055 and its companion is greater than this threshold, we determine that the pair is unlikely to be physically associated. It should be emphasized that this boundary is target dependent, as it varies based on the surrounding on-sky population.
Download figure:
Standard image High-resolution image

4.3. Distance Estimates

One of the unique advantages of Gaia DR2 is that it allows for unparalleled accuracy in stellar distance estimates for large swaths of sources. We would expect the separation between any physically associated components to be less than 1 pc given the extant measured frequencies of very wide binaries (e.g., Lépine & Bongiorno 2007; Raghavan et al. 2010; El-Badry et al. 2019). Thus, distance estimates based on Gaia DR2 can be used to provide a final discriminatory check against spurious results of the isochrone and proper motion criteria. When making use of Gaia DR2 data to determine distances to sources, we use the Bailer-Jones (2015) distance estimates, which benefit from a full Bayesian treatment; accordingly, they take into account negative (and other spurious) parallax values to provide more accurate distance measurements.

4.4. Radius Correction Factors

As stated above, the effect of flux dilution depends on the orbital configuration of the system. Determining which star a planetary candidate orbits is out of scope for this paper; as such, we present our calculations of radius correction factors for each valid orbital configuration. While the stars in our target sample are not currently known to host exoplanet candidates, the radius correction factors that we calculate would need to be applied to any exoplanets discovered to orbit these stars through Kepler data products. Importantly, these correction factors also affect estimates of vetting completeness required for occurrence rate calculations (see Section 5).

Similarly to Ciardi et al. (2015), we segment our sample of stars with companions into the following groups: (1) pair (either unbound or bound) in which the planet is orbiting the primary, (2) bound pair in which the planet is orbiting the secondary, (3) unbound pair in which the planet is orbiting the secondary, (4) triple in which the planet is orbiting the primary, (5) triple in which the planet is orbiting the (bound) secondary, and (6) triple in which the planet is orbiting the (unbound) secondary, (7) triple in which the planet is orbiting the (bound) tertiary, and (8) triple in which the planet is orbiting the (unbound) tertiary. Ciardi et al. (2015) do not analyze unbound systems; Hirsch et al. (2017), however, note that the radius correction factor cannot be reliably calculated for groups in which the planet orbits a companion star that is not bound to the primary (groups 3, 6, and 8). This is due to the fact that the color of the unbound companion would be influenced by unknown stellar type, a factor complicated by unknown degrees of interstellar extinction. Although Gaia provides line-of-sight extinction and reddening for sources with Gaia G < 17 (Gaia Collaboration et al. 2018), none of our unbound companions recovered by Gaia fall within that magnitude range.

Per Furlan et al. (2017), our radii conversion equations are as follows:

$\begin{eqnarray*}\begin{array}{rcl}{R}_{{\rm{p}},\mathrm{corr}1} & = & {R}_{{\rm{p}}}\sqrt{1+{10}^{-0.4{\rm{\Delta }}m}}\\ {R}_{{\rm{p}},\mathrm{corr}2} & = & {R}_{{\rm{p}}}\displaystyle \frac{{R}_{\sec }}{{R}_{\mathrm{prim}}}\sqrt{1+{10}^{0.4{\rm{\Delta }}m}}\\ {R}_{{\rm{p}},\mathrm{corr}4} & = & {R}_{{\rm{p}}}\sqrt{1+{10}^{-0.4{\rm{\Delta }}{m}_{\sec }}+{10}^{-0.4{\rm{\Delta }}{m}_{\mathrm{tert}}}}\\ {R}_{{\rm{p}},\mathrm{corr}5} & = & {R}_{{\rm{p}}}\displaystyle \frac{{R}_{\sec }}{{R}_{\mathrm{prim}}}\sqrt{1+{10}^{0.4{\rm{\Delta }}{m}_{\sec }}(1+{10}^{-0.4{\rm{\Delta }}{m}_{\mathrm{tert}}})}\\ {R}_{{\rm{p}},\mathrm{corr}7} & = & {R}_{{\rm{p}}}\displaystyle \frac{{R}_{\mathrm{tert}}}{{R}_{\mathrm{prim}}}\sqrt{1+{10}^{0.4{\rm{\Delta }}{m}_{\mathrm{tert}}}(1+{10}^{-0.4{\rm{\Delta }}{m}_{\sec }})}.\end{array}\end{eqnarray*}$

As noted by Furlan et al. (2017), these equations are only valid with respect to magnitude measurements in the Kepler bandpass. As such, we make use of our predicted Δm_Kp in calculating radius correction factors.

Our target sample's radius corrections are presented in Table 4; the mean radius correction factor is 1.10, assuming that the planet candidate (if it exists) orbits the primary star. For stars that have a bound companion, the mean radius correction factor is 1.78, assuming that the planet candidate orbits the companion. In accordance with the care observed by Furlan et al. (2017), we caution against applying these average correction factors to a population at large, given the high standard deviations (0.17 for planets orbiting the primary and 0.31 for planets orbiting the secondary) in the calculations.

Table 4. Likelihood of Physical Association and Derived Planet Radius Corrections

	Bound Criteria
Object	Proper Motion	Isochrone	Distance	Designation	R_p,corr, Orbit Prim	R_p,corr Orbit Sec
K03640967	N/A	True	N/A	Bound	1.2738	1.7781
K03648376	N/A	True	N/A	Bound	1.641263	1.3926
K08108098	N/A	True	N/A	Bound	1.0996	1.68932
K08284195	N/A	True	N/A	Bound	1.01099	2.26115
K03442846	False	False	False	Unbound	1.0793	N/A
K08108098	False	False	False	Unbound	⋯	N/A
K08301912	N/A	False	N/A	Unbound	1.00198	N/A
K08479329	False	True	False	Unbound	1.007930	N/A
K09766587	False	False	False	Unbound	1.0154	N/A
K11294748	False	False	False	Unbound	1.00204	N/A
K11550055	False	False	False	Unbound	1.0170	N/A
K04175216	False	True	True	Ambiguous	1.0411	1.3417
K06468660	True	True	True	Ambiguous	1.0027	1.5962
K09883594	N/A	True	N/A	Ambiguous	1.1809	2.1463

Download table as: ASCII Typeset image

4.5. Physical Association of AO-observed Stars

Of our detected companions, two within 1'' are bound (67%) and three within 4'' are bound (21%). See Table 3 for the properties of these companions as estimated using the procedure described in Sections 4.1–4.3. For our purposes, we designate a pair of stars as bound if their distance estimates lie within 1σ of one another and the pair meets at least one other criterion for physical association: isochrone analysis or proper motion comparison (see Table 4). In all cases, the isochrone criterion is available for analysis. If Gaia cannot recover both the target star and the candidate stellar companion (and thus proper motion comparison and distance estimate comparisons are not possible), our bound designation is based solely on the isochrone criterion; corresponding cells in Table 4 are denoted by "N/A." None of our target stars with bound companions are noted as binaries based on the Berger et al. (2018) "binary = 1" or "binary = 3" flags, both of which take into account the position of a star on a Gaia stellar radius–T_eff plot. This is expected for our more distant companions that Gaia resolved, as the stellar radius–T_eff analysis applies to binaries that are represented as single sources within Gaia data.

5. Impact on Occurrence Rate Calculations

Using our data, we conduct a preliminary investigation of the role of contamination from nearby stars on previous estimates of the frequency of potentially Earth-like planets (Figure 1). To do so, we turn to the recently released Kepler DR25 completeness and reliability products. Here, constraining "completeness" requires understanding what fraction of exoplanets in nature have been detected, while constraining "reliability" requires understanding to what extent every entry in an exoplanet catalog represents a true exoplanet. We first follow the method described in Bryson et al. (2020) to produce a uniform stellar catalog that is amenable to completeness and reliability characterization, with updated stellar radius data based on Gaia DR2. Following Bryson et al. (2020), we clean the set of GK stars observed by Kepler, removing stars on a variety of conditions: those with poor Gaia DR2 astrometric fits; those that are likely binaries as determined by Berger et al. (2018), not this study; evolved stars; noisy stars—i.e., those identified as exhibiting systematic brightening events on the timescale of a planetary transit (Burke & Catanzarite 2017); those without provided limb-darkening coefficients; those with a duty cycle <60% or without duty cycles altogether; those indicating an undesirable amount of data removal in the process of completing a transit search and removing known transit signals; those with a data span of less than 1000 days; and those with the timeoutsumry flag ≠ 1—i.e., those for which the Kepler pipeline was unable to search for transits at all desired transit durations without a timeout issue. We refer to the resulting stellar catalog as our "cleaned" version of the Kepler Input Catalog.

After refining its input stellar catalog, the Bryson et al. (2020) pipeline then constrains the vetting completeness of the automated planet search conducted using the Robovetter code (Thompson et al. 2018). This tool was developed to promote periodic signals in Kepler data from "threshold-crossing events," or "TCEs," to the more stringent designation of "planet candidates," or "PCs." Vetting completeness is expressed as the fraction of TCEs that the Robovetter correctly classifies as PC. This fraction is constrained by modeling a successful vetting as a draw from a binomial distribution with parameters that are fitted for in a Bayesian framework. The second aspect of completeness—detection completeness, or the degree to which the Kepler spacecraft itself was able to detect exoplanets—is addressed next. This function is dependent on per-star expected MES for exoplanets of differing properties.

Subsequently, the pipeline constrains the per-planet reliability, with respect to both astrophysical false positives and their instrumental analogs, "false alarms" (i.e., thermal responses of CCDs, cosmic-ray hits). With respect to astrophysical false positives (e.g., grazing eclipsing binaries), the Bryson et al. (2020) pipeline makes use of the probabilistic framework described in Morton et al. (2016), which in turn relies on the vespa code (Morton 2015) to encapsulate the probabilities of a variety of false-positive scenarios with a single score. Centroid shifts during transit-like events, which have been used to determine whether a transit-like event originates from the primary star or a companion (e.g., Bryson et al. 2013; Dressing et al. 2014; Barclay et al. 2015), are not considered by vespa. Constructing a holistic, quantitative estimate of reliability involves parameterizing the Kepler pipeline's false-positive effectiveness (its ability to correctly identify false positives) and its observed false-positive rate. From these constraints, the number of true exoplanets per star is computed by Bayesian inference, expressed in terms of orbital period and exoplanet radius.

Each major intermediary step in the Bryson et al. (2020) pipeline (detection completeness, vetting completeness, reliability) involves fitting parameters of models that describe each function. This is done within a Markov Chain Monte Carlo (MCMC) formalism, resulting in probability distributions for each parameter of each constraining function. For example, during the vetting completeness portion of the Bryson et al. (2020) pipeline, the best-fitting parameterization of how the rate function describes the probability of successful vetting is determined with MCMC. In the interest of making a clearer comparison to the previous estimates by Bryson et al. (2020), we do not make significant adjustments to the pipeline. We do, however, lengthen most MCMC chains used; all chains were initially run for 2000 steps, but we run all chains for 200,000 steps, except for one chain that was run for 10,000 steps. We choose these chain lengths to ensure that the integrated autocorrelation time (Goodman & Weare 2010) and Gelman-Rubin statistic (Gelman & Rubin 1992) imply chain convergence.

We assess the relevant impact of our study as follows. For each star in the "cleaned" Kepler Input Catalog, we must determine if the star can host an Earth-like planet that is above the detection threshold—determined to be 7.1σ (Jenkins et al. 2010). We define "Earth-like" planets as meeting the criteria R_p = R_⊕ ± 0.2R_⊕ and P_p = P(I = I_⊕) ± 0.2P(I = I_⊕), with I_⊕ being the insolation received by Earth from the Sun. This results in a range of insolation fluxes of 1.35I_⊕ to 0.78I_⊕. To calculate the detectability of a planet with a given period, radius, and host star, we make use of the KeplerPORTs code. With KeplerPORTs, we access the MES of simulated planets that were injected into the light curves of stars observed by Kepler and recovered by the Robovetter during testing. These simulated planets were injected at regular intervals of planetary parameters. To determine how a hypothetical Earth-like planet would have fared in these injection and recovery tests, we then interpolate along existing points in MES/orbital period space. If the estimated MES exceeds 7.1, we flag the primary star as amenable to the detection of potentially Earth-like planets (Figures 10–13).

5.1. Impact on Completeness

When radius correction factors are taken into consideration, some Kepler targets that initially appeared amenable to the detection of Earth-like planets will no longer be. We quantify this effect by computing MES values with respect to the period and radius of potentially habitable planets in the vicinity of Earth's region of parameter space, then repeating the process using planet radii corrected for flux dilution. This essentially penalizes the injection depth—simulating the diluting effect of a stellar companion. In our sample, three stars—K03442846, K03640967, and K09883594—were initially assumed to be amenable to the detection of an Earth-sized planet with a 1 yr orbital period but are not when the corresponding radius correction factor is taken into account (see Figure 12). However, these low-mass stars are still amenable to the detection of an Earth-sized planet receiving the same insolation flux as Earth: K03442846 in the period range of 111 days to 166 days, K03640967 between 73 days and 109 days, and K09883594 between 78 and 117 days.

**Figure 10.** Hertzsprung–Russell diagram of the Kepler DR25/Gaia DR2 overlap stars used to estimate this study's impact on occurrence rates (red, black) and our target sample (beige). Contours are overplotted in dense regions for the DR25/Gaia DR2 overlap stars, and normalized marginal histograms for each group in each parameter are presented as well. The contour threshold is 30 counts per 2D bin, of which there are 40; contours are based on the log₁₀ number of stars per bin. We adopt the same radius cut as Bryson et al. (2020) in the interest of preventing a transit duration from exceeding 15 hr, which was the longest transit duration searched by the Kepler pipeline.
Download figure:
Standard image High-resolution image

**Figure 11.** Temperature–magnitude plot of the DR25/Gaia DR2 overlap stars used to estimate this study's impact on occurrence rates (red, black) and our target sample (beige). We make use of contour parameters as in Figure 10 for the Kepler DR25/Gaia DR2 stars, overplotting scatter points in less dense regions; for the amenable Kepler targets, we only plot scatter points.
Download figure:
Standard image High-resolution image

**Figure 12.** MES levels for K03640967 with respect to planetary orbital period and radius with this paper's analysis taken into account. Mathur et al. (2017) list the target as having an effective temperature of ${4398}_{-87}^{+78}$ K (a K-type star, per Pecaut & Mamajek 2013), a radius of ${0.579}_{-0.022}^{+0.034}{R}_{\odot }$ , a mass of ${0.609}_{-0.035}^{+0.026}{M}_{\odot }$ , and [Fe/H] of $-{0.44}_{-0.15}^{+0.15}$ dex. Each detection limit contour is drawn at an MES level of 7.1. Without a radius correction factor applied to system, the MES of an Earth-sized planet orbiting this star with a period of 365 days would be 10.85—above the detection threshold. With the primary radius correction applied, a planet with those same parameters would have an MES of 6.95—below the detection threshold. With the secondary radius correction factor applied, an Earth-sized planet receiving Earth's insolation flux (a period of 91 days, the dotted gray line in the figure) would be below the detection threshold, as well.
Download figure:
Standard image High-resolution image

**Figure 12.** MES levels for K03640967 with respect to planetary orbital period and radius with this paper's analysis taken into account. Mathur et al. (2017) list the target as having an effective temperature of ${4398}_{-87}^{+78}$ K (a K-type star, per Pecaut & Mamajek 2013), a radius of ${0.579}_{-0.022}^{+0.034}{R}_{\odot }$ , a mass of ${0.609}_{-0.035}^{+0.026}{M}_{\odot }$ , and [Fe/H] of $-{0.44}_{-0.15}^{+0.15}$ dex. Each detection limit contour is drawn at an MES level of 7.1. Without a radius correction factor applied to system, the MES of an Earth-sized planet orbiting this star with a period of 365 days would be 10.85—above the detection threshold. With the primary radius correction applied, a planet with those same parameters would have an MES of 6.95—below the detection threshold. With the secondary radius correction factor applied, an Earth-sized planet receiving Earth's insolation flux (a period of 91 days, the dotted gray line in the figure) would be below the detection threshold, as well.
Download figure:
Standard image High-resolution image

To assess the broader impact of our study on the planet occurrence rate estimates based on Kepler data, we begin with those stars that were amenable to the detection of Earth-like planets and probabilistically flag a subset based on our observed companion rate. For these flagged stars, we apply a radius correction factor. As noted by Ciardi et al. (2015), applying the median radius correction factor to all flagged stars would not be prudent, as the distribution of radius correction factors is decidedly non-Gaussian. Consequently, to gain a more comprehensive picture of how the presence of nearby stellar companions affects the transit search depth for Kepler target stars, we perform Gaussian kernel density estimation on our distribution of radius corrections for the case in which an exoplanet orbits the primary.

For each flagged star, we next draw radius correction samples from the estimator. A benefit of using this method is that artificial gaps in the distribution due to binning and small number statistics are essentially smoothed over, allowing us to sample a continuous distribution that is likely to be more representative of what a larger, deeper, and higher-resolution survey might yield. One pitfall of the immediate application of this method, though, is that the Gaussian nature of the estimator causes some draws to be unphysical. Namely, because the distribution peaks near 1.0, draws from the estimator will often produce values less than 1.0, as Gaussian distributions are symmetric about their peak—the opposite of the intended flux dilution effect. This issue is ameliorated by enforcing a reflective boundary at 1.0, such that any draw from the estimator less than 1.0, D, is mapped to 2 − D. Given the aforementioned symmetry of Gaussian distributions, this adjustment produces the desired effect of a peak near the minimum value while retaining the benefits.

Once a radius correction factor is assigned to a star, we use the KeplerPORTs code to determine whether the star is amenable to the detection of an Earth-like planet. This process is repeated for all stars in the "cleaned" stellar sample, producing a list of stars that, even when considering flux dilution effects, would be amenable to the detection of Earth-like planets. Finally, we compare this list of stars to the list of stars amenable to the detection of Earth-like planets without accounting for flux dilution to determine the overall effect of ignoring flux dilution on estimates of planet occurrence.

Averaging the results of this approach over 10 iterations, we find that 7.8% ± 0.2% of stars that were amenable to the detection of Earth-like planets would not be after applying representative radius correction factors based on our observed companion rate.

5.2. Impact on the Frequency of Potentially Habitable Planets

To assess our study's impact on the estimates of the frequency of potentially habitable planets, we begin with the "cleaned" set of Kepler DR25 stars. We next draw from a binomial distribution informed by our observed companion fraction to determine whether to assign a companion to a given star in that smaller subset—and, therefore, whether to withhold this star from the stellar sample used for our occurrence rate calculation. This smaller, refined stellar sample is then passed through the Bryson et al. (2020) pipeline discussed in Section 5, and its computed occurrence rates are compared to the baseline results.

For various definitions of the region of parameter space corresponding to "Earth-like" planets (e.g., Hsu et al. 2019; Zink & Hansen 2019; SAG13;¹² see Figure 13), we arrive at an average of a 6% increase in occurrence around GK stars (see Table 5). For example, we increase the estimate for the SAG13 occurrence rate from ${0.15}_{-0.07}^{+0.11}$ to ${0.16}_{-0.07}^{+0.12}$ —an increase of 6% (0.9σ), but consistent with the previously reported errors and likely not discernible between occurrence rate studies.

**Figure 13.** A depiction of the approximations to η_⊕ considered in this work with respect to Earth's flux from the Sun and its radius. All confirmed planets with known radii (green squares) are shown to provide a sense of scale. Regions used to compute estimates of η_⊕ are also shown. The "×" marks ${{\rm{\Gamma }}}_{\oplus }=\tfrac{{\partial }^{2}f}{\partial \mathrm{log}P\partial \mathrm{log}R}$ , evaluated at R = R_⊕, P = 365 days. The teal box represents the habitable zone adopted by Hsu et al. (2018), equivalent to 237 < P(days) < 500, 1.0 < R(R_⊕) < 1.75. The black box represents the habitable zone adopted by Zink & Hansen (2019), equivalent to 222.65 < P(days) < 808.84, 0.72 < R(R_⊕) < 1.7. The maroon box represents the SAG13 definition of η_⊕, equivalent to 237 < P(days) < 860, 0.5 < R(R_⊕) < 1.5. The gray box represents the region within 20% of Earth's orbital period and radius, considered by Burke et al. (2015). The gold box represents the region 50 < P(days) < 200, 1 < R(R_⊕) < 2, also considered by Burke et al. (2015).
Download figure:
Standard image High-resolution image

**Figure 13.** A depiction of the approximations to η_⊕ considered in this work with respect to Earth's flux from the Sun and its radius. All confirmed planets with known radii (green squares) are shown to provide a sense of scale. Regions used to compute estimates of η_⊕ are also shown. The "×" marks ${{\rm{\Gamma }}}_{\oplus }=\tfrac{{\partial }^{2}f}{\partial \mathrm{log}P\partial \mathrm{log}R}$ , evaluated at R = R_⊕, P = 365 days. The teal box represents the habitable zone adopted by Hsu et al. (2018), equivalent to 237 < P(days) < 500, 1.0 < R(R_⊕) < 1.75. The black box represents the habitable zone adopted by Zink & Hansen (2019), equivalent to 222.65 < P(days) < 808.84, 0.72 < R(R_⊕) < 1.7. The maroon box represents the SAG13 definition of η_⊕, equivalent to 237 < P(days) < 860, 0.5 < R(R_⊕) < 1.5. The gray box represents the region within 20% of Earth's orbital period and radius, considered by Burke et al. (2015). The gold box represents the region 50 < P(days) < 200, 1 < R(R_⊕) < 2, also considered by Burke et al. (2015).
Download figure:
Standard image High-resolution image

Table 5. Updated Occurrence Rates for Planets Orbiting GK Stars

Statistic	Position in Parameter Space ^a	Initial Estimate ^b	Updated Estimate	% Difference
Γ_⊕	$\tfrac{{\partial }^{2}f}{\partial \mathrm{log}\ P\partial \mathrm{log}\ R}{\| }_{R=1,P=365}$	${0.11}_{-0.05}^{+0.08}$	${0.12}_{-0.05}^{+0.08}$	+6
F₁	50 < P < 200, 1 < R < 2	${0.16}_{-0.03}^{+0.03}$	${0.17}_{-0.03}^{+0.04}$	+5
ζ_⊕	292 < P < 438, 0.8 < R < 1.2	${0.018}_{-0.008}^{+0.012}$	${0.019}_{-0.008}^{+0.013}$	+6
SAG13 HZ	237 < P < 860, 0.5 < R < 1.5	${0.15}_{-0.07}^{+0.11}$	${0.16}_{-0.07}^{+0.12}$	+6
Hsu et al. (2018) HZ	237 < P < 500, 1.0 < R < 1.75	${0.078}_{-0.03}^{+0.03}$	${0.082}_{-0.026}^{+0.036}$	+5
Zink & Hansen (2019) HZ	222.65 < P < 808.84, 0.72 < R < 1.7	${0.17}_{-0.06}^{+0.10}$	${0.18}_{-0.07}^{+0.10}$	+6
Super-Earth	0.355 < P < 384.843, 1.0 < R < 1.75	${3.23}_{-0.17}^{+0.18}$	${4.06}_{-0.21}^{+0.23}$	+26
Super-Earth^c	0.5 < P < 50, 1.0 < R < 1.75	${1.92}_{-0.10}^{+0.10}$	${2.40}_{-0.12}^{+0.12}$	+25
Sub-Neptune	0.928 < P < 509.997, 1.75 < R < 3.5	${1.95}_{-0.08}^{+0.08}$	${2.46}_{-0.10}^{+0.10}$	+26
Sub-Neptune^d	0.5 < P < 50, 1.75 < R < 3.5	${0.94}_{-0.04}^{+0.04}$	${1.18}_{-0.05}^{+0.05}$	+26

Notes.

^aRadii are expressed in R_⊕, and periods are expressed in days.^bWe calculate initial and final estimates of each statistic with the same pipeline, modified from Bryson et al. (2020; see Section 5).^cWe consider super-Earths once more, now within a tighter period range.^dWe consider sub-Neptunes once more, now within a tighter period range.

Download table as: ASCII Typeset image

5.3. Impact on the Frequency of Super-Earths

We next repeat our steps above for a second planetary population: super-Earths. In making a distinction between super-Earths and sub-Neptunes, we refer to Teske et al. (2018), who sought to understand how both detected and undetected stellar companions affect estimates of the exoplanet radius distribution—with particular focus paid to the "gap" in the observed radius distribution between super-Earths and sub-Neptunes (Fulton et al. 2017). Per Teske et al. (2018), we here define super-Earths as falling within the radius range of 1.0R_⊕ < R_p < 1.75R_⊕. For our period range, we first adopt the extremes of confirmed super-Earths discovered from Kepler data: Kepler-845b, with a period of 0.355 days (Dai et al. 2019), and Kepler-452b, with a period of 384.843 days (Morton et al. 2016). Again restricting ourselves to GK stars, we see an increase in super-Earth occurrence from ${3.23}_{-0.17}^{+0.18}$ to ${4.06}_{-0.21}^{+0.23}$ —an increase of 26% (4.5σ). Tightening the period range to 0.5 < P(days) < 50 yields a similar increase in occurrence, from ${1.91}_{-0.10}^{+0.10}$ to ${2.40}_{-0.12}^{+0.12}$ (+25%, 4.6σ). As above, both of these estimates are calculated by the same pipeline within this paper. Unlike our analysis of potentially habitable planets, we here make no cuts informed by planetary insolation flux.

5.4. Impact on the Frequency of Sub-Neptunes

Finally, we perform our occurrence rate comparison for the sub-Neptune class of exoplanets. Once more adopting size categories from Teske et al. (2018), we take this class of planets to occupy the parameter space region 1.75R_⊕ < R_p < 3.5R_⊕. We initially bound our periods with those of Kepler-845b, with a period of 0.928 days (Gajdoš et al. 2019), and Kepler-1630b, with a period of 509.997 days (Morton et al. 2016). We determine an increase in sub-Neptune occurrence from ${1.95}_{-0.08}^{+0.08}$ to ${2.46}_{-0.10}^{+0.10}$ , which corresponds to a 26% (6.4σ) increase. As with the super-Earth calculation, we find that constricting the period range to 0.5 < P(days) < 50 results in a similar change in occurrence, from ${0.94}_{-0.04}^{+0.04}$ to ${1.183}_{-0.046}^{+0.046}$ (+26%, 6.7σ). Once more, we perform these calculations without taking into account planetary insolation flux.

5.5. Comparison to Theory

Bouma et al. (2018) constructed models of increasing complexity to explore the effect of undetected stellar companions on exoplanet occurrence rates determined from transit survey data. With their most sophisticated model, they predicted that for small (R_p < 2R_⊕) planets, calculations ignoring stellar binarity will overestimate occurrence by up to 50%. This is in contrast to our study, which finds that transit surveys have been underestimating exoplanet occurrence by ignoring the effects of nearby stellar companions.

This discrepancy is reconciled by noting that Bouma et al. (2018) and our study account for the presence of undetected binaries in data in different manners. In our study, we remove some sets of stars that are no longer amenable to the detection of our planet of interest because of the presence of a nearby stellar companion. Bouma et al. (2018) did this in addition to accounting for the higher number of stars that have been searched if targets previously thought to represent one stellar component actually represent multiple stars. The former adjustment increases occurrence estimates, while the latter in general decreases them; these decreases, however, are sensitive to assumptions about planet occurrence in stellar binaries. For example, Bouma et al. (2018) found an increase in hot Jupiter occurrence when assuming that they were as likely to form around secondaries as primaries.

Accordingly, we may attribute the difference between the two methods to Bouma et al. (2018) including an additional effect. They noted that the suppressive effect of stellar binarity on exoplanet occurrence (from an observational, not formation-based, perspective) is also more pronounced for smaller planets; this may account for the occurrence increase in our study being less for smaller planets than for larger planets. Further studies will be required to better link the two approaches, namely by determining whether the process of considering more searched stars (i.e., stellar companions) in the calculation, as opposed to removing them altogether, can be incorporated into estimates of transit search completeness, candidate reliability, and planet occurrence. Furthermore, the Bouma et al. (2018) study considered planetary occurrence primarily as a function of planet radius, complicating direct comparisons to our work.

Finally, the discrepancy between our results and the simulations by Bouma et al. (2018) may be caused in part by biases in our modestly sized target sample. Upcoming analyses of a larger sample of Kepler target stars using AO observations at Palomar Observatory (J. Christiansen et al. 2020, in preparation) and speckle observations (K. Hardegree-Ullman et al. 2020, in preparation) will explore the effects of target selection on our conclusions by probing complementary regions of parameter space and expanding the target sample.

6. Conclusion

By observing 71 Kepler target stars using ShARCS on the Shane 3 m telescope, we detect 14 companions within 4'' of 13 Kepler objects. Our detection process, when compared to the companions revealed in Gaia DR2 data, is shown to be robust down to the predicted magnitude limits of ShARCS. We ascertain the effects of these companions on the completeness of Kepler's search for transiting planets and the resulting estimates of planet occurrence rates. First, we determine the physical association of our detected companions by making use of Gaia-derived distance estimates, proper motion comparisons, and color–magnitude analysis based on our ShARCS data (the canonical method for determining physical association of stars in AO imaging surveys). These methods yield a bound fraction of 21% within 4'' and 67% within 1''. In the process, we find that the canonical method produces results that are commensurate with analysis of data from Gaia DR2. With respect to our Gaia analysis, we find that our target stars with bound companions are not noted as binaries per the "binary = 1" or "binary = 3" flag in Berger et al. (2018), indicating the complementary approach that our study takes to identifying bound stellar pairs.

Utilizing our bound/unbound designations, we calculate reliable radius correction factors necessary to compute the radii of any transiting planets revealed in the Kepler light curves of our target stars. The mean radius correction factor if planets were to orbit their primaries is 1.10, while the mean radius correction factor if planets were to orbit their stellar companions is 1.78 (1.74 including those with stars with ambiguous physical association). This increase in radius affects the detectability of Earth-like planets—7.8% ± 0.2% of amenable Kepler stars would no longer be amenable if our radius correction distribution were sampled and applied to all such stars.

Finally, we utilize our detected companion rate to probabilistically filter out a fraction of Kepler stars without planet candidates in computing the occurrence rate of Earth-like planets around Sun-like stars. With the open-source Kepler DR25 completeness and reliability products, we find an average fractional increase of 6% across numerous estimates of η_⊕. This discrepancy is well within the error bars produced by our pipeline. This implies that assumption of stellar singularity does not greatly impact estimates of η_⊕ due to the large errors on current estimates. When we repeat our calculations for the super-Earth and sub-Neptune populations, we find increases in occurrence of 26% (a 4.5σ and 6.4σ increase, respectively). While these increases in occurrence when incorporating unresolved binaries run contrary to the decreases predicted by Bouma et al. (2018), this discrepancy may be due to differing treatments of stellar binarity on occurrence and the small size of the stellar sample considered in this pilot study. Our occurrence rate step of probabilistic filtering of stars based on observed multiplicity rates can be reproduced after factoring in our observed companion rate to refine an input stellar sample to occurrence rate calculations by other methods.

To improve the robustness of these estimates, future analysis should take into account the Poisson errors on our observed multiplicity rates and better address our observational biases. With respect to the latter point, for instance, while Figure 10 indicates that our target stars are distributed over a representative range of the amenable stars' temperature and log radius, our target stars were systematically brighter than the broader set of amenable Kepler stars (Figures 1, 11). The observations of this pilot study will also be complemented by ongoing AO observations at Palomar Observatory (J. Christiansen et al. 2020, in preparation) and speckle observations (K. Hardegree-Ullman et al. 2020, in preparation) of a larger sample of Kepler targets without detected planet candidates.

Future, comprehensive metastudies will be able to combine estimates of stellar multiplicity in the Kepler field—ultimately painting a fuller picture of how the components of stellar systems impacts exoplanetary occurrence.

This research has made use of the NASA Exoplanet Archive and ExoFOP, which are operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

A.W.M. is supported by the NSF Graduate Research Fellowship grant no. DGE 1752814. We acknowledge funding support from the Hellman Family Faculty Fund, the Sloan Foundation, and the David and Lucile Packard Foundation.

We thank Ellianna S. Abrahams for helpful conversations and insights. We also thank the anonymous reviewer for their thoughtful and detailed comments.

Facilities: ADS - , Exoplanet Archive - , Gaia - , Shane (ShARCS infrared camera). -

Software: astropy (Price-Whelan et al. 2018), emcee (Foreman-Mackey et al. 2013), IPython (Pérez & Granger 2007), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), pandas (McKinney et al. 2010), photutils, (Bradley et al. 2019), SciPy (Virtanen et al. 2020), tqdm (da Costa-Luis 2019).

A Closer Look at Exoplanet Occurrence Rates: Considering the Multiplicity of Stars without Detected Planets

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Target Sample