Wide Binaries in Tycho-{\it Gaia}: Search Method and the Distribution of Orbital Separations

TL;DR

This study develops a Bayesian method to identify wide stellar binaries in the Tycho-Gaia catalog, revealing their distribution and properties, and demonstrating Gaia's capability to detect orbital velocities in very wide systems.

Contribution

We introduce a Bayesian search method for wide binaries in Gaia data and analyze their orbital separation distribution without assuming Galactic structure.

Findings

Identified 6196 high-confidence wide binary candidates.

At least 0.6% of TGAS stars have wide binary companions.

Distribution follows a power-law with index -1.6 at large separations.

Abstract

We mine the Tycho-{\it Gaia} astrometric solution (TGAS) catalog for wide stellar binaries by matching positions, proper motions, and astrometric parallaxes. We separate genuine binaries from unassociated stellar pairs through a Bayesian formulation that includes correlated uncertainties in the proper motions and parallaxes. Rather than relying on assumptions about the structure of the Galaxy, we calculate Bayesian priors and likelihoods based on the nature of Keplerian orbits and the TGAS catalog itself. We calibrate our method using radial velocity measurements and obtain 6196 high-confidence candidate wide binaries with projected separations $s\lesssim1$ pc. The normalization of this distribution suggests that at least 0.6\% of TGAS stars have an associated, distant TGAS companion in a wide binary. We demonstrate that {\it Gaia}'s astrometry is precise enough that it can detect projected orbital velocities in wide binaries with orbital periods as large as 10$^6$ yr. For pairs with $s\ \lesssim\ 4\times10^4$~AU, characterization of random alignments indicate our contamination to be $\approx$5\%. For $s \lesssim 5\times10^3$~AU, our distribution is consistent with \"{O}pik's Law. At larger separations, the distribution is steeper and consistent with a power-law $P(s)\propto s^{-1.6}$; there is no evidence in our data of any bimodality in this distribution for $s \lesssim$ 1 pc. Using radial velocities, we demonstrate that at large separations, i.e., of order $s \sim$ 1 pc and beyond, any potential sample of genuine wide binaries in TGAS cannot be easily distinguished from ionized former wide binaries, moving groups, or contamination from randomly aligned stars.

Figures

Tables

Table 1. Table 1: The astrometric parameters of the 12 open clusters from which we remove candidate pairs. RA, Dec, and R 𝑅 R provide the position and angular radius around which we identify stars associated with these clusters. μ RA subscript 𝜇 RA \mu_{\rm RA} , μ Dec subscript 𝜇 Dec \mu_{\rm Dec} , and Δ ​ μ Δ 𝜇 \Delta\mu define the corresponding proper motion constraint, and ϖ italic-ϖ \varpi and Δ ​ ϖ Δ italic-ϖ \Delta\varpi define the parallax of each cluster. Note that these are “by-eye” approximations to Gaia astrometry rather than precisely measured characteristics.

Cluster	RA	Dec	$R$		${\mu }_{\mathrm{RA}}$		${\mu }_{\mathrm{Dec}}$		$\Delta \mu$		$\varpi$		$\Delta \varpi$
	[deg.]	[deg.]	[deg.]		[mas yr-1]		[mas yr-1]		[mas yr-1]		[mas]		[mas]
Pleiades	03:46:00	$+$24:06:00	8.	6	20.	10	$-$45.	39	10.	0	7.	4	2.0
Coma Ber	12:24:00	$+$26:00:00	7.	5	$-$11.	75	$-$8.	69	10.	0	11.	53	3.0
Hyades	04:27:00	$+$15:52:12	18.	54	110.	0	$-$30.	0	30.	0	21.	3	5.0
Praesepe	08:40:00	$+$19:42:00	4.	5	$-$35.	81	$-$12.	85	10.	0	5.	49	2.0
$\alpha$ Per	03:30:00	$+$49:00:00	8.	0	22.	73	$-$26.	51	5.	0	5.	8	2.0
IC 2391	08:40:00	$-$53:06:00	8.	0	$-$24.	69	22.	96	10.	0	6.	90	2.0
IC 2602	10:40:48	$-$64:24:00	2.	9	$-$17.	02	11.	15	5.	0	6.	73	2.0
Blanco I	00:04:24	$-$30:06:00	2.	5	20.	11	2.	43	5.	0	4.	83	2.0
NGC 2451	07:41:12	$-$38:30:00	2.	0	$-$21.	41	15.	61	5.	0	5.	45	2.0
NGC 6475	17:53:36	$-$34:48:00	3.	5	2.	06	$-$4.	98	10.	0	3.	70	2.0
NGC 7092	21:31:36	$+$48:24:00	3.	5	$-$8.	02	$-$20.	36	5.	0	3.	30	2.0
NGC 2516	07:57:36	$-$60:42:00	3.	5	$-$4.	17	11.	91	7.	0	2.	92	2.0

Equations

$$P(a)\propto a^{-1}$$

$$\displaystyle P(e)\propto e\ \ \$$

$${\vec{x}}_i$$

$${\vec{x}}_j$$

$$\theta \approx \sqrt{({\alpha }_A-{\alpha }_B{)}^2\cos {\delta }_A\cos {\delta }_B+({\delta }_A-{\delta }_B{)}^2}.$$

$$\Delta \mu \approx \sqrt{({\mu }_{\alpha ,A}^{\ast }-{\mu }_{\alpha ,B}^{\ast }{)}^2+({\mu }_{\delta ,A}-{\mu }_{\delta ,B}{)}^2},$$

$$P(C_2|{\vec{x}}_i,{\vec{x}}_j)=\frac{P({\vec{x}}_i|C_2,{\vec{x}}_j)P(C_2|{\vec{x}}_j)}{P({\vec{x}}_i|{\vec{x}}_j)}.$$

$$P({\vec{x}}_i|{\vec{x}}_j)=\sum _{k=1,2}P({\vec{x}}_i|C_k,{\vec{x}}_j)P(C_k|{\vec{x}}_j).$$

$$P(C_1|{\vec{x}}_j)$$

$$P(C_2|{\vec{x}}_j)$$

$$\frac{N(N-1)}{2}$$

$$N{f}_{\mathrm{bin}}$$

$$A_1$$

$$A_2$$

$$P(C_2|{\vec{x}}_i,{\vec{x}}_j)=\frac{P({\vec{x}}_i|C_2,{\vec{x}}_j)f_{\mathrm{bin}}}{P({\vec{x}}_i|C_2,{\vec{x}}_j)f_{\mathrm{bin}}+P({\vec{x}}_i|C_1,{\vec{x}}_j)A_1\rho ({\vec{x}}_j)}.$$

$$P(C_1|{\vec{x}}_i,{\vec{x}}_j)=1-P(C_2|{\vec{x}}_i,{\vec{x}}_j)$$

$$P({\vec{x}}_i|C_2,{\vec{x}}_j)$$

$$P(\theta ,\Delta \mu ,\varpi |C_2)=P(s,\Delta V|C_2)\frac{1}{{\varpi }^2},$$

$$P({\vec{x}}_i|C_2,{\vec{x}}_j)$$

$$\Delta V\le \sqrt{\frac{G{M}_{\mathrm{tot}}}{s}}.$$

