#  paired: A Statistical Framework for Determining Stellar Binarity with Gaia
RVs. I. Planet Hosting Binaries

Quadry Chance Department of Astronomy, University of Florida, 211 Bryant Space
Science Center, Gainesville, FL, 32611, USA  Center for Computational
Astrophysics, Flatiron Institute, Simons Foundation, 162 5th Ave, New York, NY
11238  Daniel Foreman-Mackey Center for Computational Astrophysics, Flatiron
Institute, Simons Foundation, 162 5th Ave, New York, NY 11238  Sarah Ballard
Department of Astronomy, University of Florida, 211 Bryant Space Science
Center, Gainesville, FL, 32611, USA  Andrew R. Casey School of Physics &
Astronomy, Monash University, Victoria, Australia  Center of Excellence for
Astrophysics in Three Dimensions (ASTRO-3D), Australia  Trevor J. David Center
for Computational Astrophysics, Flatiron Institute, Simons Foundation, 162 5th
Ave, New York, NY 11238  Adrian M. Price-Whelan  Center for Computational
Astrophysics, Flatiron Institute, Simons Foundation, 162 5th Ave, New York, NY
11238

###### Abstract

The effect of stellar multiplicity on the formation and evolution of planetary
systems is complex. At a demographic level, campaigns with both high-
resolution imaging and radial velocity observations indicate that planet
formation is strongly disrupted by close binaries, while being relatively
unaffected by wide companions. However, the magnitude and distance-limited
nature of those tools mean that large ranges of mass ratios and separations
remain largely unexplored. The Early Data Release 3 (EDR3) from the Gaia
Mission includes radial velocity measurements of over 6.5 million targets,
which we employ to explore the effect of binary companions within a
statistical framework called `paired`. These companions present as a source of
excess radial velocity noise in the Gaia catalog, when compared to the typical
noise for stars of similar spectral type and magnitude. Within this framework,
we examine the evidence for stellar multiplicity among the stars surveyed by
NASA’s Kepler and TESS missions. We use radial velocity errors published in
Gaia EDR3 to estimate the probability of an unresolved stellar companion for a
large subset of the Kepler and TESS Input Catalog stars, where possible
benchmarking our inferred radial velocity semi-amplitudes against the those
from ground-based radial velocity surveys. We determine that we are typically
sensitive out to several AU and mass ratios $>0.1$, dependent upon the stellar
magnitude. We aim for `paired` to be a useful community tool for the
exploration of the effects of binarity on planets at a population level, and
for efficient identification of false-positive transit candidates.

Binary stars (154), Close binary stars (254), Spectroscopic binary stars
(1557), Planet hosting stars (1242), Exoplanet astronomy (486),
Astrostatistics tools (1887), Astrostatistics (1882)

††journal: ApJ††software: astropy (Astropy Collaboration et al., 2013, 2018),
numpy (Harris et al., 2020), matplotlib (Hunter, 2007), pandas (pandas
development team, 2020; Wes McKinney, 2010), numpyro (Phan et al., 2019)

##  1 Introduction

The rapid growth in the number of known exoplanets has enabled subsequent
detailed studies of planetary demographics. Investigations of the key
relationship between stellar properties and planet occurrence emerged within a
decade of the discovery of the first planet around a main-sequence star (e.g.
Fischer & Valenti 2005). The addition of thousands more planets to the known
sample has allowed exoplanetary astronomers to map this relationship with
increasing detail. The knowledge that planets exist in binaries followed only
a year after the discovery that planets reside around other main-sequence
stars: the fourth detected exoplanet orbits the K dwarf member of the binary
pair 55 Cnc A and B. Twenty years after the first detection of a planet
orbiting a single member of a binary, the discovery of Kepler-16b showed that
some planets even orbit both members. Circumbinary planets and planets in
close S-type systems alone now number in the dozens. About half of the FGK
stars in the Milky Way have main-sequence companions (Offner et al., 2022),
making this question a critical one for the theory of planet formation.

Previous studies with both adaptive optics imaging (AO) and radial velocity
(RV) campaigns have established that binarity plays an important role in
shaping planetary systems e.g. (Holman & Wiegert, 1999; Paardekooper et al.,
2008; Marzari & Gallina, 2016; Kraus et al., 2012; Su et al., 2021). Wang et
al. (2014) showed that planet occurrence drops to a third for stars with
companions at separations of 100 AU, when compared to single stars. This
figure drops to a tenth for stars with companions at 10 AU. Indeed, Moe &
Kratter (2021) showed that all but 15% of planetary systems are suppressed due
to the presence of companion stars interior to 10 AU. Ziegler et al. (2021)
found that suppression at the 15% survival level holds for separations as wide
as 50 AU. With an adaptive optics imagining campaign sensitive to larger
separations, Kraus et al. (2016) found that 30% of planetary systems survive,
when the host star has a companion interior to 50 AU (relative to stars
without a companion). This suppression effect, per Moe & Kratter (2021) seems
to weaken or disappear at “wide” separations of $>$ 200 AU.

The physical mechanism responsible for the effect of binarity on the formation
of planets may occur at the disk stage, in some cases entirely preventing that
disk from forming (Silsbee & Rafikov, 2015). The literature on the impact of
stellar multiplicity upon planets often employs an “all-or-nothing” approach,
with a binary companion either “suppressing” planet formation altogether or
allowing it to proceed completely undisturbed. However, works such as Martin &
Triaud (2014) and Martin & Lubow (2017) have shown that if the disk survives,
a binary companion may affect the orientation and characteristics of that
protoplanetary disk. Such an effect could manifest as different planet
demographics around stars with binaries at intermediate separations: far
enough to allow the disk to survive but close enough to modify the disk or
resulting planets. Planet suppression changes steeply in the 1-50 AU
separation range, so that it stands to reason that the planetary systems that
do reside in these types of binaries illuminate the conditions under which
planets can survive.

At present, incompleteness in our understanding of the rate at which binary
systems host planets, as a function of their separation, hinders detailed
study of planets detected by large transit surveys. ESA’s Gaia mission, with
both astrometric and radial velocity observations, presents a novel
opportunity to study binarity on a large scale. For example, El-Badry et al.
(2021) identified 1.3 million spatially-resolved binaries from the Gaia
catalog at the $>$90% confidence level. This type of analysis is timely for
exoplanet science, as the increasing number of planets places increasing
pressure on traditional time- and resource-intensive campaigns with adaptive
optics (Ziegler et al., 2018, 2021) or radial velocity follow-up (Teske et
al., 2021; Chontos et al., 2022). Out of necessity, observations to identify
binary companions often focus upon known planet host stars, to eliminate
false-positives or identify the effects of stellar contamination on inferred
planet radius (Ciardi et al., 2015). The Sisyphean task of observing even a
fraction of the stars in a transit survey like Kepler often means that
campaigns must winnow this planet host star sample down yet further. To
maximize scientific return, there is a strong incentive to focus resources
upon the most scientifically interesting systems, such as those with multiple
planets. For the study of exoplanets, Gaia presents a twofold opportunity to
extract binarity information: first, it observes planet host stars in uniform
way. Searches for stellar companions by ground-based means are shaped
individually by survey completeness, which is itself informed by spectral
resolution, stellar magnitude, and observing cadence. This naturally renders
some portions of companion parameter space (mass ratio and orbital separation)
much better studied for some stars than others, for which comparatively little
is known. Secondly, Gaia observes non-planet-hosts in the same manner (a
sample at least an order of magnitude larger still for transit surveys),
enabling a precious field star comparison. Such comparisons are ultimately
crucial for discerning the effects of stellar multiplicity upon planet
occurrence at a demographic level.

To this end, we describe in this manuscript a technique to leverage available
Gaia radial velocity data to investigate stellar binarity. As a complement to
survey work with ground-based imaging and radial-velocity surveys, we have
devised a probabilistic model to extract binary likelihood for every star with
radial velocity observations from Gaia. Presently, only two quantities from
these radial velocity observations are publicly available. These are (1) the
number of observations and (2) the radial velocity “error”. The latter value
is derived under the assumption of a constant radial velocity model for the
star, and thus will (we argue in this manuscript) increase when that
assumption fails: namely, when the stellar velocity is not constant, but
rather varying due to a companion. As yet, Gaia has not released the
individual radial velocity time series measurements for most targets, and our
understanding of binarity will ultimately be refined with future public data
releases. Yet, we demonstrate here that binarity information is encoded in
these two (at present) publicly available quantities, at a sufficient level
for probabilistic measurements of binarity. On an individual stellar level, a
probability extracted from two reported RV quantities has perhaps limited
utility, but the same probabilities for every member of a large-scale survey
of exoplanet host stars can be illuminating. While Gaia has gathered radial
velocities for more than 7 million stars, in this paper we focus specifically
upon the subsample of these stars observed by NASA’s Kepler and TESS missions.
A total of 29,698 stars in the Kepler Input Catalog (about 15% of the total
stars observed) possess Gaia radial velocity data, along with 854,652 stars in
the TESS Input Catalog (about 9% of the total). This method is complimentary
to astrometric binary detection methods using RUWE (Wood et al., 2021; Penoyre
et al., 2022a, b; Andrew et al., 2022). Radial velocity methods are, by
nature, more sensitive to close companions and equal mass binary systems. This
manuscript is organized in the following manner. In Section 2 we describe our
machinery for extracting the likelihood of stellar binarity from the publicly
available radial velocity Gaia data. We assess our sensitivity to stellar
companions as a function of mass ratio and orbital separation with a series of
injection-and-recovery tests. In Section 3 we apply this machinery to the
sample of Kepler and TESS Input Catalog stars with Gaia radial velocity data.
We go on to validate our findings against important benchmarks. More
specifically, we investigate how our sample of “likely binaries” overlaps with
samples of binaries identified by other means. These means include the
presentation of eclipses in photometry, the location on the HR diagram, and
excess radial velocity noise from ground-based radial velocity surveys. In
Section 3 we describe findings from the population of computed binarity
probabilities in our sample. We touch upon a number of future scientific
applications for these probabilities, including:

  * •