$$P(s,\Delta V|C_2)=\{\begin{array}{ll}0& \Delta V>\Delta V_{\max }(s)\\ KDE(s,\Delta V)& \Delta V\le \Delta V_{\max }(s)\end{array}$$

$$\Delta V_{\max }(s)=2.11{\left(\frac{s}{10^3\mathrm{AU}}\right)}^{-1/2}\mathrm{km}{\mathrm{s}}^{-1}$$

$$P(\varpi)=\left\{\begin{array}[]{lc}\frac{1}{2L^{3}\varpi^{4}}\exp\left[-\frac{1}{L\varpi}\right]&{\rm for}\ 0<\varpi,\\ 0&{\rm otherwise}.\end{array}\right.$$

$$\int \mathrm{d}\vec{y}P({\vec{z}}_j,|\vec{y},{\vec{z}}_i)P(\vec{y}|{\vec{z}}_i)\approx \frac{1}{N}\sum _k^{N}P({\vec{z}}_j|{\vec{y}}_k,{\vec{z}}_i),$$

$$P({\vec{x}}_i|C_2,{\vec{x}}_j)\approx \frac{1}{N}\sum _{k}P(s,\Delta V|C_2)\frac{P({\varpi }_k)}{{\varpi }_k^2}P({\varpi }_2^{\prime }|{\varpi }_k,\Delta {\mu }_k),$$

$$P({\varpi }_1^{\prime }|\varpi ,\Delta \mu )P(\Delta {\mu }^{\prime }|\varpi ,\Delta \mu )$$

$${\mu }_{\alpha ,1,k}^{\ast },{\mu }_{\delta ,1,k},{\varpi }_{1,k}$$

$${\mu }_{\alpha ,2,k}^{\ast },{\mu }_{\delta ,2,k},{\varpi }_{2,k}$$

$$P({\varpi }_2^{\prime }|{\varpi }_k,\Delta {\mu }_k)=\mathcal{N}\left({\mu }_{\alpha ,2,k}^{\ast },{\mu }_{\delta ,2,k},{\varpi }_k;{\mu }_{\alpha ,2}^{\prime },{\mu }_{\delta ,2}^{\prime },{\varpi }_2^{\prime },{\mathbf{\Sigma }}_2\right).$$

$$P({\vec{x}}_i|C_1,{\vec{x}}_j)=\int \mathrm{d}{\varpi }_1\mathrm{d}{\varpi }_2\mathrm{d}\Delta \mu P({\varpi }_1,{\varpi }_2,\Delta \mu ,{\vec{x}}_i|C_1,{\vec{x}}_j).$$

$$P({\vec{x}}_i|C_1,{\vec{x}}_j)$$

Full text

††pubyear: 2017††pagerange: Wide Binaries in Tycho-Gaia: Search Method and the Distribution of Orbital Separations–A

1 Introduction

Roughly half of all main-sequence stars are found in binary systems with orbital separations extending from $\sim$10 ${R_{\odot}}$ to $\sim$pc (Duquennoy & Mayor, 1991; Fischer & Marcy, 1992; Raghavan et al., 2010; Duchêne & Kraus, 2013). Unresolved binaries at the smallest separations are typically found by identifying radial velocity (RV) variations indicative of orbital motion, while binaries at intermediate separations can be identified by observing the stars’ motions through their orbits over many years (Raghavan et al., 2010), or by identifying in decades of archival data dips in a star’s luminosity caused by eclipses due to an unseen stellar companion (Rodriguez et al., 2016). At the largest separations, stars are typically associated by watching them move through the solar neighborhood in the same direction and at the same rate; these common proper motion pairs have separations from 102 to $\geq$104 AU (Luyten, 1971).

The orbits of wide binaries may represent a significant fraction of the size of the clusters within which most stars are born (Lada & Lada, 2003). Such widely separated binaries are likely to form during the dissolution of a single-age stellar cluster (Kouwenhoven et al., 2010; Moeckel & Clarke, 2011) or from the evolution of higher-order systems (Reipurth & Mikkola, 2012; Elliott & Bayo, 2016). Since these pairs may be dynamically removed or slowly separated from the cluster (Theuns, 1992; Fujii et al., 2012), the population of wide binaries can probe the physics of cluster dissolution (Kroupa et al., 2001; Kouwenhoven et al., 2010). Additionally, triples or higher order systems formed via the fragmentation of circumbinary disks (Bonnell & Bate, 1994) may evolve dynamically into independent binaries (Bate et al., 2002; Moeckel & Bate, 2010) or stable, hierarchical triple systems (Reipurth & Mikkola, 2012). The observed distribution of orbital periods is therefore the result of the physics of formation and of subsequent dynamical reprocessing (Kroupa, 1995; Jiang & Tremaine, 2010; Moeckel & Clarke, 2011).

Beyond $\approx$103 AU, the components of wide binaries are not expected to interact significantly over their lifetimes. These pairs can be considered mini-open clusters (Binney & Tremaine, 1987; Soderblom, 2010) and are uniquely powerful for answering questions that cannot be addressed with other stellar systems. For example, wide binaries consisting of white dwarfs, subdwarfs, or M-dwarfs with FGK main sequence companions have been used to calibrate metallicity scales and age-rotation relations (Lépine et al., 2007; Rojas-Ayala et al., 2010; Garcés et al., 2011; Chanamé & Ramírez, 2012; Li et al., 2014). Other white dwarf-main sequence pairs and double white dwarfs have been used to determine the initial-to-final mass relation for white dwarfs (Finley & Koester, 1997; Zhao et al., 2012; Andrews et al., 2015). At the widest separations, binaries are affected by interactions with other stars, giant molecular clouds, and the Galactic tide, each of which acts to disrupt them stochastically (Weinberg et al., 1987; Jiang & Tremaine, 2010). These binaries have been used as probes of the gravitational potential of their host, including the contribution of dark objects and of dark matter in general (Bahcall et al., 1985; Yoo et al., 2004; Hernandez & Lee, 2008; Quinn et al., 2009; Monroy-Rodríguez & Allen, 2014; Peñarrubia et al., 2016).

We define a bona fide wide binary as a system of two stars whose properties are consistent with being gravitationally bound to each other. We do not claim that the components of wide binaries presented in this work are demonstrably orbiting their common center of mass, as that would at least require measurement of the accelerations of the components or of their orbit on the plane of the sky (e.g., Tokovinin & Kiyaeva, 2016). Instead, we regard a pair as a genuine binary when, to the limits of the data, their phase-space coordinates are consistent with the hypothesis that they are gravitationally bound to each other.

This definition becomes crucial when considering how wide one can go in a search for wide binaries. Within the Milky Way, a pair of stars cannot remain gravitationally bound for long beyond the separation at which they become sensitive to the Galactic tidal field. This tidal limit (or Jacobi radius) is $\approx$2 pc for \approx$${M_{\odot}} stars within the solar neighborhood (Yoo et al., 2004; Jiang & Tremaine, 2010). Pairs wider than this limit do exist: they may be wide binaries that have been ionized, members of moving groups, the tidal remains of larger structures, etc. However, as we demonstrate later in this work, they overlap in observational space with a sizable population of (false) pairs that occur simply by chance. Given the large variety of astronomical problems that can benefit from the assembly of reliable samples of wide binaries, we argue that pairs beyond the Galactic tidal limit should not be confused with genuine wide binaries, and we limit our search to pairs with projected separations $\lesssim$1 pc.

The exact form of the distribution of bona fide wide binary separations can help to constrain their formation scenarios. Although some groups find that this distribution is flat in log space out to the largest separations (following what is known as Öpik’s Law; e.g., Öpik, 1924; Poveda et al., 2007; Allen & Monroy-Rodríguez, 2014), the current consensus favors a distribution that decreases at larger separations (Chanamé & Gould, 2004; Lépine & Bongiorno, 2007; Tokovinin & Lépine, 2012; Halbwachs et al., 2017). As of yet, however, only a few attempts have been made to link directly the observed orbital separation distribution with models: Kouwenhoven et al. (2010) and Moeckel & Clarke (2011) have both compared $N$-body simulations of cluster dissolution to the observed orbital separation distribution. One reason for this scarcity of comparisons between data and simulations is that the observed samples are typically small. Identifying large numbers of wide binaries is difficult, and samples often suffer from significant and hard-to-characterize observational biases.

Pioneering work identifying wide binaries focused on small, specialized catalogs of stars (Luyten, 1971; Wasserman & Weinberg, 1991; Garnavich, 1991; Gould et al., 1995). More recent efforts have employed methods of varying complexity to systematically search larger, astrometric catalogs such as those produced by the New Luyten Two-Tenths survey (NLTT; Luyten, 1979) and the revised NLTT (rNLTT; Gould & Salim, 2003; Salim & Gould, 2003), e.g., Allen et al. (2000); Chanamé & Gould (2004), the Sloan Digital Sky Survey (SDSS; York et al., 2000), e.g., Sesar et al. (2008); Quinn & Smith (2009); Dhital et al. (2010); Longhitano & Binggeli (2010); Andrews et al. (2012); Baxter et al. (2014); Andrews et al. (2015); Dhital et al. (2015), the Lépine-Shara Proper Motion-North catalog (LSPM; Lépine & Shara, 2005), e.g., Lépine & Bongiorno (2007), and Hipparcos (Perryman et al., 1997), e.g., Shaya & Olling (2011); Tokovinin & Lépine (2012). These searches typically fall in one of two categories: there are those where stars are matched using only position (for instance using a two-point correlation function) and those where proper motion is included as well (to identify common proper motion pairs), corresponding to matching two or four dimensions, respectively, of the six-dimensional phase space.111Note that several groups have included photometric distance estimates in their matching algorithm. Although not used within the matching algorithm itself, there have also been efforts to include RVs (e.g., Latham et al., 1984), astrometric parallaxes (e.g., Lépine & Bongiorno, 2007), or both (Halbwachs et al., 2017) as extra phase-space dimensions to confirm the binary nature of wide pairs.

Shaya & Olling (2011) performed the first large-scale search for binaries that included Hipparcos parallaxes in a matching algorithm. The additional phase-space dimension adds a powerful constraint for separating genuine binaries from randomly aligned stellar pairs. However, because of its non-linear dependence on distance and non-trivial prior probabilities, the inclusion of parallax also adds a significant degree of complexity to the matching algorithm. Nevertheless, using a Bayesian formalism, these authors identify over 1000 binaries, including several high-probability pairs with projected separations of $\sim$pc.

The first data release from the Gaia satellite (Gaia Collaboration et al., 2016a, b; Lindegren et al., 2016) includes $2\times 10^{6}$ Tycho-2 stars (Høg et al., 2000), for which a Tycho-Gaia astrometric solution (TGAS; Michalik et al., 2015) is provided. The TGAS catalog is nearly complete to $V=11$ mag, with stars out to $V\approx 12$. Each of these stars has positions (with uncertainties of $\approx$0.4 mas), proper motions (uncertainties $\approx$1.3 mas yr*-1*), and parallaxes (uncertainties $\approx$0.3 mas, although a systematic bias of $\approx$0.25 mas has been identified in Gaia parallaxes; Stassun & Torres, 2016).

Oelkers et al. (2016) and Oh et al. (2016) performed searches for wide binaries and moving groups in the TGAS catalog, including for subsets with separations $>$1 pc. The initial Oelkers et al. (2016) catalog includes 8516 candidate pairs in which both stars are in TGAS. Although these authors use a sophisticated Galactic model to determine the random alignment likelihood for any particular wide binary, they do not calculate a corresponding likelihood that any binary could be a genuine pair. Furthermore, while their method combines proper motions, parallaxes, and their associated uncertainties to identify real pairs, it neglects correlations between the proper motion and parallax uncertainties, which can be considerable.

Oh et al. (2016) use an improved method to identify co-moving stars. These authors calculate the likelihood that a particular pair is a genuinely associated binary and that it is a random alignment of unassociated stars. In calculating these two likelihoods, these authors assume a model for the stellar components of the Milky Way.

We describe a Bayesian method for identifying wide binaries that combines positions, proper motions, astrometric parallaxes and their correlated uncertainties. Similar to Oh et al. (2016), we calculate two likelihoods for each stellar pair. However, our calculations use a data-based approach and are independent of any assumptions about the structure of the Milky Way.222Furthermore, Oh et al. (2016) select only those pairs with a parallax signal-to-noise (S/N) ratio $>$8, whereas we search the entire TGAS catalog. We separate genuine binaries from randomly aligned pairs using a posterior probability that combines these two likelihoods with prior probabilities based on the local stellar density. We further improve on previous work by verifying and calibrating our algorithm using two independent methods. First, using our algorithm, we identify and characterize a set of random alignments by matching the catalog to a shifted version of itself, which produces a sample containing only random alignments. Second, we use the subset of our pairs with radial velocity (RV) measurements to determine the rate of contamination due to randomly aligned stars. Our discussion is focused on identifying pairs within the TGAS catalog, but our method is designed to be adaptable to any stellar catalog with five-dimensional astrometry and is scalable to the expanded catalogs expected from future Gaia data releases.

In Section 2, we outline the challenges posed by the size of the TGAS catalog and the quality of Gaia data. Our statistical method for combining positions, proper motions, and parallaxes is described in Section 3. Before applying it to the full TGAS catalog, we describe in Section 4 our efforts to predict the characteristics of the contamination due to false binaries. In Section 5, we assemble our candidate TGAS pairs, calibrate our method to set a probability threshold for identifying real binaries, and estimate our contamination by false binaries. We present the resulting catalog of wide binaries, and use a low-contamination subset to characterize the projected separation distribution in Section 6. In Section 7, we compare our sample to pre-TGAS catalogs of wide binaries, discuss the impact of our assumptions on our search, address the question of whether we can reliably identify pairs at separations $\buildrel>\over{\sim}$1 pc, and provide some guidance to users of our catalog. We conclude in Section 8 with a few thoughts about identifying wide binaries in future Gaia catalogs. Appendix A is a test of our method on the subset of the rNLTT catalog with Hipparcos parallaxes; we compare the resulting catalog with a previously identified catalog of rNLTT wide binaries.

2 Challenges

Hipparcos provided the only data comparable to the TGAS data, for $\approx$105 stars. Since the number of possible pairs of a catalog of size $N$ is the arithmetic mean $N(N-1)/2$, identifying binaries in Hipparcos data implies the examination of \approx$$5\times 10^{9} individual pairs of stars. While the TGAS catalog is only 20 times larger than the Hipparcos catalog, the total number of pairs it contains is \approx$$2\times 10^{12}, a factor of 500 increase, and the final Gaia catalog will contain \approx$$5\times 10^{17} individual pairs. Of course, not every stellar pair need be checked (for instance, one can select a maximum angular search radius around each star), but the \approx$$N^{2} scaling with catalog size is likely unavoidable. Clearly the algorithms developed to search previous astrometric catalogs, even those designed to include Hipparcos parallaxes, need to be rethought.

Furthermore, Gaia’s superb astrometric precision creates new challenges for identifying wide binaries. Given this precision, the traditional approach of selecting pairs of stars with identical proper motions will miss a large number of true binaries: Gaia is sensitive to the differential orbital velocities of many wide binaries.

To illustrate, we generate 104 binaries with separations $a$ randomly drawn from a flat distribution in log space333Although evidence has indicated that the distribution of wide binary separations follows a power-law with a somewhat steeper index than $-$1 (Chanamé & Gould, 2004; Lépine & Bongiorno, 2007), adopting this flat distribution more clearly demonstrates the region in $s-\Delta V$ space populated by binaries. and eccentricities $e$ drawn from a thermal distribution:

[TABLE]

We further adopt random binary orientation angles to determine the distribution of the projected separations $s$ and of the projected differences in orbital velocities $\Delta V$.

Figure 1 compares the distribution of our binaries in $\log s-\log\Delta V$ space to Gaia’s astrometric limits. The TGAS data can be used to identify pairs of stars with angular separations as small as $\approx$2*′′* and proper motions as small as $\approx$1.3 mas yr*-1*. We apply these limits to pairs at three representative distances, 10, 100, and 500 pc, which allow us to translate the angular units observed by Gaia into physical units. In the top panel of Figure 1, we show the corresponding areas in $\log s-\log\Delta V$ space for which Gaia can measure $\Delta V$, the difference between the two velocity components of a binary due to orbital motion.

By the final Gaia data release, the astrometry will be substantially more precise: on-ground processing may allow identification of stellar pairs with separations as small as 0.1*′′* (Harrison, 2011), and proper motions may reach a precision of $\approx$0.25 mas yr*-1*. Using these limits, we generate the rectangular areas in the bottom panel of Figure 1. For most of wide binaries within 100 pc, and for a substantial fraction of those within 500 pc, $\Delta V$ must be taken into account.

Figure 1 also highlights one of the primary motivations for characterizing the population of wide binaries: a substantial fraction ($\approx$50%, assuming Öpik’s Law) have orbital separations $>$102 AU. Those binaries are too widely separated for a full orbit to be observed, and the only currently robust method to associate these pairs involves matching their position in phase space. Previous searches for such pairs were limited to nearby stars with large proper motions. With Gaia, we now have the ideal data to further explore this region of binary parameter space with significantly larger samples.

Finally, Gaia presents an additional, unprecedented opportunity to access all six dimensions of phase space. Its final RV catalog will achieve a precision of $\approx$10 km s*-1*, enough to aid substantially in separating genuine binaries from randomly aligned pairs. To our knowledge, RVs have been used exclusively in the past to check the binarity of wide binary candidates (e.g., Latham et al., 1984; Halbwachs et al., 2017); no RV-related algorithms have yet been included in common proper motion pair searches. For the first time, large numbers of binaries will be identified by their full three-dimensional positions and velocities.

3 Statistical Method

Figure 2 shows the positional (top panel) and proper motion (bottom panel) distributions of stars in the TGAS catalog. The density of stars is clearly a function of position and proper motion, and the probability of identifying two randomly aligned, unassociated stars as a binary depends strongly on the pair’s position in phase space, so not taking these biases into account may lead to serious problems. At the same time, the nature of Keplerian orbits combined with an assumption about the orbital separations of wide binaries provides us with an expectation about the properties of a distribution of real binaries. In this section, we quantify the likelihood that any stellar pair is both a random alignment and a true binary. The combination of these likelihoods with their associated prior probabilities provides the Bayesian posterior probability that any pair is a true binary.

As the details of our model are fairly technical, readers may wish to read the overview in Section 3.1 and then skip to the discussion of contaminants in Section 4. A test of our method on an existing catalog of wide binaries drawn from the rNLTT is presented in Appendix A.

3.1 Overview of the Method

Figure 3 illustrates part of our analysis for two separate pairs. The top row is for a pair we identify as a true binary, while the bottom row is for a pair that we determine to be a random alignment. The left column compares the astrometric parallaxes (based on their uncertainties) of the two stars in each pair as red and blue Gaussian curves, while the middle column shows the corresponding proper motion uncertainty ellipses for both stars in each pair. Concentric ellipses indicate correlated 1, 2, and 3$\sigma$ regions. The system in the top row clearly has matching astrometry, while in the bottom row the parallaxes and proper motions are only consistent at the $\approx$2$-$3$\sigma$ level.

To generate the plots in the last column of Figure 3, we transform each pair’s angular separation and each star’s parallax and proper motion to $s$ and $\Delta V$. The distribution in these two parameters for both pairs of stars is indicated by the contours, which correspond to 1, 2, and 3$\sigma$ confidence levels, while the blue background in these panels is the expected distribution of binaries simulated in Section 2. The convolution in $s-\Delta V$ space of the observed with the expected distribution provides the likelihood that the difference in a pair’s proper motions is consistent with binary motion. This likelihood, along with a statistic measuring the consistency of the two stars’ parallaxes, forms the probability that a pair’s astrometry matches that of a genuine binary. We construct this term mathematically in Section 3.3.

This is one part of the calculation in a Bayesian formalism; we also need to determine the likelihood that any given pair is consistent with being formed from a random alignment of stars. The likelihood function is based on the probability that any two stars can have the observed angular separation, proper motion difference, and parallax; random alignments are more likely to have larger angular separations, larger proper motion differences, and smaller parallaxes. This is one of the crucial aspects for the reliable identification of wide binaries, and we describe how we obtain this term in Section 3.4.

Finally, we need to calculate the prior probabilities of any given pair being a genuine binary and a random alignment of stars. Since the number of stellar binaries in a sample scales with the size of the sample, $N$, while the number of randomly aligned pairs scales with $N(N-1)/2$, any arbitrary pair has a very strong prior of being a random alignment. The exact value of this prior, which depends on the local density of a pair’s position and proper motion, needs to be determined individually for every pair. We describe how these priors are calculated in the following section.

3.2 Constructing the Bayesian Formalism

We consider any two stars occupying similar positions in five-dimensional phase space (position: $\alpha$, $\delta$; proper motion: $\mu_{\alpha}$, $\mu_{\delta}$; and parallax: $\varpi$) to be formed from one of two classes: random alignments ($C_{1}$) or genuine binaries ($C_{2}$). This probability will, in general, depend on all five parameters. For example, random alignments are more likely to be in dense stellar regions and to involve stars with small proper motions. To account for variations in this probability, we first assume that we are only interested in pairs with similar enough positions and proper motions that there is no difference in the average stellar density in the vicinity of each star. This is a reasonable assumption for all but the most widely separated stars. This assumption allows us to separate the two sets of astrometric parameters (one for each star) into $\vec{x}_{j}$, the pair’s position in four-dimensional phase space, and $\vec{x}_{i}$, which contains the scalar difference between the two stars’ positions and proper motions as well as each stars’ measured parallax:

[TABLE]

We use primes to indicate observed quantities that have some non-negligible uncertainty associated with them: uncertainties in the angular separation, $\theta$, can be ignored while uncertainties in $\mu$ (and hence $\Delta\mu$) and $\varpi$ cannot.

For computational efficiency, we take the components of $\vec{x}_{j}$ to be the primary star’s position and proper motion. The components of $\vec{x}_{i}$ will be determined by the difference between the two measured values; each component (save for $\theta$) will have an uncertainty associated with it. For closely separated stars $\theta$ can be determined from the two stellar coordinates:

[TABLE]

$\Delta\mu$ can be similarly calculated:

[TABLE]

where $\mu^{*}_{\alpha,i}=\mu_{\alpha,i}\cos\delta_{i}$.

Using these sets of parameters and the two classes $C_{1}$ and $C_{2}$, we can now use Bayes’s theorem to construct the generalized probability that any pair forms a true binary:

[TABLE]

The first term in the numerator of Equation 7 is the likelihood; the second term is the prior probability of a binary with $\vec{x}_{j}$. The denominator is the evidence, which can be determined by summing the probability that the particular pair can be produced by both a binary ($C_{2}$) and a random alignment ($C_{1}$):

[TABLE]

We determine the prior probabilities, $P(C_{1}\,|\,\vec{x}_{j})$ and $P(C_{2}\,|\,\vec{x}_{j})$, based on our understanding of binaries and random alignments: $P(C_{1}\,|\,\vec{x}_{j})$ should scale with the square of the stellar density in position and proper motion phase space, $\rho^{2}(\vec{x}_{j})$, while $P(C_{2}\,|\,\vec{x}_{j})$ should scale with $\rho(\vec{x}_{j})$ if we include both resolved and unresolved pairs:

[TABLE]

We determine the coefficients of these equations by integrating the relations over all four dimensions of $\vec{x}_{j}$, which will provide us with the total number of random alignments in the sample, $N(N-1)/2$, and the total number of binaries in the sample (including unresolved pairs), $f_{\rm bin}N$, where $f_{\rm bin}$ is the binary fraction:

[TABLE]

The integral to determine $A_{1}$ must be calculated numerically for a specific sample, while $A_{2}$ immediately becomes $f_{\rm bin}$:

[TABLE]

where we have used the substitution $\rho(\vec{x}_{j})=NP(\vec{x}_{j})$. We calculate the integral over $A_{1}$ using Monte Carlo random sampling over kernel density estimates (KDEs) of the distribution in position space and in proper motion space. We describe these density estimators in detail in Section 3.4.

We can now substitute the priors into Equation 7 to obtain our posterior distribution:

[TABLE]

We note here that $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ is normalized in such a way that the corresponding posterior probability for a particular binary being a random alignment can be expressed as:

[TABLE]

3.3 Binary Likelihood: $P(\vec{x}_{i}\,|\,C_{2},\vec{x}_{j})$

We now determine the probability that a true binary could produce the observations, $P(\vec{x}_{i}\,|\,C_{2},\vec{x}_{j})$. We begin by accounting for quantities with observational uncertainties. We marginalize over the true proper motion difference $\Delta\mu$, and the true binary parallax, $\varpi$. Both components of a true binary can be approximated as having the same distance. This is because the widest binaries have separations $\sim$1 pc, less than the uncertainty on distances derived from parallax (using Gaia’s nominal parallax uncertainty of 0.3 mas) for stars at distances beyond $\approx$60 pc. Since, as we will show in Figure 5, most binaries are found at significantly smaller separations and larger distances, we can simplify the problem by only marginalizing over the true parallax of the binary rather than each star’s parallax independently:

[TABLE]

We have separated out the observables $\varpi^{\prime}_{1}$, $\varpi^{\prime}_{2}$, and $\Delta\mu^{\prime}$ since they are dependent only on their true underlying values and their correlated uncertainties.

We can factor the first term in the integrand of Equation 18 so $P(\theta,\Delta\mu,\varpi\,|\,C_{2})=P(\theta,\Delta\mu\,|\,C_{2},\varpi)\ P(\varpi)$. Our understanding of binaries does not directly yield distributions in $\theta$ and $\Delta\mu$. However, binary theory and observations provide $P(s,\Delta V\,|\,C_{2})$, the expected distribution of the projected physical separation, $s$, and the tangential velocity difference, $\Delta V$. These are the physical forms of the angular variables $\theta$ and $\Delta\mu$: $\theta=s\varpi$ and $\Delta\mu=\Delta V\varpi$. Calculating the first term in Equation 18 then amounts to performing a Jacobian transformation from $(\theta,\Delta\mu)$ to $(s,\Delta V)$:

[TABLE]

where one factor of $\varpi$ enters in the denominator for each of $s$ and $\Delta V$ due to the Jacobian transformation. Care must be taken to ensure that these two factors of $\varpi$ are in the correct units: one factor should convert from $\Delta V$ to $\Delta\mu$ and the other should convert from $s$ to $\theta$. We can now reduce Equation 18 to:

[TABLE]

We now focus on the first term in the integrand in Equation 20, $P(s,\Delta V\,|\,C_{2})$, which expresses the likelihood that a random binary would produce the observed $s$ and $\Delta V$. In general, this function depends on assumptions made about populations of binary stars. We assume a binary is completely determined by four parameters: the two stellar masses, $M_{1}$ and $M_{2}$, the orbital separation, $a$, and the eccentricity. We do not include any binary evolution interactions between the two stars.

We determine the density of true binaries in $s-\Delta V$ space by randomly generating 104 binaries from two different distributions.444Our tests indicate that this is a sufficient sample size to balance accuracy with computation time. The first distribution uses a prior on $a$ that is flat in log space ($P(a)\propto a^{-1}$), following Öpik’s Law. We allow $a$ to range between 10 $R_{\odot}$ and 1 pc, so that half of all binaries have $a>$ 102 AU. The priors for the remaining orbital parameters required to produce distributions in $s$ and $\Delta V$ are described in Section 2; see also Figure 1.

We also test a model with a steeper power-law prior distribution: $P(a)\propto a^{-1.6}$. This is motivated by the results of Chanamé & Gould (2004, hereafter CG04) and Lépine & Bongiorno (2007), who find a steeper power-law-like decline in the distribution of binaries at large separations. As in the log-flat model, $a$ extends to 1 pc; the inner boundary of $\approx$33 AU is chosen so that half of all binaries have $a>$ 102 AU. Throughout the remainder of the paper, we use both orbital separation models to derive our results.

The differences between the models are substantial: the log-flat model places relatively more weight at larger separations, increasing the possibility of identifying binaries at wider separations. The cost is an increased rate of contamination at these large separations. We test samples of binaries obtained from both methods to obtain our conclusions on the orbital separation distribution of wide binaries in Sections 6.2 and 6.3 to ensure that our results are not dependent upon our adopted model assumptions.

From both populations of binaries, we use a KDE with a top hat kernel to create a normalized probability density function in $s-\Delta V$ space. Evaluating the KDE at a particular $s$ and $\Delta V$ provides $P(s,\Delta V\,|\,C_{2})$.

Circular binaries satisfy the inequalities that $s\leq a$ and $\Delta V^{2}\leq GM_{\rm tot}a^{-1}$, where $G$ is the gravitational constant and $M_{\rm tot}$ is the combined binary mass. These can be combined:

[TABLE]

Eccentric binaries do not obey this inequality since, for a fraction of the orbit near pericenter, the relative velocity of the two stars can substantially exceed the circular approximation. Furthermore, the mass of an arbitrary binary is unknown, so a fiducial total binary mass must be adopted, and the restriction on $\Delta V$ must be relaxed.

Nevertheless, we find that when non-circular binaries and a range of masses are taken into account, the observed $\Delta V$ for a given $s$ is typically lower than in the circular case. Guided by our population of randomly generated binaries, we set the maximum binary mass to 10 ${M_{\odot}}$ to account for flexibility in the eccentricity and total binary mass. We show the boundary of this inequality as a red line in Figure 3. Including this as an additional constraint greatly decreases computation time, as there are fewer calls to the KDE of $P(s,\Delta V\,|\,C_{2})$:

[TABLE]

where $\Delta V_{\rm max}(s)$ comes from Equation 21, where we set $M_{\rm tot}$ to 10 ${M_{\odot}}$:

[TABLE]

Returning to Equation 20, the $P(\varpi)$ term is the prior on the parallax measure. This probability distribution is not flat, which is the origin of the Lutz-Kelker bias (Lutz & Kelker, 1973). However, Smith (1987) demonstrated that the exact form of this distribution for a magnitude-limited sample differs from the original formulation in Lutz & Kelker (1973). And recently, Astraatmadja & Bailer-Jones (2016) showed that the distances to stars in the TGAS catalog are reasonably approximated using an exponentially decaying prior distribution with a scale length, $L$, of 1.35 kpc. This prior distribution, when represented in terms of $\varpi$, is:

[TABLE]

We use this form for $P(\varpi)$ throughout this work. An improved model could include the directional dependence of the parallax prior. For instance, the distribution of distances to stars is likely substantially different depending on whether we are looking out of the Galactic Plane or toward the Galactic Center. In practice, our results are largely unaffected by our choice of a parallax prior (see discussion in Section 7.2).

We now return to Equation 20. The integral on the right hand side of the equation is a multidimensional integral of the form: $\int\mathrm{d}\vec{y}\ P(\vec{z}_{j},\,|\,\vec{y},\vec{z}_{i})\ P(\vec{y}\,|\,\vec{z}_{i})$, where $\vec{z}_{i}$ and $\vec{z}_{j}$ are arbitrary variables. This integral can be approximated using Monte Carlo importance sampling:

[TABLE]

where $\vec{y}_{k}$ is drawn randomly $N$ times from the joint conditional probability: $\vec{y}_{k}\sim P(\vec{y}\,|\,\vec{z}_{i})$.

Adopting this approximation, Equation 20 reduces to:

[TABLE]

where $\varpi_{k}$ and $\Delta\mu_{k}$ are random samples drawn from their joint conditional probabilities. To obtain appropriately sampled $\varpi_{k}$ and $\Delta\mu_{k}$ based on $P(\varpi^{\prime}_{1}\,|\,\varpi,\Delta\mu)$ and $P(\Delta\mu^{\prime}\,|\,\varpi,\Delta\mu)$, we combine these two terms into one multivariate normal distribution ($\mathcal{N}$) and recognize its symmetric nature:

[TABLE]

where $\Delta\mu^{\prime}$ is dependent upon the proper motions for both stars.

We can now obtain random samples of $\varpi$ and $\Delta\mu$ by drawing random samples from both stars’ observed multivariate normal distributions and their corresponding covariance matrices ($\bm{\Sigma}$):

[TABLE]

Random samples of $\varpi$ are the set of $\varpi_{1,k}$ and random samples of $\Delta\mu$ are the scalar difference between each star’s proper motions calculated using the definition in Equation 6; we discard random samples of $\varpi_{2,k}$. Finally, we can calculate $P(\varpi^{\prime}_{2}\,|\,\varpi_{k},\Delta\mu_{k})$ in Equation 26 by recognizing that it can be inverted into a multivariate normal distribution, as was done in Equation 27. Now, we evaluate the multivariate normal distribution for star 2 at $\varpi^{\prime}_{2}$ conditional on the set of already generated proper motion samples, as well as $\varpi_{k}$ and $\Delta\mu_{k}$:

[TABLE]

The sum in Equation 26 provides the likelihood function. Convergence tests indicate that 105 Monte Carlo random samples of $\varpi$ and $\Delta\mu$ provides sufficient accuracy for our purposes.

3.4 Random Alignment Likelihood: $P(\vec{x}_{i}\,|\,C_{1},\vec{x}_{j})$

The likelihood that a pair of stars with $\vec{x}_{i}$ and $\vec{x}_{j}$ is the product of a random alignment is $P(\vec{x}_{i}\,|\,C_{1},\vec{x}_{j})$. We begin by marginalizing over the true individual parallaxes $\varpi_{1}$ and $\varpi_{2}$ and the true proper motion difference $\Delta\mu$ to account for observational uncertainties in these quantities:

[TABLE]

Now we can substitute for $\vec{x}_{i}$ and factor out $\varpi^{\prime}_{1}$, $\varpi^{\prime}_{2}$, and $\Delta\mu^{\prime}$, as observed quantities are dependent only on their underlying values (and their associated, correlated uncertainties):

[TABLE]

Our assumption that $\theta$ and $\Delta\mu$ are independent allows us to split the last term in the integrand. Of course, Galactic structure implies that position and proper motion are not entirely independent. However, accounting for this joint dependence adds a significant degree of complexity and computational expense to the problem. In the test of our method in Appendix A, we recover all the previously identified pairs in the rNLTT catalog, indicating any reduction of the effectiveness of our method due to this approximation is minor. Equation 32 becomes:

[TABLE]

This integral over $\varpi_{1}$, $\varpi_{2}$, and $\Delta\mu$ is of the form shown on the left hand side of Equation 25 and can therefore be calculated using Monte Carlo importance sampling:

[TABLE]

where $\varpi_{1,k}$, $\varpi_{2,k}$, and $\Delta\mu_{k}$ are $N$ samples generated from the probability distribution $P(\varpi^{\prime}_{1},\varpi^{\prime}_{2},\Delta\mu^{\prime}\,|\,\varpi_{1},\varpi_{2},\Delta\mu)$. As was done in Equation 27, this distribution is equivalent to a multivariate normal distribution, the arguments of which can be flipped:

[TABLE]

We generate these random samples in the same way as was done in Section 3.3. First, we generate random proper motions and parallaxes for each star using the observed multivariate normal distribution and its covariance matrix using Equations 28 and 29. $\Delta\mu_{k}$ is then determined from the scalar difference of randomly selected proper motions of each of the two stars using Equation 6.

To determine the $\theta$ term in Equation 34, we recognize that as $\theta$ increases, the probability of random alignments increases linearly. This is because the probability of a star having a randomly aligned companion at an angular separation $\theta$ is found from the integrated stellar density around an infinitesimally thin annulus of radius $\theta$. For a population with a linearly increasing density in space (or a uniform density), this probability is proportional to $2\pi\ \theta$ and the stellar density at the central star’s position. So long as changes in the stellar density does not substantially deviate from linearity on spatial scales larger than the typical binary separation, this approximation is appropriate.

Using an analogous argument, the $\Delta\mu$ term in Equation 33 should scale linearly with both $\Delta\mu$ and the density in $(\mu_{\rm\alpha^{*}},\mu_{\delta})$ space. However, Figure 2 shows that, contrary to what is seen in position space, the density in proper motion space changes substantially on relatively smaller scales, particularly at smaller $\mu$ where the distribution is strongly peaked. We tested our prescription by finding the typical $\Delta\mu$ at which the number of stars surrounding a point significantly deviates from our linear approximation; while most wide binaries have $\Delta\mu\sim 1-3$ mas yr*-1*, our approximation is typically valid for $\Delta\mu<5$ mas yr*-1*. Nevertheless, our locally linear approximation may become inaccurate in specific cases. We leave a higher-order extension of this probability for future work.

The resulting conditional probabilities for $\theta$ and $\Delta\mu$ can now be expressed as:

[TABLE]

where $\rho(\alpha,\delta)$ and $\rho_{\mu}(\mu^{*}_{\alpha},\mu_{\delta})$ are the positional and proper motion-dependent densities, respectively, and $N$ is the number of stars in the catalog and serves as a normalization constant.

To calculate these densities, we generate two separate KDEs with top hat kernels for stars in the TGAS catalog: one in $(\alpha\cos\delta,\delta)$ space, and one in $(\mu^{*}_{\alpha},\mu_{\delta})$ space. The KDEs act as smoothing functions that interpolate over the distributions in the two panels in Figure 2 and produce a normalized probability density function in position and $\mu$ space.

We use the KDE implementation in the Python package scikit-learn. This algorithm offers a substantial efficiency gain as it is based on a tree method, which has a computational time scaling with $O(N\log M)$, where $N$ is the number of function calls and $M$ is the number of elements in the tree, rather than the $O(NM)$ scaling of standard KDEs. We further optimize our algorithm by building the KDE tree based on a randomly chosen subset of 105 TGAS stars rather than for the full TGAS catalog. The size of this subset and the KDE bandwidth are optimized to minimize computation time while maintaining structure in the position and proper motion distributions.

The two parallax terms of form $P(\varpi)$ in Equation 34 represent the prior probability on parallax in the TGAS catalog. We use the exponentially decaying prior provided in Equation 24.

Using our Monte Carlo sampling from Equations 28 and 29, and substituting for the functions in Equations 24, 36 and 37, we can now use the sum in Equation 34 to calculate $P(\vec{x}_{i}\,|\,C_{1},\vec{x}_{j})$, the likelihood that any particular pair of stars is the product of a random alignment. Our convergence tests indicate that 105 random samples provide sufficient accuracy for our purposes here.

4 Characterizing Contamination

We address the problem of separating genuine wide binaries from randomly aligned stars in the TGAS catalog. Our approach is to first characterize the population of contaminating objects, which may be co-moving stars in kinematic structures such as open clusters, or, more likely, randomly aligned, unassociated stars.

4.1 Identifying and Removing Open Clusters

In addition to applying our statistical method constructed in Section 3, we remove any pairs consistent with being members of 12 open clusters: the Pleiades, Coma Ber, Hyades, Praesepe, $\alpha$ Per, IC 2391, IC 2602, Blanco I, NGC 2451, NGC 6475, NGC 7092, and NGC 2516. We rely largely on van Leeuwen (2009) to define the astrometric positions of these clusters, although we update the values based on Gaia’s improved astrometry (we adjusted the distance to the Pleiades, for instance). We remove stars that simultaneously satisfy a position, proper motion, and parallax constraint: they are within the cluster radius $R$ and have a proper motion and parallax within the cluster proper motion width $\Delta\mu$ and a parallax width of $\Delta\varpi$. Table 1 shows the parameters used to identify stars in these 12 open clusters. Note that these are not precisely measured cluster characteristics; they are “by-eye” approximations to the Gaia astrometry of several known open clusters.

Additional clusters and moving groups are likely to be found in the TGAS catalog; indeed, Oh et al. (2016) identify dozens of such co-moving groups composed of three or more elements in the TGAS catalog. Such structures typically span many square degrees on the sky, and are found with projected separations $>$1 pc. Our prior on the projected separation of wide binaries excludes co-moving stars with these separations from being identified by our algorithm, and we do not expect co-moving groups to be a strong source of contamination for our sample. We discuss this possibility further in Section 7.3.

4.2 Constructing a Sample of Random Alignments

To characterize the contamination from random alignments, we generate a sample of false pairs by taking the entire TGAS catalog and displacing it in both position and proper motion space, by $+$1∘ in $\delta$ and $+$3 mas yr*-1* in both $\mu^{*}_{\alpha}$ and $\mu_{\delta}$. We then match each star in this shifted catalog to the original TGAS catalog; the resulting matches are a catalog in which every pair is a random alignment. The size of this shift is chosen to be large enough that genuine pairs are not accidentally matched in our false catalog, but small enough that the local densities in position and proper motion space around each star are not changed substantially. This technique is an adaptation of that employed by Lépine & Bongiorno (2007).

Figure 4 shows the distances and projected separations for our sample of random alignments; the top row is the model with a log-flat prior on $a$, while the bottom row is for a power-law prior. We use the average parallax, weighted by the parallax uncertainties, to calculate the distance to each pair. The three columns show an increasing posterior probability from left to right; that random alignments populate the panel with a posterior probability of 99% indicates that contamination exists even at the highest of posterior probabilities. Importantly, random alignments may exist with separations down to the angular resolution of the instrument. The green and blue regions in Figure 4 define regions in distance and $\theta$ that are relatively free of contamination; we discuss these regions further in Section 5.1.

From Equation 16 the posterior probability of a particular pair being a random alignment is equal to $1-P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$. Our model assigns a probability $<$1% of being a random alignment to all those pairs in the rightmost panels of Figure 4, yet by construction none of them are binaries. Future improvements to our model may help address these imperfections. For instance, our model does not currently account for geometric effects from the non-Euclidean equatorial coordinates which can lead to slight problems in the proper motion matching for the most widely separated pairs, particularly at large declinations.555To illustrate this effect, imagine two stars with widely separated right ascensions at high declinations moving directly North. The two stars will have non-parallel proper motions. We assume that these two distributions are independent, but the proper motion distribution may vary with position. To account for this, a four-dimensional KDE could be generated based on both position and proper motion.

5 Constructing and Validating the Catalog

5.1 Identifying Wide Binaries

We now apply our method to the 2,057,050 stars in TGAS. Matching the entire catalog takes $\approx$300 CPU hours. Our algorithm is designed to be embarrassingly parallel, and future searches with larger catalogs can take advantage of large, high-performance computing clusters.

After removing several hundred binaries consistent with being in open clusters, our algorithm identifies 33,169 and 17,895 stellar pairs with posterior probabilities above 1% using the log-flat and power-law priors, respectively, on the orbital separation distribution, as described in Section 3.3.666Although our prior probabilities are calculated for pairs with $a\leq 1$ pc, non-zero eccentricities may allow us to detect some binaries with projected separations $>$1 pc.

Figure 5 shows the distribution of posterior probabilities for the two sets of pairs. The posterior probabilities tend to be either very close to unity or very close to zero. This indicates that although we must choose a critical posterior probability above which to define a candidate wide binary, the resulting catalog is relatively insensitive to that choice.

5.1.1 Distances, Angular Separations, and Physical Separations

Figure 6 shows the distance, angular separation, and projected physical separation of the candidate wide binaries identified by our method. Distances are determined for each pair from the average of the two stars’ parallaxes, weighted by their uncertainties. As in Figure 4, the different columns correspond to samples with different posterior probabilities.

Figure 4 and Figure 6 show remarkable similarities in their left columns, where binaries have posterior probabilities between 1% and 90%: in all four panels there is a locus of pairs whose density increases at both larger $s$ and larger distances. This suggests that the vast majority of pairs in this column in Figure 6 are due to random alignments.

Restricting the samples of false and real binaries to 90% $<P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})<$ 99% substantially reduces the total number of pairs, but their distributions in the middle panels of Figures 4 and 6 are still very similar. This similarity between the false and real binaries implies that vast majority of these are also random alignments of stars, although there may an excess of pairs at smaller $s$ in the real data, suggesting that some of these may be genuine binaries.

The locus of random alignments in the right column in Figure 4 at larger $s$ and larger distances also exists in our sample of candidate binaries in the right column of Figure 6. However, the distribution of our candidate binaries with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$% includes an additional locus of points with $s\lesssim 4\times 10^{4}$ AU. Since these data do not exist in our sample of random alignments in Figure 4, they are the genuine binaries in our sample. Furthermore, comparing Figures 4 and 6 shows that few genuine binaries are lost by requiring that pairs have $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$%.

Grey lines in Figure 6 show lines of constant angular separation. The TGAS catalog has a minimum angular separation cutoff at $\approx$2*′′, but there appears to be a substantial decrease in the population at separations $<$10′′. This is clearly seen by the sharp transition at 10′′* in Figure 7, which shows the angular separation distribution of our samples from the two methods. Since this decrease scales with $\theta$ rather than $s$, it cannot be an intrinsic property of wide binaries. Instead, this transition is most likely due to the input characteristics of the Tycho-2 catalog, rather than our search method. Future Gaia data releases will be sensitive to double stars within a fraction of 1*′′* (de Bruijne et al., 2015).

An additional test of our wide binary sample can be made by comparing the scalar proper motion difference of the two stars in each of the pairs. The left panels of Figure 8 show this difference, normalized by the proper motion uncertainties, estimated using the quadrature sum of the uncertainties on $\mu^{*}_{\alpha}$ and $\mu_{\delta}$ for each star, as a function of the projected separation. The set of binaries in our sample with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$% are shown both as points and the greyscale contours. The bulk of the pairs are consistent within 3 $\sigma$ (black horizontal line). However, it is immediately apparent that a substantial fraction of binaries (particularly those at smaller projected separations, where contamination by random alignments is negligible; see Figure 13) have component proper motions that deviate from each other at a greater than 3 $\sigma$ significance. This confirms our claim in Section 2 and predicted by Figure 1: Gaia’s astrometric precision is fine enough that it can detect the orbital motion of a wide binary with a separation of 104 AU, or equivalently an orbital period of 106 yr. Binaries with the most significant proper motion differences are typically Hipparcos stars with proper motion uncertainties smaller than $\approx 0.25$ mas yr*-1*. A wide binary sample produced from a method that does not account for this differential proper motion will fail to identify a large number of wide binaries.