Vetting of planet candidates and identification of false positives

  * •

Comparison between distributions of planet- and non-planet-hosting stars

  * •

Identification of circumbinary planet candidates

  * •

Direct detection of sub-stellar-mass companions in the brown dwarf and giant
planet mass range, and

  * •

Examination of confirmed planetary systems with anomalous radial velocities.

This list is not complete, and we are hopeful that the catalog will be broadly
useful to the community. While a detailed consideration of the effects of
stellar multiplicity upon planet occurrence is out of the scope of this paper,
we hope that our findings will be applicable to that end. We summarize and
conclude in Section 5.

##  2 Methods

The connection between radial velocity uncertainty and stellar multiplicity is
well-known (Maxted & Jeffries, 2005; Clark et al., 2011; Maoz et al., 2012).
As RV surveys have evolved to observe hundreds of thousands to millions of
stars within the last two decades, the method of using noisy, sparse radial
velocity data to statistically characterize stellar multiplicity has evolved
alongside them (Price-Whelan et al., 2017, 2020). The basic method begins with
a simple assertion: stars have some intrinsic distribution of radial velocity
noise that is created by internal processes. If this is true, single stars
should form a Gaussian distribution of RV noise centered on a typical value.
In the case where the source is an unresolved binary, semi-amplitudes for
binary stars are generally in the 1-100 km/s range. Most internal processes do
not create an RV signal on that scale. If the criteria for identification as a
binary is set conservatively at 3-$\sigma$, few single stars should be. The
details of this operation are expanded upon in Appendix A; in short, we have
reduced the question of stellar binarity to this: does a source have RV noise
in excess of the amount expected for similar-looking single stars?

The publicly available (as of EDR3) measurements for targets of the Gaia
radial velocity survey comprise three values: (1) the median radial velocity,
(2) the uncertainty on this median radial velocity, and (3) the number of Gaia
radial velocity observations (called “transits”) used to quantify (1) and (2).
The first two values are derived under the assumption that the measured radial
velocities are attributable to a single star. Our procedure for extracting
information about stellar binarity hinges upon a key assumption: that a radial
velocity model for a single star will furnish a poor fit to a bona fide binary
radial velocity curve. In the simplest case, this problem is akin to
characterizing a sine wave with a single “median” value. The estimated
uncertainty on the median value will correspondingly be inflated, since in
this case the variability due to the underlying sine wave is
(mis)characterized as “noise.” The degree to which a single star model will
furnish a particularly “poor” fit, rendering the binary identifiable, is
dependent upon several factors. These include the time sampling of the radial
velocity curve and the properties of the binary itself, including mass ratio,
orbital separation, eccentricity, and alignment relative to the line of sight.
In this section, we describe our methodology for identifying likely binaries
via the reported radial velocity uncertainty. In Section 2.1 we detail our
procedure for identifying anomalously high radial velocity errors, which is
based upon a comparison among similar stars. In Section 2.2 we characterize
our sensitivity to binaries as a function of mass ratio, orbital separation,
and stellar magnitude.

###  2.1 Binary detection & characterization

Appendix A and Foreman-Mackey, Chance, et al. (in prep.) describe our method
of binary detection and characterization in more detail, but we also summarize
the key results here. As discussed above, the Gaia catalog includes an
estimate of the “RV error” for each RV target. Since this error is estimated
by computing the sample variance of the RV time series, this value encodes
information about both the per-transit RV measurement uncertainty and the
excess RV variability introduced by companions.

To identify probable binaries in this sample, we first build a data-driven
model to estimate the per-transit RV measurement uncertainty as a function of
color and apparent magnitude (see Appendix A for more details), and then use
this estimate to identify targets with statistically significant excess RV
error. The key quantity that we compute is a $p$-value for each RV target in
the sample under the null hypothesis of a single star time series. Informally,
for a given target, this $p$-value quantifies the probability that such an
extreme RV error would have been measured if the system were in fact a single
star system with measurement uncertainty. Somewhat more formally, for a target
$n$, this $p$-value can be written as

$p_{n}=\mathrm{Pr}(\sigma>\sigma_{n}\,|\,\mathrm{target}\,n\,\mathrm{is\,single})$
(1)

where $\sigma_{n}$ is the observed RV error for target $n$. It’s important to
note that this is not the same as quantifying the probability that a given
target is a binary but, as we show below, it is a useful quantity for
selecting candidates. In the following sections, we identify probable
candidates by selecting targets with $p$-values below $0.001$, and then
calibrate the sensitivity of this method using injection and recovery tests.

Beyond identifying likely binary candidates, we also place probabilistic
constraints on the RV semi-amplitude for each target, under the assumption
that the RV variations can be described by Keplerian motion. The details for
this procedure are described in Appendix A, and in Section 3.3 we demonstrate
the accuracy and precision of these measurements when compared to
spectroscopic binaries detected by targeted RV surveys.

###  2.2 Injection and Recovery Analysis

We summarize here an injection-and-recovery analysis to investigate our
sensitivity to stellar binarity, within the framework we outline in the
previous section. To this end, we create a synthetic population of binaries,
which we “observe” at a typical Gaia cadence and with typical Gaia per-
measurement radial velocity uncertainty. We then calculate the median radial
velocity and radial velocity error of these “observations”, to obtain a
synthetic data set that resembles the reported radial velocity Gaia data.
Given a ground truth knowledge of the population, we can then assess the
ability of our machinery to correctly identify sources of anomalously high
radial velocity error. In Section 2.2.1, we describe the impact of Gaia sample
selection on this empirical sensitivity map.

First, we draw our underlying stellar binary properties with uniform log
spacing in semi-major axis, $a$, and in mass ratio, $q$. Our range in $q$
spans $10^{-4}$ to $1$, which for a Sun-like star brackets a range with a
Jupiter-mass exoplanet at one extreme and an equal mass stellar companion at
the other. Our semi-major axis range spans $10^{-2}$ \- $10^{3}$ AU, but our
sensitivity falls off substantially for $a>1$ AU. Within this $\\{q,a\\}$
parameter space, we generate one million binary star pairs. We then assign
orbital inclination relative to the line of sight, orbital eccentricity, and
longitude of periapse from uniform distributions of $\cos{0}-\cos{\pi/2}$,
0-1, and 0 to 360 degrees, respectively. Next, we calculate the binary’s
orbital period and employ Kepler’s equations to solve for the positions and
velocities of the two stars. The number of synthetic observations $n$ is drawn
from the distribution of RV transits shown in Figure 1. This depicts the
distribution in the number of transits for the observed Gaia stars in the
Kepler Input Catalog and the TESS Candidate Target List. We “measure” the
radial component of the velocity at uniform random times $t_{n}$ over the 668
day observing time span of EDR3. We also calculate the semi-amplitude of the
heavier (and presumably brighter) component.

Figure 1: Top: A histogram of the number of RV observations (transits) of
sources in the KIC crossmatch. Bottom: Identical plot for the TESS Community
Target List. The CTL sources tend to have more observations on average, along
with a longer tail toward more observations.

We add realistic per-measurement radial velocity error as follows. First, we
partition the real observed Gaia catalog stars into bins in color-magnitude
space. For each of these bins, we tabulate the average radial velocity error.
We randomly assigning each simulated stellar pair to a color-magnitude bin,
and draw our per-measurement error from the observed underlying distribution
for that bin. This procedure has an important caveat. In the idealized case,
the radial velocity noise distribution in each color-magnitude bin would
correspond to a population of single stars, with individual error budget
dominated by the measurement noise from Gaia’s spectrograph. However, in
practice each bin also contains bona fide binaries. In these cases, the
observed Gaia RVs may in fact be dominated by the reflex motion of the stars,
and associated “noise” may in fact be inflated. In this sense, our sample of
“control” stars is contaminated by some binaries. Other stars will have
intrinsically higher amounts of RV noise for reasons unrelated to a binary
companions. We seek to mitigate this contamination by using a very
conservation definition to detect “outliers”: at the $>$3-$\sigma{}$ level
from the mean in each bin. While this decision leaves us with a significant
false-negative rate, it also means that we can be reasonably certain that our
detections represent significant levels of noise.

Figure 2: Detection efficiency using RV errors as a function of injected
separation and mass ratio. The 99% recovery rate is shown for different mean
apparent magnitude contours. The color of each hexbin indicates the median
p-value for pairs in that bin.

We then run `paired` (as described in Appendix A) on this simulated sample of
binary stars, with our results depicted in Figure 2 as a function of semi-
major axis and mass ratio. The recovery rates with overplotted 99% recovery
rate contours are depicted by white contours corresponding to 5 different Gaia
magnitudes. We are most sensitive to binaries near equal mass, recovering
$>80\%$ for all magnitude bins. We find that our recovery rate drops of
steadily with mass ratio, but more quickly with separation.

The results of this injection-and-recovery test indicate that we have some
sensitivity to brown dwarfs around bright stars ($m_{g}<7$) out to a few AU
($logq$ of -1.1 for a Sun-like star). Around those same stars, we are
sensitive at the 99% confidence level to mass ratios of $q$=0.1 at 10 AU. It
is outside the scope of this work to compare raw binary occurrence to other
works, but we note (per Duquennoy & Mayor 1991) that 10 AU for a G dwarf
corresponds to an orbital period of $P\sim 10^{4}$ days, where 90% of
companions to G dwarfs have $q>0.1$. We therefore possess modest completeness
to stellar companions to bright G dwarfs. Specifically with respect to known
planet hosts, while the recovery rate drops off steeply with magnitude, we
retain the ability to detect stellar mass companions ($q>0.1$) interior to 1
AU around Sun-like stars down to $m_{g}$ of 12.6 with 99% confidence. While
Gaia only gathers radial velocity observations for a fraction of its targets,
we comment briefly on our empirical sensitivity as it applies to the
population of transit candidate host stars. This fiducial magnitude cutoff of
$m_{g}\geq$12.6, where we reliably recover stellar companions to Sun-like
stars interior to an AU, corresponds about half of the host stars to Kepler
Objects of Interest (Coughlin et al. 2016, hereafter KOIs). Roughly 80% of the
TESS Observations of Interest (TOIs) from the primary mission and 20% of TOIs
from the Faint Star catalog (Guerrero, 2021; Kunimoto et al., 2022) are
brighter than this magnitude threshold (here we have converted TESS and Gaia
magnitudes per Stassun et al. 2019). We conclude from our injection-and-
recovery analysis that, for the fraction of transit host stars with radial
velocity observations from Gaia (described further in Section 3), `paired` is
capable of detecting close stellar companions and occasionally substellar
companions.

####  2.2.1 Completeness to equal-mass, short-period binaries

There is an additional consideration for our injection-and-recovery analysis,
related to sensitivity in the portion of $\\{a,q\\}$ space corresponding to
short period binaries with mass ratios $q$ near to 1. While we would a priori
expect our sensitivity to peak in this corner of binary parameter space, we
note that sources with radial velocity errors $>$20 km/s are culled from the
Gaia sample of radial velocity stars and not reported. This down-selection
necessitates a careful consideration of the `paired` sensitivity map when
comparing to known samples of binaries, in the following way. Given a
hypothetical sample of 1000 Kepler stars of equivalent brightness and an
underlying rate of binarity of 50%, we might expect 10% of this sample to be
Gaia radial velocity targets. These would be selected randomly, so that we
could employ `paired` to interpret the Gaia RV data of 100 stars out of the
original 1000. In the simplified case of uniform sensitivity to all binaries,
we might expect `paired` to report that approximately 50 stars of the 100 are
“likely” binaries, reflecting the true underlying rate. However, in the case
of down-selection of the stars with highest apparent RV noise, perhaps the 10
most obvious binaries of the 50 are culled from the table of the stars with
radial velocity data from Gaia. Among the 90 in the subsequent public sample,
`paired` would detect the remaining 40 binaries, but the interpretation of a
40/90 (45%) underlying binarity rate would be incorrect. We anticipate that a
comparison of `paired` against known binary catalogs will always underestimate
the true rate for this reason. The degree of that underestimation depends upon
the fraction of binaries “obvious” enough to be culled from the catalog.

##  3 Results

In the previous section, we assess our theoretical sensitivity to binaries
with the machinery of `paired`. In this section, we describe our findings from
running `paired` on the 9% of the Kepler Input Catalog (29,698 stars) and 15%
of the TESS Input Catalog (854,652 stars) with EDR3 radial velocity
observations. Below, we compare our sample of “likely” binaries from this
exercise against existing benchmark indicators of binarity, to assess the
extent to which the outliers in radial velocity noise we are detecting are
bona fide stars in binary systems. We first consider `paired`’s performance on
the highest fidelity sample of individually identified eclipsing and
spectroscopic binaries (Section 3.1). We then turn to the results of large
stellar surveys, where indicators of binarity include position on the HR
diagram (Section 3.2) and radial velocity semi-amplitude (Section 3.3).

###  3.1 Validation Against Catalogs of Known Binaries

We construct a “likely binary” catalog by selecting all cross-matched KIC and
CTL targets with a $p$-value (described in Section 2) below $0.001$ (6,630 and
329,712 sources, respectively). We emphasize that the $p$-value is not the
same quantity as the false positive rate of binarity, but rather the
likelihood that the source’s radial velocity noise is indeed anomalously high
compared to stars of similar color and magnitude. A calibration of the latter
quantity with the rate of “false positives” would depend on the true fraction
of both binaries and astrophysical sources of false positives, which is beyond
the scope of this work. Figure 3 shows the cumulative $p$-value distributions
for all of our cross-matched sources. We find that approximately 22% of the
KIC and 38% of CTL stars exhibit excessive radial velocity noise at the
$p<0.001$ level.

Known eclipsing binaries (EBs) in the Kepler and TESS catalogs present an
ideal sample for testing `paired`’s performance. They comprise the types of
edge-on and short period systems to which we claim high sensitivity (per
Section 2.2), and can be confirmed as binaries independent of spectroscopic
measurements. We expect our rate of recovery to be similar when comparing to
the Kepler and TESS eclipsing binary catalogs. This is because, While TESS
stars are brighter on average than Kepler stars, eclipsing binaries reside in
the part of $\\{a,q\\}$ space where we retain high recovery rates even for dim
($m_{g}>15$) stars. Our recovery rate will depend on a combination of factors:
the distribution in $\\{a,q\\}$ space of the detected eclipsing binaries
(itself a function of the underlying transit survey sensitivity), our
detection efficiency over that range, and the magnitude distribution of the
stars.

Indeed, we find that for the 331 EBs in the Kepler EB catalog crossmatch
(containing 2922 sources) that we have at least 3 RV transits for (Kirk et
al., 2016), $73.8\%$ exceed the $p<0.001$ threshold. For the 458 EBs in the
TESS EB catalog crossmatch (containing 4584 sources) (Prša et al., 2022) we
recover $73.1\%$ as binaries.

Figure 3: Left: CDF of p-values of KIC crossmatch and Villanova Kepler
Eclipsing binary catalog (Kirk et al., 2016). The red dashed line is our
$p<0.001$ threshold. The intersection of this line with the CDFs is the
proportion of the sample that are flagged as binaries. Right: Identical plot
for the TESS Community Target List and Villanova TESS Eclipsing binary catalog
(Prša et al., 2022). Eclipsing binaries from both catalogs have excess RV
noise that meets our threshold in $\sim 70\%$ of stars compared to the ”field”
rate of $\sim 20\%$.

Our recovery fraction is modestly higher for spectroscopic binaries. As a
point of comparison, we perform a cross-match with the SB9 catalog (Pourbaix
et al., 2009). Of the 1669 single-lined spectroscopic binaries with more than
3 RV transits included in (Pourbaix et al., 2009), we recover $82.8\%$. We
estimate that this higher recovery rate is attributable to the similarities
between our sensitivity curve (reliant upon radial velocities) and a radial
velocity survey, as compared to the relationship between our sensitivity curve
and that of a transit survey.

###  3.2 Location on the Hertzsprung Russell Diagram

Figure 4 shows the median p-value of targets in bins on the color-magnitude
diagram for our KIC and CTL crossmatches. The yellow bins indicate an area
that (as a rough estimation) has a high fraction of binaries. From this
diagram, it can be seen that the expected features of an equal-mass binary
main sequence and increasing binary fraction with mass are present in both
catalogs.