In the right panels of Figure 8, we show the corresponding plot for the parallax difference between the two components of the wide binaries, normalized by the uncertainty in the parallax differences, as a function of the projected separations. As expected, most binaries have parallaxes that differ by less than 1 $\sigma$. There is evidence that the parallax difference increases somewhat at separations larger than 4$\times 10^{4}$ AU (vertical gray line), where the sample is expected to be dominated by random alignments. Stellar pairs with parallaxes inconsistent at the greater than 3 $\sigma$ level (but are nevertheless identified by our model with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$%) are indicated as red points in both the left and right panels.

5.1.2 RVs as a Test and Calibration of Our Method

While our method for identifying candidate wide binaries relies on five of the six dimensions of phase space, the components of genuine wide binaries should also stand out observationally in the remaining phase space dimension. They should share essentially the same RVs, whereas randomly aligned pairs should have a random distribution of RVs. Differences in the RVs of the components of a genuine wide binary should be of order $\Delta V$, which, as Figure 1 shows, is $<<$km s*-1* for all but the closest wide binaries in this sample.

The Fifth Data Release (DR5) of the RAdial Velocity Experiment (RAVE; Kunder et al. 2016) provides RVs for $>$250,000 Southern stars in TGAS. The typical accuracy of these RVs is better than 2 km s*-1*. We download the crossmatch between TGAS and RAVE included in RAVE DR5 and eliminate duplicate entries by retaining the observation with highest spectroscopic S/N. This produces a catalog of 210,368 stars, which we search for RV information for the random alignments constructed in Section 4.2.