Figure 4: _Top_ : A color-magnitude diagram of our cross match with the KIC
covering 29698 stars. _Bottom_ :A color-magnitude diagram of our cross match
with the CTL covering 854652 stars.The color corresponds to the median
$p$-value for that hexbin, with yellow (zero) indicating a high fraction of
binaries. We can see both the binary main sequence and the increasing binarity
with decreasing color.

###  3.3 Comparison to radial velocity catalogs

For roughly 35% of sources, we calculate the probability distribution for the
radial velocity semi-amplitude as described in Appendix A. Inferred semi-
amplitude posteriors from sparsely sampled data generally have long tails,
accommodating the possibility of larger eccentricity values (Price-Whelan et
al., 2017).

Predictably, agreement between our method and other RV surveys increases with
semi-amplitude. The vast majority of the crossmatched RVs are in agreement to
1-$\sigma$ (see Figure 5)(Price-Whelan et al. (2020), Price-Whelan, in prep.).
Our ability to infer a semi-amplitude degrades near the _Gaia_ RV noise floor
of about 1 km/s. Lower than that, the RV signal is indistinguishable from
noise. We also lose sensitivity at higher semi-amplitudes because these
objects are generally removed from the Gaia catalog because of the large RV
error.

Figure 5: Here we compare with the ground based mid-resolution RV surveys
APOGEE (top) and SB9 (bottom) for all the sources where we have overlap with
our catalog and can constrain semi-amplitudes. There is tight agreement for
high semi-amplitudes that stars to loosen as K approaches Gaia’s noise floor
of $\sim 1$ km/s.

##  4 Analysis

###  4.1 Exoplanet Candidate Vetting and Identification of False Positives

One of the anticipated uses of our catalog is in vetting exoplanet candidates.
Currently, one of the procedures for transforming a transiting planet
“candidate” to a confirmed planet is a follow-up RV campaign. Where possible,
a measured planetary mass is ideal, but such campaigns can also simply confirm
(by a lack of sizeable RV variations) that the transit signal is not being
caused by a grazing eclipsing binary.

The ability to pre-sort planet candidates by p-value into likely false
positives and bona fide candidates would enable follow-up resources to
prioritize the most promising candidates without taking a single RV
measurement. There also exists a retroactive application. Candidates which
were statistically validated my warrant a revisit if they have low p-values
that indicates a high-mass companion.

###  4.2 Potential sources of false positives

We emphasize here that a p-value calculation cannot be interpreted as the
false positive rate for our sample; these quantities are distinct and should
not be conflated. A known source of false positives is rapidly-rotating stars.
The rotational-broadened spectral lines of stars with $v\sin{i}>10$ km/s can
lead to an over-estimation of the RV measurement error from a source. There is
no way to distinguish this from RV noise from a companion, therefor rapid
rotators are likely to be flagged as binaries. A careful reexamination of the
assumption that $\sigma_{n}$ is constant at a given color and magnitude my
disentangle this noise from that more likely to be cause by companions, but
that is out side of the scope of this paper.

Figure 6: Top: CDF of p-values of the KIC crossmatch and KOI catalog
(exclusive of known false positives). The red dashed line is the threshold for
determining a star is a binary. The intersection of this line with the CDFs is
the proportion of the sample that are flagged as binaries. The most striking
feature of this plot is the reduction in the proportion of low p-value systems
hosting planets can be seen starting near the cutoff value of 0.001. Bottom:
Identical plot for the TESS Community Target List and TOI catalog (Guerrero,
2021). The TOI catalog (with known false positives removed) appears to be
substantially more contaminated with false positives than the KOI catalog.

####  4.2.1 Differences between the Kepler and TESS catalog

The cross match with Kepler contains 29,698 matches out of nearly 200,000
Kepler target stars, down to a Gaia g magnitude of 14.84. The relatively low
fraction is a result of removing stars without a minimum of 3 RV transits. The
cross match with TESS contains 854,652 out of approximately 9.5 million stars
in the Candidate Target List, down to a Gaia g magnitude of 15.81. The sources
in the KIC crossmatch on average have fewer RV transits that those in the CTL
list. Recall that Equation A1, $\epsilon_{n}$ is proportional to this number.
In essence, while the median number of transit is is similar, some of the CTL
sources are far more well sampled than the KIC sources. The Kepler Object of
Interest catalog has been repeatedly and thoroughly (statistically, in most
cases) vetted for false positives. The TESS Object of Interest catalog is
still being added to with follow-up observations and vetting on going.
Combined with a more lax policy for disposition as a planet candidate
(Guerrero, 2021), the TOI planet candidates are much more likely to be an
astrophysical false positive. The side-by-side comparison of the two catalogs
in Figure 6 is revealing in this aspect. There is a clear difference between
the p-values of planet host and non-hosts for Kepler, but the populations look
quite similar for TESS.

###  4.3 Binary probabilities of planet host stars and field stars

Calculating occurrence rates is a non-trivial problem, especially when using a
collection of surveys all with varying sensitivity constraints. What we will
present in this paper serves as more of a crude guide to what we should expect
to see in a full treatment. Although our results are only a statistic and not
the actual binary fraction of stars, it should be close enough for quiet main-
sequence stars. Previous studies of close binaries hosting exoplanets (Wang et
al., 2014; Moe & Kratter, 2021) have established that the presence of a close
binary is deleterious to planet formation. Particularly in the case of systems
the have stellar companions within 10 AU, the fraction of planet hosts
decreases precipitously. With that information in hand, we expect to see that
reproduced using our method of establishing stellar binarity. In the left
panel of Figure 6 we present the cumulative distribution functions of our
crossmatch with the Kepler Input Catalog (KIC) and the Kepler Objects of
Interest (KOIs). We have removed candidate planets and known false positives
in order to use the KOIs as our sample of planet hosts. Figure 6 shows the
distribution of the selected stars binary p-values where values closer to 0
indicate increasing likelihood of being a binary. These plots can also be read
as a kind of percentile ranking of the excess RV noise from each target. Our
relatively conservative limit of $10^{-3}$ or 99.9th percentile of RV noise
for binary classification should exclude the vast majority of targets that
have excess RV noise by chance. We note 2 important features. First, while the
two CDFs are quite similar approaching the cutoff as they descend in p-value,
after the cutoff their behavior diverges. The number of KOI targets classified
as binaries sharply decreases relative to the KIC targets. This decrease is
also insensitive to our choice of classification limit. We interpret this
feature as the expected decrease in the number of planet hosting systems in
close binaries. Second, there is a slope as the p-values approach zero.
Smaller p-values (and therefore larger signals) generally indicate higher mass
ratios and smaller separations. We therefore interpret this feature as systems
hosting larger and/or closer companions are less likely to host a KOI. These
findings are consistent with the observation that close binaries suppress
planet formation. This suppression is not, however, complete. Investigating
systems that appear to host a massive companion while also hosting planets has
the potential to more precisely determine what mechanism this suppression
operates by examining the surviving systems.

###  4.4 Direct detection of sub-stellar-mass compansions

Our machinery, when applied to the Gaia data, has marginal sensitivity to both
gas giants and brown dwarfs at short periods. These objects, while not meeting
the $p<0.001$ threshold we have set for stellar binaries do create signal at
the p $\sim 0.01-0.1$ level. This group of host stars will be 10 to 100x more
contaminated than our stellar binary even if our simple assumptions hold. Our
interpretation of this peak for sub-stellar objects at p= .1-.01 relies upon
our assumption for how p-values scale with our recovery rate as shown in
Figure 2. If recovery rate scales with p-value as shown, we can estimate, for
example, that $\sim 50\%$ of pval = 0.01 signals are generated by companions.

Figure 7: The distribution of p-values for bright massive hot Jupiter hosts.
Despite the low number of host stars in the sample (21), there is a clear
excess of 90th percentile-level RV noise, possibly from their previously known
planets.

To test this hypothesis, we used a population of known massive hot Jupiters.
For a selection of hot Jupiters with the largest semi-amplitudes, there is a
noticeable excess when compared to the rate for the entire sample Figure 7, as
we would expect if we are detecting sub-stellar objects. We also see a steep
decline in the fraction of lower p-values, which we would expect given that
hot Jupiters are unlikely to have close stellar companions creating clear
signals(Knutson et al., 2014; Piskorz et al., 2015; Ngo et al., 2015). The
juxtaposition of these two features implies that they are produced by sub-
stellar companions. Even if we cannot confidently call this a detection
individually from the Gaia RVs alone, we suspect that an RV search of these
targets would prove fruitful.

###  4.5 Identification of circumbinary planets

This catalog has the potential to uncover examples of one of the more elusive
planet formation outcome: misaligned circumbinary (or P-type) planets. Despite
the thousands of known exoplanets and indications that they should be just as
common as planets orbiting single stars (Martin & Triaud, 2014), only 12
transiting circumbinary planets are currently known (Kostov et al., 2020). The
fact that we only know of a small number of systems therefore must be
reflective of some combination of intrinsic and observational biases. For the
most part, we only know that the host stars are in a close binary because they
eclipse one another. Examinations of protoplanetary disks around young binary
stars indicates that $\sim 70\%$ of these disks are aligned with the binary
for for periods of less than 20 days, but the range widens all the way to
polar configurations for periods greater than 30 days (Czekala et al., 2019).
Circumbinary systems where this is not the case have to date not been
identified as such. For misaligned circumbinary systems, epochs where no
transit occurs are the majority (Chen & Kipping, 2021). The case of two or
more transits per epoch is one way to unambiguously identify a circumbinary
planet, but this only occurs $\sim 3\%$ of the time. These systems may make up
some portion of systems with single observed transits. Using this catalog is a
viable way to identify systems in this configuration. Conditioning a search
for these planets on evidence that a given planet host is a close binary
drastically reduces the amount of target light curves that need to be
searched.

Planets hosts that are flagged as binaries, especially those with inferred
semi-amplitudes indicating a close binary, represent the best candidates for
being a misaligned CBP. Worthy of particular scrutiny is KOIs meeting those
criteria that have only one or two planet transits present in their light
curves. The vast majority of CBPs that have the potential to transit do so
only in a “sparse” configuration (Martin & Triaud, 2014), where non-
consecutive transits are the norm. Many searches of the kepler catalog for
these single or sparsely transiting planet candidates have already been
conducted (Wang et al., 2015; Foreman-Mackey et al., 2016; Uehara et al.,
2016). The combination of this catalog and known single-transit planet
candidates may prove key in finding one of this as-yet unseen type of
planetary system.

Table 1: Kepler planet host candidates of particular interest with semi-amplitude ranges inferred from Gaia RVs. KOIs that meet the p-value threshold but have effective temperatures greater than 6250K are removed.  KIC | KOI | # of RV transits | $K_{16th}$ | $K_{50th}$ | $K_{84th}$ | KOI Disp. | Gaia m | $T_{eff}$
---|---|---|---|---|---|---|---|---
3440118 | K03876 | 9 | 11.828 | 15.255 | 25.527 | CP | 12.605 | 5620
4056616 | K02013 | 7 | 6.665 | 9.03 | 16.42 | CP | 13.197 | 4507
5450893 | K02970 | 9 | 13.244 | 17.522 | 28.727 | CP | 12.769 | 6046
6289257 | K00307 | 8 | 7.181 | 10.87 | 19.08 | CP | 12.757 | 6023
8240797 | K01809 | 6 | 11.258 | 15.67 | 29.812 | CP | 12.746 | 5772
8711794 | K00105 | 6 | 5.835 | 9.846 | 19.335 | CP | 12.853 | 5666
9489524 | K02029 | 11 | 6.118 | 8.562 | 13.435 | CP | 12.95 | 5208
3733628 | K00387 | 13 | 9.098 | 11.534 | 18.12 | CP | 13.364 | 4538
7742408 | K03478 | 8 | 11.557 | 15.131 | 26.124 | CP | 13.351 | 4928
5511081 | K01930 | 12 | 4.042 | 5.837 | 8.89 | CP | 12.076 | 5951
6021275 | K00284 | 6 | 5.157 | 8.466 | 16.42 | CP | 12.32 | 5931
7419318 | K00313 | 7 | 19.866 | 25.682 | 46.373 | CP | 12.931 | 5188
7584650 | K02631 | 5 | 26.003 | 37.206 | 77.673 | CP | 13.437 | 5871
9450647 | K00110 | 11 | 5.111 | 8.457 | 13.517 | CP | 12.584 | 6215
7749773 | K02848 | 7 | 4.169 | 7.087 | 12.983 | CP | 12.423 | 5691
8424002 | K03497 | 6 | 6.083 | 8.994 | 17.637 | CP | 13.317 | 3419
9886361 | K02732 | 12 | 6.152 | 9.117 | 14.103 | CP | 12.752 | 6170
10212441 | K02342 | 6 | 13.355 | 19.138 | 36.643 | CP | 13.041 | 5798
11086270 | K00124 | 12 | 4.873 | 8.545 | 13.538 | CP | 12.891 | 5869
11764462 | K01531 | 7 | 11.699 | 16.585 | 28.956 | CP | 13.122 | 5646
11760231 | K01841 | 10 | 7.892 | 11.024 | 18.273 | CP | 13.242 | 5187
9579641 | K00115 | 7 | 6.793 | 10.981 | 20.165 | CP | 12.794 | 5779
9881662 | K00327 | 7 | 13.075 | 18.197 | 31.932 | CP | 12.932 | 5945
10063208 | K04292 | 11 | 8.925 | 11.879 | 19.089 | CP | 12.833 | 5525
10528068 | K01162 | 6 | 6.11 | 10.292 | 20.193 | CP | 12.722 | 5751
11446443 | K00001 | 5 | 3.999 | 5.982 | 12.741 | CP | 11.283 | 5820
11017901 | K01800 | 10 | 4.361 | 6.428 | 10.578 | CP | 12.375 | 5494
11551692 | K01781 | 7 | 2.681 | 4.153 | 7.424 | CP | 12.226 | 4920
12301181 | K02059 | 7 | 23.879 | 30.858 | 58.983 | CP | 13.225 | 4997
12785320 | K00298 | 5 | 12.129 | 18.065 | 37.831 | CP | 13.208 | 5250
1025986 | K07621 | 6 | 3.363 | 4.57 | 9.213 | PC | 10.06 | 5604
1161345 | K00984 | 5 | 3.169 | 5.761 | 12.179 | PC | 12.391 | 5296
3129238 | K07645 | 11 | 3.51 | 5.273 | 8.401 | PC | 12.321 | 5341
3545478 | K00366 | 6 | 3.626 | 6.516 | 12.631 | PC | 11.556 | 6201
6309307 | K05262 | 9 | 2.981 | 4.175 | 6.824 | PC | 12.342 | 4973
4276716 | K01619 | 6 | 3.457 | 4.92 | 9.726 | PC | 11.801 | 4830
5802205 | K08107 | 8 | 1.808 | 3.095 | 5.392 | PC | 11.022 | 6180
7901016 | K04257 | 6 | 6.65 | 12.748 | 24.783 | PC | 13.355 | 5745
8702537 | K07905 | 7 | 6.238 | 10.191 | 18.693 | PC | 13.045 | 5616
10287723 | K01174aaSingle transit | 10 | 9.423 | 12.102 | 20.392 | PC | 13.428 | 4500
10352333 | K03991 | 8 | 12.876 | 17.388 | 29.455 | PC | 13.183 | 5519
10163067 | K07291 | 8 | 11.201 | 14.659 | 25.408 | PC | 12.574 | 6167
4139254 | K06108 | 10 | 35.969 | 44.797 | 67.957 | PC | 12.072 | 5571
4141593 | K07685 | 7 | 21.097 | 27.676 | 50.725 | PC | 13.251 | 5999
4450844 | K03168 | 7 | 3.167 | 4.268 | 7.642 | PC | 10.455 | 5843
5857656 | K05205 | 5 | 13.355 | 19.689 | 41.405 | PC | 12.771 | 5861
5858919 | K08109 | 11 | 2.086 | 2.688 | 4.307 | PC | 11.882 | 3297
8009811 | K05459 | 7 | 2.631 | 4.275 | 7.669 | PC | 11.419 | 6039
7661893 | K04296 | 8 | 31.578 | 39.644 | 67.002 | PC | 12.997 | 5615
8873090 | K02968 | 13 | 3.407 | 4.677 | 7.071 | PC | 11.422 | 6133
9330740 | K05656 | 17 | 6.516 | 8.645 | 12.801 | PC | 12.652 | 5717
9072639 | K04198 | 6 | 16.049 | 21.63 | 41.245 | PC | 12.55 | 6017
8005002 | K05457 | 8 | 3.298 | 5.352 | 9.242 | PC | 12.422 | 5146
9469494 | K07938 | 11 | 13.172 | 16.385 | 26.321 | PC | 11.655 | 5989
10736489 | K07368 | 9 | 8.596 | 11.287 | 19.17 | PC | 12.831 | 5227
11182260 | K07415 | 9 | 35.532 | 44.791 | 69.626 | PC | 13.247 | 4995
11182706 | K08046 | 12 | 3.158 | 4.58 | 7.123 | PC | 12.619 | 4859
11600299 | K05916 | 7 | 2.184 | 2.984 | 5.342 | PC | 12.052 | 3661
11502218 | K00970 | 17 | 35.866 | 41.432 | 56.48 | PC | 6.948 | 4813
12170767 | K05960 | 8 | 23.423 | 29.814 | 52.671 | PC | 13.3 | 5539

###  4.6 Confirmed planetary systems with anomalous radial velocities

Here we highlight a few systems that may warrant a closer look in light of
their known system characteristics and high binary probabilities.

####  4.6.1 Kepler-411