In Figure 9, we compare the RVs of each star in the sample of random alignments generated using the power-law prior for the orbital separation distribution. The typical RV uncertainties are smaller than the plot symbols; clearly these pairs have inconsistent RVs. Figure 9 shows some non-zero covariance in the distribution of samples.777One possible reason for this may be due to our matching algorithm selecting members of the same kinematic component of the Galaxy, e.g. thin or thick disk. Also, we note that matching RVs is a necessary but not sufficient indicator of any particular pair being a genuine binary: Figure 9 shows that $\approx$20% of randomly aligned stellar pairs with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>$ 99% happen to have RVs consistent at the 3$\sigma$ level.

In Figure 10 we show, plotted against each other, the RAVE RVs for 583 pairs from our 17,895 candidate wide binaries obtained using the power-law prior. We divide the posterior probability $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ into five bins and give the number of pairs in each bin. The total number of pairs with RAVE RVs for both components is more than an order of magnitude larger than that identified by Oh et al. (2016) for systems with separations $<$1 pc.

The sequence of panels in Figure 10 shows that our candidate pairs progressively align along the line of equal RV as $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ increases, unambiguously telling us that genuine wide binaries are being identified with increasingly higher confidence, and thus validating our method. The top middle and top right panels, with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})<$ 50% and 90%, respectively, still look very similar to the plot showing random alignments. When restricting to $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>$ 90%, however, the behavior is very different, and a large majority of the 299 candidate pairs in this range have the RVs of the two components falling very close to the one-to-one line.

However, from our comparison between the distributions in Figures 9 and 10, candidate pairs with posterior probabilities between 90% and 99% may still include some random alignments, in agreement with our conclusions when comparing Figures 4 and 6. Therefore, to the limits of the available data, selecting candidates with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$% is likely to result in the purest sample of genuine wide binaries.

Figure 11 is the same as Figure 10, but for 863 pairs from our 33,169 candidate wide binaries obtained using the log-flat prior. Again, this sample is an order of magnitude larger than that found by Oh et al. (2016) for the same range of separations. In general terms, the same conclusions obtained from Figure 10 can be obtained from Figure 11, the main difference being the different relative numbers of candidate pairs falling in each probability range.

In the top left panel of Figure 12, we show the difference in RVs of the two stellar components from our sample of false binaries with RAVE RVs described in Section 5.1.2. These show that the population of random alignments increases at large $\theta$, and is typically of order tens of km s*-1*.888Given that use of the log-flat prior produces more contamination at low probabilities, we focus on the catalog obtained with the power-law prior for the distribution of orbital separations. Our conclusions are not altered if inspecting instead the catalog obtained using the log-flat prior.

Here again, the panels of Figure 12 corresponding to the catalog from the power-law prior show the progressive narrowing of the distribution of pairs around the equal velocity line (in this case, $\Delta$ RV = 0) as the posterior probability increases, with the overall behavior being markedly different on either side of $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})=90$% probability, i.e., a broad distribution of pairs for $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})<$ 90%, and a significantly narrower one for $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>$ 90%.

An additional difference is evident when looking at the distributions as a function of angular separation: candidate pairs are more broadly distributed at low probabilities than they are at the highest ones. Visual inspection of the panels in the upper and lower rows in Figure 12 reveals this without ambiguity. This is to be expected if the majority of pairs with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})<$ 90% are contamination from random alignments: the probability that two unassociated stars appear to be co-moving increases with the square of the angular separation. In this way, Figure 12 also argues that the catalogs are dominated by random alignments below 90% probability, and by genuine binaries at higher probabilities.

Finally, while the assumption that the components of a true wide binary share the same RV holds for most systems, there is one exception: the possibility that one of the components of a genuine wide binary is a (closer) binary system itself. Such hierarchical triples or multiples are not uncommon among wide binaries (Makarov et al., 2008; Tokovinin, 2014a, b, 2015; Halbwachs et al., 2017), and the orbit of the inner subsystem can make the RVs of the components differ from each other when the comparison is based on a single epoch. We further discuss this possibility, as well as other caveats to using RAVE RVs as a calibrator for our method, in Section 5.2. For now, we note only that these exceptions are a minority of cases, and matching the RVs of the components of a wide binary provides a satisfactory calibration method for our catalogs.

5.1.3 The Critical Posterior Probability

From all the evidence just discussed, it seems adequate to require $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>$ 90% for a candidate pair to be reliably regarded as a genuine wide binary. However, our goal is to identify as clean a sample of genuine wide binaries as possible, with no strong selection effects, and to characterize them as a population. Our comparison of sample distances, angular separations, and projected separations for the random alignments and our candidate binaries in Figures 4 and 6 demonstrate that most of the pairs with 90% $<P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})<$ 99% are random alignments. With this in mind, therefore, we require $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>$ 99% for the wide binary samples we report and analyze in this paper.

Even at $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>$ 99%, Figures 10, 11, and 12 show a small number of pairs that deviate significantly from the one-to-one and $\Delta$RV = 0 lines. As pointed out earlier, these pairs could represent contamination from chance alignments or be hierarchical triple systems. These should be followed up, ideally with multi-epoch RV measurements, but this is beyond the scope of the present work.

5.2 Estimating the Contamination Rate

The subset of pairs with RAVE RVs provides us with an estimate for the rate of contamination due to random alignments. In the left two panels of Figure 13, we show the distribution of projected separations for candidates with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$% as a function of their RV difference. This difference is normalized by the quadrature sum of the individual RV uncertainties, $\sigma_{\Delta{\rm RV}}^{2}=\sigma_{\rm RV,1}^{2}+\sigma_{\rm RV,2}^{2}$, so that deviations from zero correspond to confidence levels. Horizontal lines separate those pairs with consistent RVs at the 3$\sigma$ confidence level (blue points) from those with inconsistent RVs (orange points). The ratio of the number of candidate pairs with discrepant RVs to the overall number of candidates with RAVE matches provides an estimate of contamination. This method assumes that the subset of our candidate binaries with RAVE matches is a representative sample.

The right two panels of Figure 13 show the distribution of projected separations for pairs with consistent RVs (blue) and those with discrepant RVs (orange). Contamination exists at all separations, but the samples are clearly dominated by pairs with consistent RVs out to $s\approx 4\times 10^{4}$ AU. At larger separations, contamination begins to dominate both samples. Returning to the left two panels of Figure 13 and selecting only those systems with $s<4\times 10^{4}$ AU (delineated by the vertical gray lines), we estimate the contamination fraction at 17% for the log-flat prior (37/224) 16% for the power-law prior (33/212) models.