Kepler 411 has an inferred RV semi-amplitude of $4.15^{+3.27}_{-1.47}$ km/s
and a p-value of $2\times 10^{-4}$ placing it in the 100th percentile of
excess RV noise of similar looking stars. The system was characterized in Sun
et al. (2019). the outermost know transiting planet, Kepler-411d is seen to
have a strong interaction with with an outer, non-transiting perturber. This
investigation constrained the mass of planet e to $10.8\pm{1.1}\oplus$,
assuming that the system was co-planar. However, TTV inversion for exterior,
non-transiting perturbers is a fraught exercise. The system also sits well
above the single-star main sequence in a color-magnitude diagram. Given the
limited quality of the TTV data, it is difficult to eliminate the possibility
of a more massive inclined companions being responsible for planet d’s TTVs.

####  4.6.2 Kepler-132

Kepler-132A has an inferred semi-amplitude of $8.46^{+7.96}_{-3.31}$ km/s and
a p-value of $3\times 10^{-4}$ placing it in the 99.9996th percentile of
excess RV noise of similar looking stars. The system was briefly mentioned in
Lissauer et al. (2014) because of a curious characteristic: It a is a double
star of two nearly equal magnitude F7 components and two of the planets cannot
possibly exist while orbiting the same star. KOI-284.02 and 284.03 have
periods of 6.41 and 6.17 days, respectively. While EDR3 only has an RV
measurement for the slightly brighter component, recon spectra has indicated
that the stars had an instantaneous difference in radial velocity of about 1
km/s (Lissauer et al., 2014). The two stars have a projected separation of 450
AU. It is unlikely that the motion of these two stars over the EDR3 observing
baseline is causing the excess RV noise.

####  4.6.3 Kepler-410A

Kepler-410A has an inferred RV semi-amplitude of $1.72^{+2.57}_{-0.87}$ km/s.
This system is unusually bright (Gaia G magnitude of 9.38) so we can expect
the Gaia per-transit RV uncertainty to be around 100-200 m/s. This signal
cannot be explained by the known companion Kepler-410B, an M dwarf with a
period $\sim 2500$ years. This system has previously been investigated and
thought to possibly be host to another stellar component (Gajdoš et al.,
2017). Follow-up RV observation have ruled out a stellar mass companion.

####  4.6.4 KOI 1174

KOI 1174 has an inferred RV semi-amplitude of $12.10^{+8.29}_{-2.68}$ km/s and
a p-value less than $10^{-6}$. The system has a single observed transit with a
depth indicating a radius of $2.99R\oplus$. Single-transit objects are of
particular interest because they possibly represent detection of circumbinary
planets misaligned with their host stars. Such planets are not expected to
transit periodically, making them difficult to identify. Despite this, non-
eclipsing binaries make up the vast majority of binaries and inclined
circumbinary planets are dynamically allowable.

###  4.7 TOI 1897

TOI 1897 is an early G dwarf with an inferred RV semi-amplitude of
$6.27^{+3.39}_{-1.15}$ km/s and a p-value less than $10^{-6}$, placing it in
the 100th percentile of excess RV noise of similar looking stars. A cursory
examination of the SPOC light curve shows one transit too deep for a planetary
object and sever of a smaller object. The probability of this system being an
eclipsing binary where the primary transits have fallen into data gaps seems
quite high.

###  4.8 TOI 4573

TOI 4573 is a late K dwarf with an inferred RV semi-amplitude of
$5.39^{+1.90}_{-0.63}$ km/s and a p-value of less than $10^{-6}$, placing it
in the 100th percentile of excess RV noise of similar looking stars. The semi-
amplitude inference is unusually confident because this source has been
observed 78 times by Gaia RVs. This star is host to a $<4R_{\oplus}$ planet in
a to day orbit. Comments on the publicly available TRES follow-up spectra
indicate the there is evidence of an RV trend in the system that is out-of-
phase with the known transiting planet.

###  4.9 TOI 2119

TOI 2119 is an early dwarf with an inferred RV semi-amplitude of
$13.15^{+7.11}_{-1.98}$ km/s and a p-value of 0 placing it in the 100th
percentile of excess RV noise of similar looking stars. This star hosts a
known eccentric brown dwarf companion (Carmichael et al., 2022; Cañas et al.,
2022)

TIC | TOI | # of RV transits | $K_{16th}$ | $K_{50th}$ | $K_{84th}$ | TFOP Disposition | Gaia m | $T_{eff}$
---|---|---|---|---|---|---|---|---
285048486 | 1728 | 29 | 1.041 | 1.369 | 1.909 | CP | 11.726 | 3907.0
468777766 | 3750 | 17 | 5.224 | 7.015 | 10.397 | PC | 13.958 | 3815.0
230387153 | 2086 | 34 | 3.538 | 4.318 | 5.957 | PC | 11.822 | 5337.0
468889418 | 5156 | 18 | 5.4 | 6.401 | 9.079 | PC | 10.67 | 5857.9
29960110 | 1201 | 3 | 1.23 | 2.78 | 8.538 | CP | 12.089 | 3948.0
322262448 | 4043 | 11 | 9.727 | 12.585 | 20.412 | PC | 12.702 | 5926.0
418012030 | 2533 | 27 | 9.669 | 11.438 | 16.342 | PC | 12.002 | 6144.3
236934937 | 2291 | 18 | 4.145 | 5.111 | 7.284 | PC | 12.946 | 3869.0
180695581 | 1807 | 15 | 2.344 | 2.85 | 4.189 | CP | 9.684 | 4612.99
236387002 | 2119 | 15 | 11.176 | 13.154 | 20.264 | CP | 11.486 | 3606.0
198212955 | 1242 | 8 | 2.247 | 3.352 | 5.785 | PC | 12.348 | 4255.0
372758077 | 4168 | 5 | 24.065 | 33.144 | 66.019 | PC | 11.865 | 6016.0
347013211 | 2045 | 10 | 8.885 | 10.974 | 18.528 | PC | 11.215 | 6125.9
374348168 | 4282 | 11 | 11.363 | 14.219 | 22.697 | PC | 12.525 | 5780.0
168405330 | 2984 | 19 | 42.72 | 50.967 | 73.418 | PC | 12.558 | 4935.8
362103298 | 4573 | 78 | 4.764 | 5.392 | 7.295 | PC | 11.802 | 4207.0
232608943 | 4600aaSingle transit | 18 | 4.633 | 5.968 | 8.593 | PC | 12.434 | 5073.0
366499151 | 5115 | 6 | 10.096 | 13.762 | 26.46 | PC | 12.651 | 5129.0
32497972 | 876 | 16 | 2.023 | 2.712 | 3.999 | PC | 12.477 | 3955.0
318812447 | 2634aaSingle transit | 5 | 2.13 | 3.776 | 8.023 | PC | 11.505 | 5349.16
163539739 | 1278 | 16 | 4.116 | 5.07 | 7.317 | CP | 12.74 | 3841.0
123664207 | 2493aaSingle transit | 46 | 1.674 | 2.092 | 2.892 | PC | 12.181 | 4242.0
Table 2: TESS planet hosts of particular interest

##  5  Conclusion

In this manuscript, we describe a statistical framework called `paired` for
identifying stellar multiplicity using radial velocity information from Gaia
EDR3. The resulting p-values are extracted from the three presently available
measurements for a subset of Gaia target stars: the mean radial velocity, the
error associated with the mean radial velocity, and number of radial velocity
observations (or “transits”). From these values, we construct a statistically
robust mechanism for determining consistency between the reported radial
velocity error for a source, and the distribution of radial velocity error for
stars of similar color and magnitude. Where inconsistent, we generally assume
that the excess noise is indicative of additional orbital reflex motion of the
source (due to a companion). We demonstrate that this assumption is borne out
in general when comparing to other benchmarks of binarity: samples of known
eclipsing binaries, the clustering of “likely” binaries on the binary sequence
of the HR diagram, and existing radial velocity catalogs. We include in this
manuscript the resulting $p$-values reported by `paired` for the sample of
Gaia radial velocity targets observed by NASA’s Kepler and TESS missions. This
sample constitutes a a sizeable fraction (9% and 15%, respectively) of all the
stars in the Kepler Input Catalog and the TESS Community Target List.

There are, however, cases of “likely” binaries identified by `paired` that are
bona fide single stars. These cases highlight the ways in which our results
ought to be interpreted with caution. While meaningful in a statistical sense,
the nature of `paired`’s probabilistic diagnosis for an individual star merits
careful consideration. Some fraction of the time (per the definition of the
$p$-value), a star which we deem likely to exhibit excess RV noise does not do
so. It may also be the case that the excess apparent RV noise is due to a
physical mechanism other than reflex motion due to a companion. We focus on a
key example here that highlights the necessity for caution.