We can obtain an additional estimate of the contamination fraction using our sample of random alignments. The sample shown in Figure 4 produces approximately the same number of overall candidate pairs with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>1$% as our search for real wide binaries (both samples are dominated by random alignments with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})<<99$%). Restricting for those systems with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$% and $s<4\times 10^{4}$ AU, we estimate the contamination fraction using this method to be 6% (245/4360) and 4% (156/4024) for the log-flat and power-law prior models, respectively.

How can the two methods for estimating the contamination fraction disagree so strongly for the fraction of subset of systems with $s<4\times 10^{4}$ AU? Four factors can potentially make the RV-based estimate inaccurate. First, any star in a wide binary that includes an unresolved binary can make a genuine system appear as contamination using RVs. Makarov et al. (2008) suggest that the fraction of common proper motion pairs that are actually higher order systems may be as high as 25%. However, if only 10% of the systems in our sample were triples with an unresolved binary component, our contamination estimate would drop to $\approx$5%, in agreement with the estimate from our random alignments.

Second, a consistent RV does not guarantee binarity; even systems with RVs consistent at the 3$\sigma$ level may, in fact, be random alignments with similar, but slightly different RVs. Third, if the sample of candidate pairs for which we found matches in RAVE is not a representative TGAS sample, any estimate based on it may be inaccurate. Finally, if the RAVE RV uncertainties are somewhat underestimated, the contamination fraction would be overestimated.

Our contamination fraction estimate using randomly aligned pairs avoids these particular uncertainties. The procedure generating these false pairs is identical to that producing our catalog of candidate binaries. Furthermore, we find that restricting for those candidate pairs with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$%, both the false catalog and our catalog of wide binaries produces a comparable number of systems with $s>10^{5}$ AU, supporting the idea that the false pair catalog provides a reasonable estimate of contamination for pairs with $s<4\times 10^{4}$ AU. With these these caveats in mind, we estimate the contamination fraction to be $\approx$5%.

Regions 1 and 2 in Figure 6, outlined by the green and blue boxes, are designed to produce subsets of $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>99$% pairs within defined ranges of distances $D$ and angular separations that have very low contamination. Limited to $10{{}^{\prime\prime}}<\theta<100$′′ and $D<500$ pc, Region 1 takes advantage of the large sample of genuine binaries at smaller separations. Region 2 is confined to a distance of 100 pc, while extending to 10*′* to test the distribution of binaries at relatively wider separations. We estimate the contamination by counting the number of random alignments falling in these regions. For the log-flat prior model, our random alignment test generates 135 false pairs in Region 1 and four false pairs in Region 2, whereas our search identifies 2704 and 461 candidate pairs in these regions. This corresponds to contamination rates of $\approx$5% and 1%, respectively. For the power-law prior, the same test results in nearly identical contaminations rates of $\approx$4% and 1% for the two regions (89/2498 4/458). It is from these minimally contaminated subsets that we derive our constraints on the orbital separation distribution of wide binaries.

{sidewaystable}

The first 10 wide binaries from our catalog. We provide the Gaia source IDs, astrometry, posterior probabilities for both the log-flat prior and power-law prior models, and angular separation. This table is available online in its entirety.

Source ID Tycho-2 ID $\alpha_{1}$ $\delta_{1}$ $\mu^{*}_{\alpha,1}$ $\mu_{\delta,1}$ $\varpi_{1}$ Source ID Tycho-2 ID $\alpha_{2}$ $\delta_{2}$ $\mu^{*}_{\alpha,2}$ $\mu_{\delta,2}$ $\varpi_{2}$ $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ $\theta$

[mas yr*-1*] [mas yr*-1*] [mas]

[mas yr*-1*] [mas yr*-1*] [mas] log-flat power-law ["]

384307763771387776 2789-1319-1 00:00:10.11 +41:51:42.22 -7.34 -4.57 1.21 384313742365860992 2789-1122-1 00:00:46.53 +42:00:37.43 -6.46 -4.31 1.72 1.0 0.0 672.0

2765217352391138688 594-134-1 00:00:13.62 +09:47:56.43 26.82 -2.38 2.8 2765236731283414784 594-81-1 00:00:08.28 +09:56:05.50 27.47 -1.72 2.8 0.0 0.9973 495.4

2305900664854805632 7526-320-1 00:00:20.55 -41:02:27.77 48.7 -8.96 6.28 2305905578297389312 7526-515-1 00:01:39.22 -40:53:26.62 49.46 -7.99 5.31 0.9998 1.0 1042.4

429342935412328704 4014-3283-1 00:00:26.06 +60:25:39.51 28.74 11.84 10.21 429343004131805568 4014-1005-1 00:00:24.83 +60:25:31.28 28.11 11.49 10.28 1.0 1.0 12.3

2880169678167998336 2271-911-1 00:00:49.66 +36:46:48.02 -24.59 -22.27 6.82 2880169884326428672 2271-1988-1 00:00:48.62 +36:46:38.94 -25.78 -23.14 7.09 1.0 1.0 15.4

393882207946562176 3254-569-1 00:00:59.43 +49:50:54.46 6.05 4.77 2.75 393882311025775616 3254-1352-1 00:00:58.97 +49:52:02.10 7.02 4.13 2.09 0.9979 0.9864 68.7

4702974002115657472 9137-1774-1 00:01:17.61 -70:11:12.28 65.45 -12.93 4.67 4702974002115657600 9137-1708-1 00:01:22.50 -70:11:26.93 65.61 -13.54 4.61 1.0 1.0 28.2

531236430303749120 4298-598-1 00:01:18.66 +70:55:44.07 -5.31 -7.7 3.19 531236670821917568 4298-574-1 00:01:03.45 +70:55:47.55 -4.8 -6.92 3.44 0.9988 0.9927 74.6

2738660435728687872 1-1016-1 00:01:42.40 +01:47:12.51 -19.21 -26.34 2.66 2738660882405286272 1-168-1 00:01:41.56 +01:48:24.41 -19.67 -25.14 2.26 0.9998 0.8704 73.0

2848236321444116352 1729-1129-1 00:02:01.76 +23:46:51.13 26.99 5.65 5.04 2848423994335483264 1729-1118-1 00:01:59.86 +23:47:01.77 24.44 6.34 5.17 1.0 1.0 28.1

6 Results

6.1 The TGAS Catalog of Wide Binaries

We define wide binaries as systems with a posterior probability above 99% for either the log-flat or the power-law prior model, and our final catalog of 6196 pairs is the union of these two samples. The posterior probabilities may differ substantially between the models, particularly at large $s$ where the differences in the prior distributions is strongest.

We note that the posterior probabilities are derived from a calibrated and validated but ultimately imperfect model; pairs with $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})>$ 99% may still be contaminating objects. Indeed, at separations $s$ $\gtrsim 4\times 10^{4}$ AU, we suspect most pairs are random alignments regardless of their probability. We discuss this possibility further in Section 7.3. Nevertheless, we include these systems in the catalog: they define the locus of random alignments, and there are likely to be genuine pairs within this sample that can be extracted with a more focused analysis.

We present our catalog in Table 5.2. When individual stars match to multiple companions, each pair is listed as a separate binary; we do not use these systems when determining the wide binary orbital separation distribution (Sections 6.2 and 6.3). We provide the Gaia source IDs, Tycho-2 source IDs, positions, proper motions, and parallaxes of the components of each system, the posterior probabilities of each binary being composed of a genuine pair, and their angular separation.

6.2 The Orbital Separation Distribution of Wide Binaries

In Figure 14, we show the distribution of projected physical separations for subsets of pairs within two different distance limits for our log-flat and power-law prior models. These distributions are normalized to the total number of TGAS stars within the listed distance ranges. The small angle limit999Although our catalog of wide binaries contains pairs with separations as small as 2*′′, Figure 7 shows there is a substantial observational bias in effect at $\theta<10$′′*. for our sample leads to an excess of systems at $s\leq 2\times 10^{2}$ AU within $10^{2}$ pc (black) compared with our more distant sample (red). Importantly, the more distant binaries show an excess at $s\approx 10^{5}$ AU that is unseen in the closer sample at distances $\leq 10^{2}$ pc. This excess at large separations has been observed in other wide binary samples (Dhital et al., 2010; Oelkers et al., 2016) and suggested to be a characteristic of the Galactic wide binary population.

Were this increase a genuine characteristic of wide binaries, it should be seen in the distributions shown in Figure 14 of both the samples within 102 pc and that at distances between 102-$$2\times 10^{2} pc. Comparison between Figures 4 and 6 as well as the RV matching shown in Figure 13 demonstrates that this increase occurs at the separation where random alignments dominate the sample. We conclude that the relative flattening and subsequent apparent excess of the distribution at $s\approx 4\times 10^{4}$ AU in systems with $10^{2}$ < D < $2\times 10^{2}$ pc corresponds to the transition from genuine binaries to random alignments seen in the right panels of Figure 6. In particular, our data imply that the previously reported bimodality within separations of 1 pc in samples of wide binaries (Dhital et al., 2010; Oelkers et al., 2016) is most likely due to contamination from random alignments.

Additionally, the distribution of projected separations of our wide binaries, regardless of our assumed orbital separation prior, is in conflict with that found by Oh et al. (2016), who report an increasing distribution of co-moving pairs with $s$ ranging from a few 103 AU to several pc. However, our sample is in agreement with CG04 and Lépine & Bongiorno (2007), both of who find samples of wide binaries decreasing at separations larger than a few 103 AU. A major advantage of our sample over that of Oh et al. (2016) is that we explicitly identify and estimate the contamination fraction; after removing potential contamination our cleaned sample is nearly two orders of magnitude larger than that of Oh et al. (2016) for pairs with projected separations smaller than $4\times 10^{4}$ AU. If these authors introduced a strong selection function that distorts the observed distribution of stellar pairs, the difference between their sample and ours can be resolved.

Oh et al. (2016) further identify a set of co-moving pairs with $s\gtrsim 1$ pc, speculating that these may be the unbound remnants of ionized binaries whose components disperse in phase space during their orbital evolution. As Jiang & Tremaine (2010) demonstrate, such co-moving pairs may be prevalent in the nearby Galaxy. We discuss the nature of the population of pairs with separations $s\gtrsim 1$ pc in Section 7.3; for now we note that the distribution of wide binaries reported by Oh et al. (2016) increases in the same way as a population of purely random pairs, a sign of potential contamination.

Our sample of wide binaries shows a sharp cutoff at separations of a few 105 AU. This corresponds to the limits imposed by our binary model defined in Section 2: based on the expectation that (gravitationally bound) binaries do not survive for long at orbital separations beyond the Galactic tidal limit, our model has a strict upper limit that $a\leq$ 1 pc. As noted earlier, binaries with separations $>$1 pc may exist and indeed have been identified (e.g., CG04; Shaya & Olling, 2011). However, our results indicate that correctly identifying such widely separated systems as genuinely associated pairs is very challenging.

The dotted lines in Figure 14 correspond to the power-law distribution observed by CG04 and Lépine & Bongiorno (2007) with $P(s)\propto s^{-1.6}$. These lines are scaled vertically to compare the distribution to the power-law model at large separations. Random alignments artificially increase the distribution at large $s$ for the pairs with $10^{2}<D<2\times 10^{2}$ pc. Because the $D<10^{2}$ pc sample has less contamination at larger $s$, the power-law provides a reasonable model from \approx$$5\times 10^{3} to 105 AU. In the following section, we demonstrate that this decrease at large $s$ cannot be due to observational biases in detecting binaries with large $\theta$. The distribution of $s$ is clearly inconsistent with Öpik’s Law for s\ {\lower 3.44444pt\hbox{\buildrel>\over{\sim}}}\ 10^{4} AU, and it is not bimodal for $s\lesssim$ 1 pc.

6.3 Comparison to Models

To go beyond our by-eye comparison to the power-law separation distribution shown in Figure 14, we wish to use our full catalog of 6196 systems with posterior probabilities above 99%. The interdependence of distance, projected separation, and angular separation of these binaries shown in the right panels of Figure 6 demonstrates the difficulty in deriving the underlying distribution of physical separations. Contamination at large $s$, bias introduced by the distance distribution of stars in TGAS, and the angular resolution limit of TGAS at small $s$ may all cause inaccuracies when determining the distribution of $s$. Our approach is to generate a sample of synthetic wide binaries from an arbitrary binary population model, convert those binaries to an angular separation distribution, $P(\theta)$, accounting for the angular separation limits and distance dependence of TGAS. In this way, we are calculating $P(\theta)$ as predicted by the models to directly compare to subsets of our sample of binaries.