TrES-2 (Kepler-1) was the the first star identified in the Kepler field to
host a transiting planet, discovered from the ground several years prior to
launch (O’Donovan et al., 2006). Its early discovery and subsequent high-
precision photometry from Kepler render it a target of high interest, and it
possesses a decade and a half of radial velocity observations. This star is a
notable case for where archival RV data does not support the existence of the
companion with a 50th-percentile semi-amplitude of $\sim 5$ km/s that is
indicated in our catalog. Despite our advertised caution, this catalog is
intended to be useful wherever some idea of the close binary rate of stars is
useful. We have illustrated some possible uses for this catalog in
constraining the binary fraction of planet hosts and identifying known
planetary systems that have evidence for stellar or large sub-stellar
companions. We are particularly interested in the prospect of using this
catalog to particularize the suppression effect that close binaries have upon
planet occurrence. We can see the suppression effect in the crude exercise of
comparing the p-value distribution of known planet hosts to the field, but
occurrence rate calculations demand a better understanding of possible false
positive scenarios.

This work was performed in part through the Pre-Doctoral program at the
Flatiron Institute’s Center for Computational Astrophysics. Q. C.’s work was
supported by the National Aeronautics and Space Administration under Grant No.
80NSSC21K1841 issued through the Future Investigators in NASA Earth and Space
Science and Technology 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. This research has made use
of the Exoplanet Followup Observation Program website (ExoFOP), which is
operated by the California Institute of Technology, under contract with the
National Aeronautics and Space Administration under the Exoplanet Exploration
Program. This research has made use of the NASA Exoplanet Archive, which is
operated by the California Institute of Technology, under contract with the
National Aeronautics and Space Administration under the Exoplanet Exploration
Program. This work made use of the gaia-kepler.fun crossmatch database created
by Megan Bedell. We are thankful for helpful discussions with Natalia
Guerrero, Veselin Kostov, Kaitlin Kratter, Elisa Quintana, and Jamie Tayar.

##  Appendix A Mathematical details of the model

The next paper in this series (Foreman-Mackey, Chance, et al., in prep.)
provides a more detailed description of our methodology, but we summarize the
key elements here in order to provide some intuition for the procedure. For
each RV target, the Gaia EDR3 catalog reports the median RV, an estimate of
the error on this median, and the number of individual RV measurements (called
“transits”) that were used in this estimate. In our analysis, we ignore the
median RV and instead only consider the RV error and number of transits. The
RV error $\epsilon_{n}$ for a target $n$ is computed as (Katz et al., 2019)

$\displaystyle{\epsilon_{n}}^{2}=\left(\sqrt{\frac{\pi}{2\,T_{n}}}\,s_{n}\right)^{2}+0.11^{2}$
(A1)

where

$\displaystyle{s_{n}}^{2}=\frac{1}{T_{n}-1}\sum_{i=1}^{T_{n}}\left(v_{n,i}-\bar{v}_{n}\right)^{2}$
(A2)

is the sample variance of the individual RV measurements $v_{n,t}$. In
Equation A2, $t=1,\cdots,T_{n}$ indexes the $T_{n}$ transits for target $n$
and

$\displaystyle\bar{v}_{n}=\frac{1}{T_{n}}\sum_{i=1}^{T_{n}}v_{n,i}$ (A3)

is the mean RV. In all that follows, it is much simpler to reason about the
sample variance ${s_{n}}^{2}$ instead of the error on the median
${\epsilon}_{n}$, but that isn’t a problem because we can compute

$\displaystyle{s_{n}}^{2}=\frac{2\,T_{n}}{\pi}\left({\epsilon_{n}}^{2}-0.11^{2}\right)$
(A4)

for any target, using only values available in the catalog: $\epsilon_{n}$ and
$T_{n}$.

Now, for a given target $n$, in the absence of any excess radial velocity
signal, and if we know it’s per-transit radial velocity measurement
uncertainty $\sigma_{n}$, we would expect the quantity

$\displaystyle{\xi_{n}}^{2}$ $\displaystyle=$
$\displaystyle\frac{(T_{n}-1)\,s_{n}^{2}}{{\sigma_{n}}^{2}}$ (A5)

to be chi-squared distributed with $T_{n}-1$ degrees of freedom.

We can even go one step further, and derive the expected sampling distribution
for ${s_{n}}^{2}$ under some model for the excess RV signal. For example, in
our case, we want to infer the orbital properties of any binary systems based
on the observed ${s_{n}}^{2}$. In this case, our model for the expected RV
time series is

$\displaystyle\mu_{n}(t;\,\theta_{n})$ $\displaystyle=$ $\displaystyle
v_{0,n}+K_{n}\,\left[\cos(f_{n}(t)+\omega_{n})+e_{n}\,\cos\omega_{n}\right]\quad,$
(A6)