We restrict ourselves to Regions 1 and 2, defined in Figure 6, to calculate $P(\theta)$ expected from the model by integrating over the two oblique regions in $D$-$s$ space and then factoring terms based on independence:

[TABLE]

where $\mathcal{Z}$ is an integration constant to normalize the probability and $s_{\rm min}$ and $s_{\rm max}$ are the minimum and maximum projected separations, respectively, to which we are sensitive. For a given distance, these are limited by our angular separation constraints: at the close end by the angular sensitivity of the Tycho-2 catalog and at the distant end by our search radius of 100*′′* for Region 1 and 10*′* for Region 2.

Since $\theta$ can be expressed as a function of $D$ and $s$, Equation 38 can be reduced by substituting a Dirac delta function for $P(\theta\,|\,D,s)$:

[TABLE]

The delta function reduces the inner integral over $s$. After dividing by the Jacobian term, $1/D$, which comes from the derivative of the argument of the delta function, Equation 39 becomes:

[TABLE]

where $s^{*}=\theta D$. $P(D)$ is the distance distribution of our sample of binaries, and $P(s^{*})$ depends on the specific binary model used. We determine $P(D)$ empirically by generating a KDE from the subsample of our binaries with $5\times 10^{3}<s<5\times 10^{4}$ AU. This sample should be free of biases imposed by the TGAS angular separation limits and of contamination from random alignments at large $s$.

To generate $P(s)$, we create 105 binaries with randomly chosen orbital parameters and calculate the projected physical separation in the same way as described in Section 3.3. A thermal eccentricity distribution is assumed. Again, we generate $P(s)$ for two models, one with a log-flat distribution for the orbital separation in the range 10 {R_{\odot}}$$-1 pc, and one with a power-law of index $-$1.6. We approximate each of these distributions of 105 random binaries with gaussian KDEs to generate $P(s)$.

The integral in Equation 40 is normalized over the entire range of $\theta$. To obtain $P(\theta)$ over the range to which we are sensitive, we need to calculate the integration constant $\mathcal{Z}$:

[TABLE]

where $\theta_{\rm min}$ and $\theta_{\rm max}$ are the angular separation limits defining the region to which we are sensitive. To avoid incompleteness at small separations and random alignments at large separations, we set $\theta_{\rm min}$ to 10*′′* and $\theta_{\rm max}$ to 10*′* for Region 1 and $\theta_{\rm max}$ to 100*′′* for Region 2.

We compare the actual $P(\theta)$ for binaries in our two samples to these distributions. For both the log-flat and power-law distributions, we find $P(\theta)$ by numerically calculating the integral in Equation 40. Each panel in Figure 15 compares $P(\theta)$ for each of these two models. The top row of panels shows the distributions for our sample using the power-law model as our binary population prior, while the bottom panel shows that for the log-flat prior. The two columns compare the distributions for Region 1 extending to 10*′* and 100 pc and Region 2 extending to 100*′′* and 500 pc.

The left panels in Figure 15, calculated using those binaries in Region 1, focus on the distributions of pairs with small $\theta$. Although the two samples are very similar at these separations, the number of binaries in our sample allow us to discern between the two models. For \theta\ {\lower 3.44444pt\hbox{\buildrel<\over{\sim}}}\ 80′′, the population starts to have a shallower slope, closer to the log-flat distribution. This indicates that the turnover observed at separations of \approx$$5\times 10^{3} AU in Figure 14 is genuine.

The right two panels, calculated for binaries in Region 2, show that, regardless of the binary population model assumed, at separations out to 300*′′* the observed population matches the power-law distribution. Beyond 300*′′, the distribution deviates slightly, which may be due to the increased importance of random alignments (however, the uncertainties at separations beyond 300′′* are large). At smaller angular separations, the two models overlap.

Importantly, our results are identical regardless of the priors we use: both models show that $P(s)$ should be roughly log-flat at small separations and scale with roughly $s^{-1.6}$ at large separations. These results are consistent with the results of CG04 and Lépine & Bongiorno (2007), both of whom, using independent samples, find a flattening of the distribution at small $s$. For instance, Lépine & Bongiorno (2007) determines that $P(s)\propto s^{-1}$ for $s<2500$ AU while at larger separations, $P(s)\propto s^{-1.6}$; for $\theta<10{{}^{\prime\prime}}$, the log-flat and power-law distributions diverge further. Future Gaia data releases will be sensitive to these differences.

7 Discussion

7.1 Comparison with Pre-Gaia Wide Binary Samples

Searches of proper-motion catalogs for wide binaries have used different astrometric catalogs and were therefore sensitive to different samples of wide binaries. We briefly summarize the work of several of these groups before comparing their samples to that identified in this work.

CG04 searched the rNLTT catalog, a reanalysis of stars in the NLTT catalog (Salim & Gould, 2003), to identify wide binaries among common proper motion pairs. The stars in this catalog extend to $V\approx 19$ mag; the relative proper motion sensitivity is $\approx$3 mas yr*-1*.

Lépine & Bongiorno (2007) searched for common proper motion companions to Hipparcos stars in the LPSM-North catalog, a catalog identified using a reanalysis of the stars in the Digitized Sky Survey and containing stars with $\mu>0.15{{}^{\prime\prime}}$ yr*-1* (Lépine & Shara, 2005). Since the primary stars are found in Hipparcos data, this catalog contains the brightest and nearest stars.

Dhital et al. (2010) searched SDSS for common proper motion pairs containing low-mass stars. Proper motions are provided by the SDSS team by comparing to photometric plates from the USNO-B and POSS-II surveys (Munn et al., 2004). The digital nature of SDSS and the long time baseline results in proper motions with precisions $\approx$few mas yr*-1* for stars as faint as $g\approx 20$ mag. However, this sample is limited by several factors: source confusion prevents pairs with angular separations $\lesssim$7*′′* from being identified and SDSS has a bright limit of $g\approx 14$ and limited sky coverage.

Finally, Shaya & Olling (2011) searched the Hipparcos catalog for wide binaries. This sample includes only the brightest pairs, since the Hipparcos catalog is limited to $V\lesssim 9$ mag; however, the resulting catalog extends to proper motions as small as a few mas yr*-1*. Shaya & Olling (2011) were able to identify such slowly moving pairs because their algorithm also matched stars based on Hipparcos parallaxes.

The left four columns in Figure 16 compares the angular separations, magnitudes, and proper motions of pairs identified in the catalogs of CG04, Lépine & Bongiorno (2007), Dhital et al. (2010), and Shaya & Olling (2011).101010We only include those pairs identified by Shaya & Olling (2011) with a probability above 50%. The right column shows the distributions for the wide binaries identified in this work. The bottom panel in the rightmost column demonstrates that, like for Shaya & Olling (2011) and Lépine & Bongiorno (2007), our pairs tend to be at brighter magnitudes (Gaia $g\lesssim 11$ mag), reflecting the relatively bright faintness limit of the Tycho-2 catalog. However, even more so than Shaya & Olling (2011), we are identifying pairs with smaller proper motions.

The number of pairs in each catalog is provided in the top panel of each column. Our catalog represents an increase in size of a factor of $\approx$5-10 over each of the previous catalogs. At the same time, the CG04 and Dhital et al. (2010) catalogs shown in Figure 16 indicate that many more pairs exist with fainter magnitudes and large proper motions. Therefore, we can infer that the sample of wide binaries identified in this work represents a small subset of a much larger population that will be identified in future Gaia data releases.

7.2 The Impact of Our Assumptions

We discuss two of the three most important assumptions we made in constructing our catalog: the prior on the parallax and the binary fraction. The other assumption, about the input binary population model, was discussed in Section 3.3.

7.2.1 The Prior on the Parallax

Through their analysis of individual stars, Astraatmadja & Bailer-Jones (2016) test several prior distributions on parallax. These authors identify a model with an exponentially decreasing density with a scale length of 1.35 kpc as a reasonably accurate description of the parallax distribution. This is the prior we use, and we provide its functional form in Section 3.3.

As discussed by Lutz & Kelker (1973) and confirmed for Gaia data by Astraatmadja & Bailer-Jones (2016), the astrometric parallax is a poor measure of distance when its fractional uncertainty, $\sigma_{\varpi}/\varpi$, is $\buildrel>\over{\sim}$20%. Stars with small $\sigma_{\varpi}/\varpi$ (which are usually closer by) are less sensitive to the assumed parallax prior. Since the typical TGAS $\sigma_{\varpi}\leq 0.3$ mas, the parallaxes should be accurate for stars within $\buildrel<\over{\sim}$700 pc, which describes the vast majority of our wide binaries.

We test this assumption by running our algorithm with an empirically motivated parallax model based on the distribution of observed parallaxes. Specifically, we generate a KDE of a subset of parallaxes in the TGAS catalog to approximate this distribution. Uncertainties will tend to make a population of stars with small parallaxes appear somewhat closer, and therefore to have larger parallaxes, than at their true distances (Lutz-Kelker bias; Lutz & Kelker, 1973). Our empirical parallax prior is therefore an extreme case.

After comparing the results of our model to the two different parallax distributions, we find that the samples at separations $\buildrel<\over{\sim}$few 104 AU, from which we draw our conclusions about the properties of wide binaries, remain essentially unchanged. Since these very different parallax distributions produce nearly identical wide binary catalogs, we conclude that our results will be unaffected by even large differences in parallax prior probabilities. However, we note that this may not be true for future Gaia catalogs with different proper motion, parallax, and magnitude properties.

7.2.2 The Binary Fraction

Estimates of the binary fraction are a strong function of mass, and range widely, from 44% for solar-like stars (Raghavan et al., 2010) to $>$80% for O stars (Chini et al., 2012). However, these fractions are typically calculated for spectroscopic binaries. For wide binaries, Lépine & Bongiorno (2007) calculated that at least 10% of nearby Hipparcos stars have companions at separations $>$103 AU, and the wide binary fraction may be somewhat larger. Here, we assume that the total binary fraction $f_{\rm bin}$ is 50%, which enters into Equation 15 as the prior probability on $C_{2}$. The prior is normalized such that half of these stars, or 25% overall, have companions with separations $>$102 AU.

In addition, we define $f_{\rm bin}$ to be the fraction of binaries visible to Gaia, rather than the fraction of all stars in binaries. Due to the relatively bright faint limit of Tycho-2 photometry, many stars with faint companions cannot be identified in our data. If the wide binary fraction in the Galaxy is 50%, we would need a somewhat lower value for $f_{\rm bin}$ to account for our inability to detect pairs with faint companions. This bias will be substantially reduced in future Gaia data releases.

Fortunately, there is some freedom in our model to account for an inaccurate value for $f_{\rm bin}$, as it scales monotonically with our posterior probability $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ in Equation 15. If our adopted value for $f_{\rm bin}$ is too high, it will be reflected in the distribution of $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ for both genuine pairs and random alignments. In separating genuine binaries from random alignments using our independent RV and random alignment tests, we calibrate the critical value for $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ and account somewhat for our ignorance of $f_{\rm bin}$. Indeed, one of the reasons we find a critical value for $P(C_{2}\,|\,\vec{x}_{i},\vec{x}_{j})$ of 99% rather than a more intuitive 50% may be because our assumed value of $f_{\rm bin}$ is too high.

7.3 Stellar Pairs at pc-scale Separations

Stellar pairs with $s\gtrsim 1$ pc are quickly dissociated by the Galactic tide and are therefore not formally bound. However, these formerly bound pairs, which we consider co-moving pairs rather than binaries, may remain in similar locations in phase space for many orbits around the Galaxy (Jiang & Tremaine, 2010). Jiang & Tremaine (2010) find two peaks in the distribution of stellar separations, one well within 1 pc and one at several 100 pc. Importantly, regardless of their input model, these authors find that the trough in the distribution of stellar separations between these two peaks is at $\approx$5 Jacobi radii ($\approx$10 pc in the local Galactic neighborhood). There is a strong motivation to identify stellar pairs at separations $>$1 pc since detection of this trough followed by the second peak in the projected separation distribution of co-moving pairs would provide confirmation of the dissociating effect of the Galactic tide and potentially a powerful probe of Galactic structure.

Previous authors have claimed to identify large numbers of binaries or co-moving stellar pairs at wide separations (Shaya & Olling, 2011; Oh et al., 2016; Oelkers et al., 2016), and our search method identifies candidate binaries at projected separations $\sim$pc. Here we list several reasons that lead us to suspect that a large fraction of these stellar pairs, and by extension potentially those in other samples, are randomly aligned stars.

First, Figure 6 shows that these pairs have positions in log D$$-log $s$ space coincident with that of our sample of randomly aligned stars, derived from applying our method to a shifted version of the TGAS catalog. Second, the population of these pairs is strongly sensitive to the assumptions of the model. Our log-flat and power-law priors produce very different numbers of candidate pairs with $s\sim 10^{5}$ AU. This implies that the identification of any individual co-moving pair at these separations is tenuous and model-dependent. Third, Figure 14 shows that, regardless of the underlying model assumptions, the distribution of pairs at $s\sim 10^{5}$ AU is distance-dependent. This population is seen in the distribution of our sample at distances between 102 and $2\times 10^{2}$ pc, but not in the distribution within the nearest $10^{2}$ pc. However, there is no clear reason why genuine co-moving pairs with pc-scale separations might exist at distances beyond $10^{2}$ pc, but not within $10^{2}$ pc. Furthermore, in Figure 13 we show the difference in RVs of pairs with RAVE data as a function of $s$. Our results indicate that stellar pairs with $s\gtrsim 4\times 10^{4}$ AU identified from TGAS, even with high posterior probabilities, typically have inconsistent RVs.

Nevertheless, co-moving pairs at these wide separations may exist, and if they do, our method and others may classify them as candidates. For instance, Oh et al. (2016) identify a sample of co-moving pairs in TGAS with $s\gtrsim 1$ pc and demonstrate using matching with RAVE that the majority with a separation $<$few pc have consistent RVs. However, at separations $>$few pc, the RVs quickly become more consistent with random alignments. Furthermore, these authors find a distribution of projected separations linearly increasing out to $s\approx 10$ pc, the separation at which Jiang & Tremaine (2010) predict a trough in the distribution of genuine co-moving pairs. This increasing distribution with projected separation can be readily explained if the Oh et al. (2016) sample is dominated at these separations by randomly aligned stellar pairs: the number of random alignments increases with phase space volume, and therefore linearly with projected separation.

Oelkers et al. (2016) find several thousand co-moving pairs at separations $>$1 pc. As these authors estimate contamination based on an assumed Galactic model rather than a data-driven, independent method such as matching RVs, we suggest that at separations {\lower 3.44444pt\hbox{\buildrel>\over{\sim}}}4\times 10^{4} AU random alignments may be a significant source of contamination in their catalog, similar to what we see in our data.

Identifying wide stellar pairs that cluster in phase space, as is done by Oh et al. (2016), may be effective for detecting co-moving stellar conglomerations such as triple or higher multiplicity systems, open clusters, OB associations, and moving groups. However, the analysis we perform indicates that at separations beyond $\sim$1 pc it becomes increasingly difficult to separate isolated co-moving pairs from randomly aligned stars. We suggest that future searches for such pairs include careful accounting for contamination and, if possible, validation with RVs.

7.4 Using this Catalog

Before using the catalog provided in Table 5.2 there are a number of caveats to consider. First, the immense majority of pairs with posterior probabilities $>$90% are likely to be random alignments, and our results are derived from only those systems with posterior probabilities $>$99%. Nevertheless, there may be some contaminating random alignments in this sample. As discussed in Sections 5.1.1 and 5.1.2, randomly aligned pairs typically have $s>4\times 10^{4}$ AU, while pairs with $s<4\times 10^{4}$ AU are mostly genuine.

However, the transition between genuine pairs and random alignments is not sharp; genuine pairs with $s>4\times 10^{4}$ AU likely exist in our catalog, as may random alignments with very small separations. In particular, Figure 4 shows that random alignments may have angular separations as small as 1*′′*. To obtain a large and pure sample, we suggest using the subset of our candidate pairs in which the power-law posterior probability is above 99% and $s<4\times 10^{4}$ AU. Based on the tests described in Section 5.1.2, we expect the contamination of such a sample to be $\buildrel<\over{\sim}$5%.

Contamination due to random alignments becomes even less problematic for the candidate pairs falling within Regions 1 and 2 in Figure 6. Contamination is nearly non-existent for pairs in Region 2 in particular. However, these data are subject to important observational biases, particularly in distance. When deriving sample statistics where sensitivity to distance may be relevant, this bias should be taken into account. In Section 6.3, we describe one method for taking this into account using a Bayesian formalism.

Although we do not expect any bias in terms of $\theta$ to exist at separations larger than 10*′′, Figure 7 shows a marked decrease in the number of stellar pairs with $\theta<10$′′*.

8 Conclusions

8.1 Summary of Results

We present a method to identify wide binaries by matching pairs of stars based on their positions, proper motions, and parallaxes. We make no effort to identify co-moving stars in stellar clusters, moving groups, or higher multiplicity systems. We use a fully Bayesian method to discern between randomly aligned double stars and genuine binaries. Uncertainties on the astrometric quantities as well as their covariances are properly included in our analysis. We argue (and confirm using our resulting sample) that Gaia’s astrometric precision is fine enough to identify small differences in the proper motions of stars in wide binaries due to orbital velocity. We account for these slight velocity differences in our search algorithm. Furthermore, our method makes no assumptions about the structure of the Galaxy, and we calibrate it using independent RV observations.

As a test, in Appendix A we apply our algorithm to the subset of the revised NLTT catalog with parallaxes measured by Hipparcos (Gould & Salim, 2003). Chanamé & Gould (2004), whose search algorithm did not include parallax, identify 44 wide binaries in this sample. We re-detect these 44 binaries and find an additional 72 new wide binaries. This test demonstrates both the discerning power of including parallax as a matching criterion and the success of our method at parsing catalogs to identify wide binaries.

After applying our method to the 2$\times 10^{6}$ stars in the TGAS catalog, we identify 6196 wide binaries with a posterior probability above 99% in at least one of the catalogs produced using our two prior probability distributions for the orbital separation. Matching stars in our pairs with the RAVE catalog provides a subsample of our catalog that have RVs for both components. These RVs provide a test for our method: widely separated but gravitationally bound binaries should have components with matching RVs (modulo the presence of hierarchical triples). This test indicates that our catalog has a high fidelity and is largely free of unassociated, randomly aligned pairs.

When compared to catalogs of previously identified wide binaries, the pairs in this catalog are typically brighter, due to the magnitude limits of the Tycho-2 catalog. Nevertheless, our pairs extend to much smaller proper motions; where previous wide binary catalogs only include binaries with large proper motions, our sample contains binaries with proper motions as small as a few mas yr*-1*.

Our sample suffers from strong biases at both small separations, due to the angular separation limits of the Tycho-2 catalog, and at large angular separations, where contamination from randomly aligned pairs dominates the sample. Nevertheless, we define a method to compare an arbitrary distribution of binary separations with our biased sample.

Our sample indicates that at large physical separations, the distribution of binaries follows roughly a power-law with an index of $-$1.6, in agreement with the wide binary distributions identified by both Chanamé & Gould (2004) and Lépine & Bongiorno (2007). At smaller separations, the distributions flatten to a power-law index closer to $-$1 (Öpik’s Law). Importantly, these results are robust, as they are independent of whether our model assumes a log-flat or a power-law distribution as a prior on the orbital separation. A larger sample at a wider range of angular separations is needed to better determine the exact shape of this distribution.

Finally, we do not find a bimodality in the separations, as reported by previous authors (Dhital et al., 2010; Oelkers et al., 2016). Using RVs and comparing our catalog to a sample of random alignments, we demonstrate that the majority of candidate pairs identified by our method with $s\gtrsim 4\times 10^{4}$ AU are composed of randomly aligned, unassociated stars. We nevertheless include these pairs in our catalog, since a more focused analysis may be able to extract genuine binaries and co-moving pairs from this sample.

Future Gaia data releases will include RVs for stars brighter than $\approx$15 mag. Such measurements allow an important additional constraint on the binarity of widely separated pairs, potentially identifying wide binaries by matching their positions in full six-dimensional phase space. However, stars with matching parallaxes and proper motions but mismatched RVs may still be associated: the pair may actually be a hierarchical triple, with one star in the pair a close, unresolved binary. Engaging in RV follow-up of matched pairs with mismatched RVs may prove to be a unique and efficient method for identifying a population of hierarchical triple systems.

8.2 Application to Future Gaia Data Releases

Our method is designed to take advantage of Gaia’s precise astrometry while remaining scalable to larger data sets. Our current algorithm has a computational cost of a few 100 CPU hours for the 2$\times 10^{6}$ stars in the TGAS catalog. With modest optimization, this algorithm should be applicable to the 109 stars expected in the upcoming Gaia data releases.

Given the magnitude limits of the first Gaia data release, the pairs identified here are likely to have mass ratios close to unity. We therefore postpone any discussion of population statistics in terms of stellar masses for this sample. Upcoming data releases will extend the sample of stars with five-dimensional astrometry to $V\approx 20$. Using the yields from previous searches for wide binaries as an indicator, we estimate that searches through these future Gaia data sets will produce 104$-$105 wide binaries. This population will be large enough to identify robustly the underlying population of wide binaries as a function of mass and ultimately test various models for the formation of wide binaries.

Acknowledgements

We thank the anonymous referee for helpful suggestions which greatly improved the manuscript. J.C. thanks Claudia Aguilera-Gómez for introducing him to TOPCAT, which was extremely useful for working with the large TGAS and RAVE databases. We used the python module corner (Foreman-Mackey, 2016) to produce several figures in this manuscript. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013). J.J.A. acknowledges funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement n. 617001. J.C. acknowledges support from the Chilean Ministry for the Economy, Development, and Tourism’s Programa Iniciativa Científica Milenio, through grant IC120009 awarded to the Millenium Institute of Astrophysics (MAS), from PFB-06 Centro de Astronomia y Tecnologias Afines, and from Proyecto FONDECYT Regular 1130373. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://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. We acknowledge use of the Metropolis HPC Facility at the CCQCN Center of the University of Crete, supported by the European Union Seventh Framework Programme (FP7-REGPOT-2012-2013-1) under grant agreement no. 316165.

Appendix A Testing our Method with the Revised NLTT Catalog

Before applying our method to the TGAS catalog, we tested it against a set of wide binaries already identified in the literature from an astrometric catalog: the rNLTT catalog (Gould & Salim, 2003; Salim & Gould, 2003), which contains 36,085 entries with a typical $\mu$ uncertainty of several mas yr*-1*. The stars in rNLTT are mostly at latitudes above and below the Galactic Plane and were selected based on their large proper motions. Because of this low density in both position and proper motion space, the probability of randomly aligned pairs is substantially smaller in rNLTT than in TGAS.

Since our method requires parallaxes, we tested it only on the 8288 rNLTT stars with Hipparcos parallaxes. We found 213 pairs with a posterior probability $>$99% of being a genuine binary. However, a fraction of these pairs have both components of the binary matched to the same Hipparcos star; after removing these, we identified 116 wide binaries in the rNLTT catalog.

Chanamé & Gould (2004) found 999 wide binaries in the rNLTT catalog; 44 of these have Hipparcos parallaxes for both components. There are some important differences between our search algorithm and the approach of CG04. CG04 made proper motion cuts, using a limit for $\Delta\mu$ of 20 mas yr*-1*, a factor of three larger than typical uncertainties on $\Delta\mu$. Our method includes uncertainties in proper motion as well as a quantification of binary probability as a function of $\theta$, $\Delta\mu$, and position in position-proper motion phase space. On the other hand, to remove potential contaminants, CG04 calculated a modified reduced proper motion to discern between disk and halo populations. Clearly, true binaries cannot contain stars from different populations, a powerful constraint that is not built into our method (matching TGAS parallaxes does provide a similar, alternative constraint).

Given these differences, we did not expect the resulting binary samples to be identical. Encouragingly, however, we recovered all 44 CG04 pairs with Hipparcos parallaxes; the other 72 pairs in our sample are newly identified wide binaries. We compare the two samples in Figure 17.

We leave discussion and follow-up of these pairs for future work. For now we note that even in catalogs such as rNLTT that have been mined for wide binaries, there may be a substantial number of true pairs still to be identified.