but the following derivation can be used more generally. We can evaluate
Equation A6, at observation times $\\{t_{n,i}\\}_{i=1}^{T_{n}}$ and orbital
parameters $\theta_{n}=\\{K_{n},\,e_{n},\,\omega_{n},\,\ldots\\}$. Of course,
in our case, we don’t actually _know_ the observation times111The Gaia
Observation Forecast Tool (https://gaia.esac.esa.int/gost/) predicts the
transit times for astrometric measurements, and we expect the RV measurements
will be a subset of these transits, but for most sources that information is
not enough to substantially improve our inferences., so in practice we will
need to marginalize those out, but for now, let’s proceed as if we do know
them. For a given set of parameters, it can be shown that the sampling
distribution for our scalar

$\displaystyle{\xi_{n}}^{2}$ $\displaystyle=$
$\displaystyle\frac{(T_{n}-1)\,s_{n}^{2}}{{\sigma_{n}}^{2}}$ (A7)

is a non-central $\chi^{2}$ distribution with $T_{n}$ degrees of freedom, and
non-centrality parameter

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{T_{n}}\left[\frac{\mu_{n}(t_{i};\,\theta_{n})-\bar{\mu}_{n}}{\sigma_{n}}\right]^{2}$
(A8)

where

$\displaystyle\bar{\mu}_{n}$ $\displaystyle=$
$\displaystyle\frac{1}{T_{n}}\sum_{i=1}^{T_{n}}\mu_{n}(t_{i};\,\theta_{n})$
(A9)

is the sample mean of the model evaluated at the observation times. This
sampling distribution defines a likelihood function for the parameters
$\theta_{n}$, and it is the key ingredient of our method.

###  A.1 Per-transit radial velocity measurement uncertainty

In the previous discussion, we assumed that the per-observation RV measurement
uncertainty $\sigma_{n}$ was known. This is not, however, a number that we
have access to. Instead, we derive a hierarchical model for determining a
calibrated data-driven estimate of $\sigma_{n}$, as a function of color and
apparent magnitude, based on the published data. Across a coarse grid in color
and apparent magnitude, we build a probabilistic model for the observed
distribution of RV error in each bin using the non-central $\chi^{2}$
likelihood defined by Equation A8 to allow for multiple star systems in the
training set. This model is high dimensional—there is a parameter for each
target in the bin!—so we use variational inference (e.g. Blei et al., 2017) as
implemented by the numpyro software package (Phan et al., 2019) to make the
model computationally tractable.

This procedure produces an estimate of the per-transit RV measurement
uncertainty $\sigma$ on a grid in color and apparent magnitude. We interpolate
this grid to produce an estimate of the measurement uncertainty for each RV
target in the Gaia catalog, and then use these results to compute $p$-values
and semi-amplitudes for this sample.

## References

  * Andrew et al. (2022) Andrew, S., Penoyre, Z., Belokurov, V., Evans, N. W., & Oh, S. 2022, Binary parameters from astrometric and spectroscopic errors, and candidate massive dark companions in Gaia eDR3, arXiv, doi: 10.48550/ARXIV.2206.04392
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. 2017, Journal of the American statistical Association, 112, 859
  * Carmichael et al. (2022) Carmichael, T. W., Irwin, J. M., Murgas, F., et al. 2022, arXiv:2202.08842 [astro-ph].  http://arxiv.org/abs/2202.08842
  * Cañas et al. (2022) Cañas, C. I., Mahadevan, S., Bender, C. F., et al. 2022, The Astronomical Journal, 163, 89, doi: 10.3847/1538-3881/ac415f
  * Chen & Kipping (2021) Chen, Z., & Kipping, D. 2021, arXiv e-prints, arXiv:2112.00966.  https://arxiv.org/abs/2112.00966
  * Chontos et al. (2022) Chontos, A., Murphy, J. M. A., MacDougall, M. G., et al. 2022, The Astronomical Journal, 163, 297
  * Ciardi et al. (2015) Ciardi, D. R., Beichman, C. A., Horch, E. P., & Howell, S. B. 2015, ApJ, 805, 16, doi: 10.1088/0004-637X/805/1/16
  * Clark et al. (2011) Clark, B. M., Blake, C. H., & Knapp, G. R. 2011, The Astrophysical Journal, 744, 119, doi: 10.1088/0004-637x/744/2/119
  * Coughlin et al. (2016) Coughlin, J. L., Mullally, F., Thompson, S. E., et al. 2016, ApJS, 224, 12, doi: 10.3847/0067-0049/224/1/12
  * Czekala et al. (2019) Czekala, I., Chiang, E., Andrews, S. M., et al. 2019, ApJ, 883, 22, doi: 10.3847/1538-4357/ab287b
  * Duquennoy & Mayor (1991) Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485
  * El-Badry et al. (2021) El-Badry, K., Rix, H.-W., & Heintz, T. M. 2021, Monthly Notices of the Royal Astronomical Society, 506, 2269, doi: 10.1093/mnras/stab323
  * Fischer & Valenti (2005) Fischer, D. A., & Valenti, J. 2005, The Astrophysical Journal, 622, 1102, doi: 10.1086/428383
  * Foreman-Mackey et al. (2016) Foreman-Mackey, D., Morton, T. D., Hogg, D. W., Agol, E., & Schölkopf, B. 2016, The Astronomical Journal, 152, 206, doi: 10.3847/0004-6256/152/6/206
  * Gajdoš et al. (2017) Gajdoš, P., Parimucha, Š., Hambálek, Ľ., & Vaňko, M. 2017, MNRAS, 469, 2907, doi: 10.1093/mnras/stx963
  * Guerrero (2021) Guerrero, N. M. 2021, The Astrophysical Journal Supplement Series, 29
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Holman & Wiegert (1999) Holman, M. J., & Wiegert, P. A. 1999, The Astronomical Journal, 117, 621, doi: 10.1086/300695
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Katz et al. (2019) Katz, D., Sartoretti, P., Cropper, M., et al. 2019, A&A, 622, A205, doi: 10.1051/0004-6361/201833273
  * Kirk et al. (2016) Kirk, B., Conroy, K., Prša, A., et al. 2016, AJ, 151, 68, doi: 10.3847/0004-6256/151/3/68
  * Knutson et al. (2014) Knutson, H. A., Fulton, B. J., Montet, B. T., et al. 2014, The Astrophysical Journal, 785, 126, doi: 10.1088/0004-637X/785/2/126
  * Kostov et al. (2020) Kostov, V. B., Orosz, J. A., Feinstein, A. D., et al. 2020, AJ, 159, 253, doi: 10.3847/1538-3881/ab8a48
  * Kraus et al. (2012) Kraus, A. L., Ireland, M. J., Hillenbrand, L. A., & Martinache, F. 2012, The Astrophysical Journal, 745, 19, doi: 10.1088/0004-637X/745/1/19
  * Kraus et al. (2016) Kraus, A. L., Ireland, M. J., Huber, D., Mann, A. W., & Dupuy, T. J. 2016, \aj, 152, 8, doi: 10.3847/0004-6256/152/1/8
  * Kunimoto et al. (2022) Kunimoto, M., Daylan, T., Guerrero, N., et al. 2022, ApJS, 259, 33, doi: 10.3847/1538-4365/ac5688
  * Lissauer et al. (2014) Lissauer, J. J., Marcy, G. W., Bryson, S. T., et al. 2014, ApJ, 784, 44, doi: 10.1088/0004-637X/784/1/44
  * Maoz et al. (2012) Maoz, D., Badenes, C., & Bickerton, S. J. 2012, The Astrophysical Journal, 751, 143, doi: 10.1088/0004-637x/751/2/143
  * Martin & Triaud (2014) Martin, D. V., & Triaud, A. H. M. J. 2014, Astronomy & Astrophysics, 570, A91, doi: 10.1051/0004-6361/201323112
  * Martin & Lubow (2017) Martin, R. G., & Lubow, S. H. 2017, The Astrophysical Journal, 835, L28, doi: 10.3847/2041-8213/835/2/L28
  * Marzari & Gallina (2016) Marzari, F., & Gallina, G. 2016, Astronomy & Astrophysics, 594, A89, doi: 10.1051/0004-6361/201628342
  * Maxted & Jeffries (2005) Maxted, P. F. L., & Jeffries, R. D. 2005, MNRAS, 362, L45, doi: 10.1111/j.1745-3933.2005.00073.x
  * Moe & Kratter (2021) Moe, M., & Kratter, K. M. 2021, Monthly Notices of the Royal Astronomical Society, 507, 3593, doi: 10.1093/mnras/stab2328
  * Ngo et al. (2015) Ngo, H., Knutson, H. A., Hinkley, S., et al. 2015, The Astrophysical Journal, 800, 138, doi: 10.1088/0004-637X/800/2/138
  * O’Donovan et al. (2006) O’Donovan, F. T., Charbonneau, D., Mandushev, G., et al. 2006, ApJ, 651, L61, doi: 10.1086/509123
  * Offner et al. (2022) Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2022, arXiv:2203.10066 [astro-ph].  http://arxiv.org/abs/2203.10066
  * Paardekooper et al. (2008) Paardekooper, S.-J., Thébault, P., & Mellema, G. 2008, Monthly Notices of the Royal Astronomical Society, 386, 973, doi: 10.1111/j.1365-2966.2008.13080.x
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Penoyre et al. (2022a) Penoyre, Z., Belokurov, V., & Evans, N. W. 2022a, Monthly Notices of the Royal Astronomical Society, 513, 2437, doi: 10.1093/mnras/stac959
  * Penoyre et al. (2022b) —. 2022b, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/mnras/stac1147
  * Phan et al. (2019) Phan, D., Pradhan, N., & Jankowiak, M. 2019, arXiv e-prints, arXiv:1912.11554.  https://arxiv.org/abs/1912.11554
  * Piskorz et al. (2015) Piskorz, D., Knutson, H. A., Ngo, H., et al. 2015, The Astrophysical Journal, 814, 148, doi: 10.1088/0004-637X/814/2/148
  * Pourbaix et al. (2009) Pourbaix, D., Tokovinin, A. A., Batten, A. H., et al. 2009, VizieR Online Data Catalog, B/sb9
  * Price-Whelan et al. (2017) Price-Whelan, A. M., Hogg, D. W., Foreman-Mackey, D., & Rix, H.-W. 2017, The Astrophysical Journal, 837, 20
  * Price-Whelan et al. (2017) Price-Whelan, A. M., Hogg, D. W., Foreman-Mackey, D., & Rix, H.-W. 2017, ApJ, 837, 20, doi: 10.3847/1538-4357/aa5e50
  * Price-Whelan et al. (2020) Price-Whelan, A. M., Hogg, D. W., Rix, H.-W., et al. 2020, The Astrophysical Journal, 895, 2, doi: 10.3847/1538-4357/ab8acc
  * Prša et al. (2022) Prša, A., Kochoska, A., Conroy, K. E., et al. 2022, ApJS, 258, 16, doi: 10.3847/1538-4365/ac324a
  * Silsbee & Rafikov (2015) Silsbee, K., & Rafikov, R. R. 2015, ApJ, 798, 71, doi: 10.1088/0004-637X/798/2/71
  * Stassun et al. (2019) Stassun, K. G., Oelkers, R. J., Paegert, M., et al. 2019, AJ, 158, 138, doi: 10.3847/1538-3881/ab3467
  * Su et al. (2021) Su, X.-N., Xie, J.-W., Zhou, J.-L., & Thebault, P. 2021, The Astronomical Journal, 162, 272, doi: 10.3847/1538-3881/ac2ba3
  * Sun et al. (2019) Sun, L., Ioannidis, P., Gu, S., et al. 2019, Astronomy & Astrophysics, 624, A15, doi: 10.1051/0004-6361/201834275
  * Teske et al. (2021) Teske, J., Wang, S. X., Wolfgang, A., et al. 2021, The Astrophysical Journal Supplement Series, 256, 33
  * Uehara et al. (2016) Uehara, S., Kawahara, H., Masuda, K., Yamada, S., & Aizawa, M. 2016, The Astrophysical Journal, 822, 2, doi: 10.3847/0004-637X/822/1/2
  * Wang et al. (2014) Wang, J., Fischer, D. A., Xie, J.-W., & Ciardi, D. R. 2014, The Astrophysical Journal, 791, 111, doi: 10.1088/0004-637X/791/2/111
  * Wang et al. (2015) Wang, J., Fischer, D. A., Barclay, T., et al. 2015, The Astrophysical Journal, 815, 127, doi: 10.1088/0004-637X/815/2/127
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Wood et al. (2021) Wood, M. L., Mann, A. W., & Kraus, A. L. 2021, The Astronomical Journal, 162, 128, doi: 10.3847/1538-3881/ac0ae9
  * Ziegler et al. (2021) Ziegler, C., Tokovinin, A., Latiolais, M., et al. 2021, The Astronomical Journal, 162, 192, doi: 10.3847/1538-3881/ac17f6
  * Ziegler et al. (2018) Ziegler, C., Law, N. M., Baranec, C., et al. 2018, AJ, 155, 161, doi: 10.3847/1538-3881/aab042