The Separation Distribution of Ultrawide Binaries across Galactic Populations

It has long been recognized that wide binaries provide a powerful tracer of the Galactic gravitational potential on small scales (e.g., Bahcall & Soneira 1981; Bahcall et al. 1985), as well as useful constraints for studies of the star formation process (van Albada 1968; Chanamé 2007; Kouwenhoven et al. 2010; Moeckel & Bate 2010; Moeckel & Clarke 2011). The orbits of wide binaries are so fragile (with escape velocities v_esc ≲ km s⁻¹) that they can be easily unbound by external gravitational perturbations arising from molecular clouds, stars, or compact objects (e.g., Retterer & King 1982; Bahcall et al. 1985; Weinberg et al. 1987; Jiang & Tremaine 2010; Allen & Monroy-Rodríguez 2014), or by internal perturbations due to evolution of the component stars (Boersma 1961; Savedoff 1966; El-Badry & Rix 2018).

One of the most important applications of wide binaries is to constrain the mass of viable MAssive Compact Halo Objects (MACHOs), such as discrete black holes, which are a proposed dark matter candidate. Studies using wide binaries to constrain MACHOs have typically assumed a simple initial separation distribution (such as a smooth power law) and then interpreted deviations from it as arising from disruption of the widest binaries by MACHOs.

A number of previous studies have used halo wide binaries—which are expected to be less vulnerable to disruption by molecular clouds—to constrain the plausible parameters of MACHOs. Chanamé & Gould (2004, hereafter CG04) built a catalog of 1147 candidate wide binaries selected via common proper motion, from which they selected 801 binaries with disk-like orbits and 116 binaries with halo-like orbits. Fitting the two subsamples, they showed that the angular separation distributions could be well described by single power laws with logarithmic slopes of −1.67 ± 0.07 for the disk binaries and −1.55 ± 0.10 for the halo binaries, over the angular separation range 35 < Δθ < 900''. In a companion paper, Yoo et al. (2004, hereafter Y04) argued that the CG04 separation distribution constraints ruled out most of the previously plausible MACHO parameter space, since there was no strong evidence of a break in the separation distribution at the widest separations, and their simulations predicted a detectable break for MACHO masses M ≳ 43 M_⊙.

When combined with constraints from microlensing surveys (which ruled out a dominant MACHO population with typical masses in the range of 10⁻⁷ ≲ M/M_⊙ ≲ 30; Alcock et al. 2001; Afonso et al. 2003) and theoretical lower limits on the long-term survivability of hydrogenous objects (which rule out MACHOs with M ≲ 10⁻⁷ M_⊙; de Rujula et al. 1992), the Y04 wide binary constraints ruled out almost all possible MACHO masses. Only a small window of MACHO masses, 30 ≲ M/M_⊙ ≲ 43, was still considered viable. However, Quinn et al. (2009) later argued that the Y04 constraints depend critically on the validity of the two widest binaries in the CG04 sample. After removing one spurious candidate binary from the CG04 sample, they repeated the analysis of Y04 and found the wide binary separation distribution consistent with a MACHO mass of 30 ≲ M/M_⊙ ≲ 500.

Precise proper motions and parallaxes from Gaia DR2 (Gaia Collaboration et al. 2016, 2018) have allowed the construction of larger and purer wide binary catalogs than were available previously. Distance measurements from Gaia parallaxes have also made it possible to measure the physical separation distribution; this is more easily interpretable than the angular separation distribution, whose physical scale depends on distance. El-Badry & Rix (2018, hereafter ER18) searched Gaia DR2 for high-confidence wide binaries within 200 pc of the Sun. Their search yielded ∼55,000 high-confidence wide binaries, including ∼3500 in which at least one component is a white dwarf (WD). A similar search was carried out by Pittordis & Sutherland (2019), who used the projected orbital velocities of Gaia wide binaries to constrain the gravitational force law in the low acceleration regime. Both of these studies limited their search to projected separations of s < 50,000 au (0.25 pc), primarily because the contamination rate from chance alignments would have become nonnegligible at wider separations given their selection criteria. Because they searched for nearby binaries without any kinematic selection, the majority of the binaries in these samples have disk-like kinematics and roughly solar metallicity (see El-Badry & Rix 2019).

Fitting the separation distribution with a broken power-law model, ER18 found that for MS/MS binaries it was nearly consistent with a single power law of logarithmic slope γ ≈ −1.6 but displayed evidence of a weak break at $\mathrm{log}(s/\mathrm{au})\approx 3.8$. They found the separation distributions for binaries containing a WD to fall off more steeply at large separations and interpreted this steepening as evidence of nonadiabatic and/or asymmetric mass loss in the end stages of stellar evolution. The effects of disruption due to external perturbations are expected to become significant only at very wide separations (s ≳ 20,000 au; Weinberg et al. 1987), so ER18 did not attempt to constrain the effects of external perturbations using their sample.

In this paper, we (a) extend the ER18 analysis to wider separations and (b) compare the separation distributions of binaries between the Milky Way disk and halo. Constructing a large sample of halo binaries requires us to search to larger distances than ER18, but we show that it is still possible to obtain a high-purity binary sample at larger distances for halo stars, which have large proper motions and thus have fewer nearby neighbors in phase space. We kinematically select three pure subsamples with different average ages, using v_⊥,tot, the total tangential velocity with respect to the Sun, as a proxy of stellar age. By investigating the three binary populations, we try to answer three basic questions: (1) Is there a slope change in the separation distribution of wide binaries at s ≈ 10,000 au, as found tentatively by ER18? (2) If so, does its strength vary between stellar populations? (3) Can it provide a meaningful constraint or clues on formation or destruction mechanisms of binaries in different populations?

Most previous works using wide binaries as dynamical tracers have implicitly assumed that the primordial separation distribution can be well approximated as a simple power law at wide separations, and that it is the same for different populations. This is a serviceable assumption, and most previous works have found a power-law separation distribution to provide a good fit for wide binaries in the range of separations where external perturbations are expected to be subdominate (Lépine & Bongiorno 2007; Andrews et al. 2017), but there is little prior motivation for it. Now that the present-day separation distribution can be constrained in more detail, we reexamine this assumption.

The remainder of this paper is organized as follows. In Section 2, we describe how to build an initial wide binary candidate catalog that is extensive but not pure. Section 3 then presents three pure subsamples with different stellar populations selected from the candidate catalog, i.e., the disk-like, the kinematically intermediate, and the halo-like binaries. In Section 4, we illustrate the observed separation distributions for the three subsamples. In order to infer the intrinsic separation distributions, we describe the method in Section 5, including the selection function of the observed sample, a smoothly broken power-law parameterization, and the likelihood for fitting the separation distribution. Section 6 lists the key results of the investigation. We discuss possible theoretical interpretations of our results in Section 7 and summarize in Section 8.

Throughout the paper, we adopt the solar motion as (U_⊙, V_⊙, W_⊙) = (9.58, 10.52, 7.01) km s⁻¹ (Tian et al. 2015) and the circular speed of the local standard of rest (LSR) as v₀ = 238 km s⁻¹ (Schönrich 2012). α* is used to denote the R.A. in the gnomonic projection coordinate system, for example, ${\mu }_{{\alpha }^{* },i}\equiv {\mu }_{\alpha ,i}\cos {\delta }_{i}$ and ${\rm{\Delta }}{\alpha }^{* }\equiv {\rm{\Delta }}\alpha \cos (\delta )$.

We first construct a catalog of wide binary candidates that is large but not pure. In Section 3, we describe how pure subsamples can then be selected from the initial candidate catalog.

We select wide binary candidates with a procedure similar to that of ER18, whose catalog was restricted to pairs of stars within 200 pc of the Sun and projected separations s < 0.25 pc. They selected stars whose positions, proper motions, and parallaxes were consistent with being gravitationally bound. For the present study, we extend the search volume from 200 pc to 4.0 kpc and the maximum projected separation from s = 0.25 to 1.0 pc. We briefly summarize the selection procedure below.

2.1. General Criteria for the Wide Binary Query

2.2. Quality Cuts

2.3. Removing Clusters, Moving Groups, and Pairs in High-density Regions

2.4. Properties of the Initial Candidates

Figure 1 (top panel) shows a color–magnitude diagram (CMD) for all 807,611 pairs of wide binary candidates. The (primary) stars with M_G < 2.75(G_BP − G_RP) + 5.75 are designated as "main-sequence" (MS) stars (a few giants are not excluded), and the objects with M_G > 3.25(G_BP − G_RP) + 9.63 are classified as likely WDs. Here ${M}_{{\rm{G}}}=G\,+5\mathrm{log}\left(\varpi /\mathrm{mas}\right)-10$. The MS and red clump are visibly smeared out, primarily because the photometry is not corrected for extinction. Deviations from the expected CMD morphology are more pronounced for secondaries (top right panel of Figure 1) because the cuts on astrometric and photometric precision (Section 2.2) are looser for secondaries. A large cloud of objects is visible between the WDs and MS for secondaries. These are primarily sources with large parallax errors; they are removed from the sample once more stringent quality cuts are applied (see lower panels of Figure 1). A secondary sequence is apparent above the MS for both primaries and secondaries, indicating that the catalog contains hierarchical triples with components that are spatially unresolved close binaries.

Zoom In Zoom Out Reset image size

Figure 2 (top panel) presents the sky distribution of all binary candidates in the b versus l plane, where b and l represent galactic latitude and longitude. Clear imprints of the Gaia scanning law are apparent, as there are more sources with well-constrained astrometry in regions of the sky that were visited more often. A higher density of candidates is also apparent near the galactic plane; this was even more pronounced prior to the removal of sources with many phase-space neighbors. A few blank patches at b ∼ 0° are very dense regions where our queries timed out and no candidates could be identified.

Zoom In Zoom Out Reset image size

Figure 3 illustrates the distributions of distance ($\mathrm{log}(d)$, first column), magnitude difference (${\rm{\Delta }}G=| {G}_{1}-{G}_{2}|$, second column), apparent magnitude (G, third column), and angular separation ($\mathrm{log}\theta$, fourth column) for all wide binary candidates (black). The distribution of physical separation is displayed separately in Figure 4. For the initial candidate sample, the distributions are dominated by contaminants (chance alignments) at large separations (θ ∼ 42 and s > 0.05 pc). This can be seen in the distributions of $\mathrm{log}\theta$ and $\mathrm{log}(s)$, where the large-separation peak is due entirely to chance alignments. We discuss chance alignment further below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.5. Contamination Rates of the Initial Sample

We use three methods to quantify the contamination rates of the candidate binaries as a function of separation: (1) A "shifted Gaia DR2 catalog" where the coordinates ($\hat{\alpha }$, $\hat{\delta }$) for each object in this catalog have been constructed by shifting 1° from its original location (α, δ), i.e., ($\hat{\alpha }$, $\hat{\delta })=(\alpha +{\rm{\Delta }}{\alpha }^{* }/\cos (\delta )$, δ + Δδ), with Δα* = Δδ = 10. We then repeat the binary candidate identification procedure, now identifying pairs that pass the binary cuts when the coordinates of the "primary" are shifted relative to candidate secondaries. This procedure removes genuine binaries but preserves chance alignment statistics; see Lépine & Bongiorno (2007) for further discussion of this method of estimating the chance alignment rate. (2) The mock Gaia DR2 catalog produced by Rybizki et al. (2018), which is based on the Besançon stellar population synthesis model (Robin et al. 2003) and is populated with the Galactic distribution function using Galaxia (Sharma et al. 2011), assuming a similar selection function and uncertainty model to Gaia DR2. (3) Comparison of radial velocities (RVs) for the subset of bright wide binaries where both stars have RVs included in Gaia DR2.

The red dotted curves in Figure 3 present the distributions of the random alignment pairs chosen from the shifted Gaia DR2 catalog in the whole sky. We matched the original Gaia DR2 with the shifted catalog and regard the matched pairs as chance alignments if they pass the selection criteria for binaries. This approach is not appropriate for the widest and most nearby binaries, which can have projected separations of 1° or larger and so continue to appear as binaries even when the coordinates of one component are artificially shifted by 1° in each direction. Therefore, we remove the chance alignments with d < 60 pc and θ > 05 from the shifted catalog. The red dotted curves suggest that the distributions for the initial candidates are indeed almost fully dominated by contaminants at large angular and projected separations, as displayed in the last two subpanels of Figure 3. On the other hand, chance alignments are subdominant at log(s) ≲ 4.5 (Figure 3). The distribution of chance alignments can be well described as N ∼ s ds, since the area within which chance alignments can be found scales as 2πs ds.

We selected in a random region of ∼4000 deg² from the mock Gaia DR2 catalog (Rybizki et al. 2018). The mock catalog does not contain any true binaries, so any pairs that pass our selection criteria in the mock catalog must be chance alignments. We only query a fraction of the sky in the mock catalog because querying it is computationally expensive and many of the queries time out. In estimating the true contamination rate, we scale up the rate predicted by the mock catalog at all separations, assuming that the contamination rate approaches ∼100% at 1 pc separations.

In defining candidate wide binaries in this study, we used the two-dimensional velocities (i.e., 2D tangential velocities from proper motions), positions, and distances. The incomplete kinematic information induces substantial contamination, particularly at large separations. Gaia DR2 provides about 7.2 million stars with measured RVs. Among them, 8220 wide binary candidates have reliable (uncertainty σ_RV < 3 km s⁻¹) RVs for both binary components. These RVs provide the third dimension of velocity and can be effectively used to check whether candidate pairs are true binaries or chance alignments. Significant RV differences (ΔRV) between member stars imply that star pairs are probably chance alignments. We regard those pairs with ΔRV/σ_RV > 5 and ΔRV > 10 km s⁻¹ as likely contaminants. Figure 5 displays the one-to-one scatter distribution (top panel) of the two component RVs and the ΔRV distribution (bottom panel) as a function of projected separation. Only 20 pairs meet the criteria of chance alignments, i.e., ΔRV/σ_RV > 5 and ΔRV > 10 km s⁻¹. Most binaries with substantially large ΔRV have larger-than-average RV uncertainties. This implies that the contamination rate of the candidate binary catalog is of order 0.24%, with most contaminants at $\mathrm{log}(s/\mathrm{au})\gt 5$, for binaries that are bright enough for both stars to have measured RVs (G ≲ 13). However, we stress that the true contamination rate for the whole catalog of binary candidates is much higher, as can be seen in Figure 3. The RVs provide an estimate of the chance alignment rate for binaries where both components are bright.

Zoom In Zoom Out Reset image size

Table 1 summarizes the contamination rates in the different projected separation bins, estimated with the above three methods.

Table 1.  Contamination Rates of the Wide Binary Candidates and Cleaned Subsamples Estimated by Different Methods

Separations	Contamination Rates
log(s/au)	Mocka				Gaia RVsb				Shifted Gaia DR2c
	All	Disk	Intermediate	Halo	All	Disk	Intermediate	Halo	All	Disk	Intermediate	Halo
(0.00, 3.40]	0%	0%	0%	0%	0%	0%	0%	0%	<0.3%	0%	0%	0%
(3.40, 4.12]	<2%	0%	0%	0%	<0.1%	0%	0%	0%	<2%	0%	0%	0%
(4.12, 4.39]	∼10%	0%	0%	0%	<0.1%	0%	0%	0%	∼10%	0%	0%	0%
(4.39, 4.65]	∼30%	0%	0%	0%	<0.2%	0%	0%	0%	∼30%	0%	0%	0%
(4.65, 4.91]	∼80%	0%	0%	0%	∼0.4%	0%	0%	0%	∼80%	∼1%	∼0.1%	∼1%
(4.91, 5.20]	>90%	0%	0%	0%	∼1.5%	0%	0%	0%	>90%	0%	∼1%	∼5%
(5.20, 5.31]	>95%	0%	0%	0%	∼11%	0%	0%	0%	>95%	0%	∼10%	∼15%

Notes.

^aUsing the mock catalog from Rybizki et al. (2018) to estimate the contamination rate, assuming that the contamination for the initial candidates reaches 100% at the widest separations. ^bUsing the reliable (σ_RV < 3 km s⁻¹) Gaia RVs to estimate the contamination rate of the brightest wide binaries, where both components have G ≲ 13. Pairs with ΔRV/σ_RV > 5 and ΔRV > 10 km s⁻¹ are regarded as chance alignments. ^cUsing the shifted Gaia DR2 to estimate the contamination rate. Each object is shifted by 1° in each coordinate, i.e., ($\hat{\alpha }$, $\hat{\delta })=(\alpha +1\buildrel{\circ}\over{.} 0/\cos (\delta )$, δ + 10), where ($\hat{\alpha }$, $\hat{\delta }$) and (α, δ) are the shifted and original coordinates, respectively.

Download table as:  ASCIITypeset image

As low contamination is crucial for our subsequent analysis of the separation distributions, we proceed to define subsamples that we can assume to be nearly pure. Specifically, we kinematically select three pure subsamples, consisting of young disk-like, intermediate, and old halo-like binaries, with different average ages. Since the phase-space density of contaminants varies between the samples (with more contaminants in the low-velocity samples), we use different quality cuts for the three samples.

Using the Gaia DR2 proper motion and parallax, we calculate the total tangential velocity with respect to the Sun for each candidate binary:

$$\begin{eqnarray}&&{v}_{\perp ,\mathrm{tot}}\equiv 4.74\,\mathrm{km}\,{{\rm{s}}}^{-1}\times ({\mu }_{\mathrm{tot}}\times \mathrm{yr})/\varpi .\end{eqnarray} \tag{ 1 }$$

Here ϖ and ${\mu }_{\mathrm{tot}}=\sqrt{{\mu }_{{\alpha }^{* }}^{2}+{\mu }_{\delta }^{2}}$ are the parallax and total proper motion of a binary, respectively.

Using v_⊥,tot as a proxy of age, we select three pure subsamples with different average ages. We determine what cuts are needed to obtain a pure subsample in each population using the shifted catalog to estimate the chance alignment rate given any subselection.

The space density of halo stars in the solar neighborhood is much lower than that of disk stars, so it is necessary to search to larger distances to obtain a large sample of halo binaries. Fortunately, the contamination rate at fixed distance is also much lower for halo stars, because they have fewer phase-space neighbors. We thus use different distance limits for the three populations and self-consistently estimate the contamination rate for each population given these limits. We choose disk-like binaries with

• 1.  
  d < 120 pc,
• 2.  
  v_⊥,tot < 40 km s⁻¹.

To select the intermediate-age binary subsample, we use the following conditions:

• 1.  
  d < 300 pc,
• 2.  
  40 < v_⊥,tot < 85 km s⁻¹.

We choose old halo-like binaries with the following criteria:

• 1.  
  d < 700 pc,
• 2.  
  v_⊥,tot > 85 km s⁻¹.

To get pure binary samples, we further impose the following cuts on all of the above three selections:

• 1.  
  Δμ ≤ Δμ_orbit + 1.0σ_Δμ and σ_Δμ ≤ 0.12 mas yr⁻¹. This more stringent cut on proper-motion uncertainty reduces the contamination from chance alignments with larger uncertainty.
• 2.  
  The number of nearby neighbors N < 2, to strictly remove contaminants at wide separation from moving groups or star clusters.
• 3.  
  We exclude binaries containing a WD from all three subsamples to remove the effect from internal orbital evolution.

The above criteria can effectively suppress the three kinds of contamination rates close to zero for the three subsamples at each (except the largest) separation bin (see Table 1). Finally, we get 4361 disk-like, 10,090 intermediate, and 4351 halo-like genuine wide binaries. Figure 4 compares the separation distributions of binaries and chance alignments in the three subsamples. With these more aggressive cuts on distance and astrometric S/N, chance alignments are subdominant out to separations of 1 pc. The separation distributions of all three samples fall off monotonically at large separations.

Figure 1 displays the CMDs for the disk-like (blue), intermediate (green), and halo-like (magenta) wide binaries. The unresolved close binary sequences (ER18; Widmark et al. 2018) are visible but already sparse above the MS in the three subsamples. Following ER18, we use the region at 1 ≲ (G_BP − G_RP) ≲ 2 of the CMD (Figure 1) to estimate the unresolved binary fraction of the disk and halo samples. In this region, the line M_G = 2.8 (G_BP − G_RP) + 2.4 divides the binary and single-star MSs. We find that about 4.0% (4.3%), 10.9% (7.4%), and 6.6% (4.8%) of the primary (secondary) disk-like, intermediate, and halo-like binaries probably have a bright unresolved companion.

Figure 2 presents the final disk-like (the 4361 blue dots), intermediate (the 10,090 green dots), and halo-like (the 4351 magenta dots) binary distribution in the b versus l plane. Figure 3 illustrates the distributions of distance ($\mathrm{log}(d)$; first column), magnitude difference (${\rm{\Delta }}G=| {G}_{1}-{G}_{2}|$; second column), magnitude (G; third column), and angular separation ($\mathrm{log}\theta$; fourth column), respectively, for the three subsamples. For demonstrating the contamination rate, the distributions of projected separations ($\mathrm{log}(s)$) for the three subsets are specifically displayed in Figure 4. In this figure, the blue, green, and magenta dotted curves show the distributions of the chance alignments selected from the shifted catalog with the disk-like, intermediate, and halo-like criteria, respectively.

We use the mock Gaia DR2 catalog (Rybizki et al. 2018) to illustrate the effect of v_⊥,tot for distinguishing the three subsamples. The mock catalog contains the full 6D phase information for each star, so the disk-like, intermediate, and halo-like stars are easily separated. We randomly select 150,000 stars in the whole sky from this mock catalog and use the criteria of v_⊥,tot to separate the mock stars into disk-like, intermediate, and halo-like populations. Figure 6 shows the Toomre diagram for disk-like stars (blue points), intermediate (green points), and halo-like stars (magenta points). As one can see, the disk-like, intermediate, and halo-like stars are clearly divided. This indicates that the criteria of v_⊥,tot work well for the selection of the three kinematic subsamples.

Zoom In Zoom Out Reset image size

We display distributions of the angular (θ) and the projected physical separations (s = (θ × (1/ϖ)) au) for the disk-like, intermediate, and halo-like binaries in Figure 3 (the last subpanels) and Figure 4. At θ < 15, there are no binaries in our sample, which is due to the Gaia angular resolution limit. At 1.5 < θ/arcsec < 6.0 or 25 < s/au < 1000 (250 < s/au < 15,000 for the halo-like binaries), the samples are significantly incomplete. ER18 investigated this incompleteness effect and concluded that the incompleteness is due to blending of the two stars, which is more severe for binaries with magnitude difference between the two components. According to Figure A1 of ER18, the separation distribution is 84% complete for ΔG < 5 at θ > 5'' and is complete for all ΔG at θ > 10'', relative to the completeness at arbitrarily large angular separations.

Consistent with ER18, we find that the MS binaries do not obey Opik's law (i.e., a flat distribution of log separation) at any s ≳ 500 au. It is worthwhile to mention that there is an obvious drop at s ∼ 10^5.0 au in the raw separation distributions of the three subsamples, as shown in Figure 4.

We note that when insufficiently stringent astrometric quality cuts are applied, the separation distribution of all binary candidates in Figure 4 appears bimodal (black histogram). However, the peak at large separations is due entirely to chance alignments, not real binaries. For this reason, it shifts toward wider separations and eventually disappears as more stringent quality cuts are applied. Some previous works (Dhital et al. 2010; Oelkers et al. 2017) have found the separation distribution of wide binaries to be bimodal and have proposed that the wider-separation population contains binaries formed by different formation mechanisms. Our analysis suggests that this apparent bimodality is simply due to chance alignments, as also suggested by Andrews et al. (2017).

As discussed in Section 4, the separation distribution at θ < 10'' is heavily incomplete. Before modeling the intrinsic separation distribution, one must account for the selection effect to compensate for the incompleteness at small separation. In this section, we directly adopt an empirical fitting function from ER18 as the selection function, and we use a smoothly broken power law with four free parameters to model the intrinsic separation distribution.

5.1. Selection Function

The selection function depends on the angular separation θ and the magnitude contrast ΔG between the two member stars. For this study, we empirically characterize the probability of detecting a companion at an angular separation θ with a fitting selection function f_ΔG(θ) found by ER18,

$$\begin{eqnarray}&&{f}_{{\rm{\Delta }}G}\left(\theta \right)=\displaystyle \frac{1}{1+{\left(\theta /{\theta }_{0}\right)}^{-\beta }}.\end{eqnarray} \tag{ 2 }$$

Here ΔG is the magnitude contrast in G band between the two member stars, θ₀ characterizes the angular separation below which the sensitivity drops to 0, and β determines how rapidly the sensitivity falls off at θ ≪ θ₀. Following ER18, we adopt β = 10 for all ΔG, θ₀ = 225 at ΔG < 1.5 mag, and θ₀ = 0.9(ΔG + 1) at ΔG > 1.5 mag. We stress that f_ΔG(θ) does not represent the absolute selection function, but rather the fraction of binaries detected at angular separation θ relative to the number that would be detected at arbitrarily large separation. See El-Badry et al. (2019), their Appendix D, for details.

5.2. Function Form of the Separation Distribution

We model the separation distribution as a smoothly broken power law:

$$\begin{eqnarray}&&\phi \left(s\right)={\phi }_{0}{\left(\displaystyle \frac{s}{{s}_{{\rm{b}}}}\right)}^{{\gamma }_{1}}{\left\{\displaystyle \frac{1}{2}\left[1+{\left(\displaystyle \frac{s}{{s}_{{\rm{b}}}}\right)}^{1/{\rm{\Lambda }}}\right]\right\}}^{({\gamma }_{2}-{\gamma }_{1}){\rm{\Lambda }}},\end{eqnarray} \tag{ 3 }$$

where ϕ₀ is a normalization parameter. The break separation s_b marks the transition between the two single power laws with index γ₁ and γ₂, respectively. The Λ parameter quantifies how abruptly or smoothly the two power laws are joined. The model has four fit parameters ${\boldsymbol{m}}=({\gamma }_{1},{\gamma }_{2},\mathrm{log}({s}_{{\rm{b}}}/\mathrm{au}),{\rm{\Lambda }})$. We suppose that ϕ(s) is independent of both the distance and the absolute magnitude of the two member stars, although the separation distribution is expected to vary somewhat with the mass (and thus absolute magnitude) of the two stars (e.g., Duchêne & Kraus 2013; Moe & Di Stefano 2017; Moe & Kratter 2018). This supposition is validated with two special data sets selected from the pure samples: one with ΔG > 2.5 mag, which contains 1633 disk-like, 3481 intermediate, and 986 halo-like binaries, and the other with M_G1 < 7 mag, which includes 2287 disk-like, 6684 intermediate, and 3542 halo-like binaries. The analysis on these data sets indicates that the large ΔG or stellar mass just slightly affects the shape of the observed separation distribution at small separations, and this selection effect can be well overcome with the selection function when the intrinsic separation distribution is derived in the following section.

5.3. Likelihood for Fitting the Separation Distribution

Given a sample of binaries with projected separations s_i, the likelihood function is

$$\begin{eqnarray}&&L=p\left(\left\{{s}_{i}\right\}| {\boldsymbol{m}}\right)=\displaystyle \prod _{i}p\left({s}_{i}| {\boldsymbol{m}}\right),\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{m}}$ is a set of free model parameters to be fitted and $p({s}_{i}| {\boldsymbol{m}})$ is the probability of detecting the ith binary given model parameters ${\boldsymbol{m}}$. For the ith binary, $p({s}_{i}| {\boldsymbol{m}})$ can be specified as

$$\begin{eqnarray}&&p\left({s}_{i}| {\boldsymbol{m}}\right)=\displaystyle \frac{\phi \left({s}_{i}| {\boldsymbol{m}}\right)}{{\int }_{{s}_{\min }}^{{s}_{\max }}\phi \left(s| {\boldsymbol{m}}\right){f}_{{\rm{\Delta }}G}\left(s| {d}_{i}\right){ds}},\end{eqnarray} \tag{ 5 }$$

where $\phi \left({s}_{i}| {\boldsymbol{m}}\right)$ is the probability that a binary with distance d_i, magnitude difference ΔG, and physical separation s_i is found in the catalog, given an intrinsic separation distribution (specified by Equation (3)) with parameters ${\boldsymbol{m}}$. The selection function, ${f}_{{\rm{\Delta }}G}\left(s/{d}_{i}\right)={f}_{{\rm{\Delta }}G}\left(\theta \right)$, is given by Equation (2). Here s_min and s_max are the minimum and maximum separations of the observed distribution in the disk and halo binary samples, and $\phi (s| {\boldsymbol{m}})$ is normalized such that ${\int }_{{s}_{\min }}^{{s}_{\max }}\phi \left(s| {\boldsymbol{m}}\right){\rm{d}}s\,=\,1$. Equation (5) does not account for the observational uncertainties in s_i or d_i. These are small because all the binaries in our clean samples have parallax errors smaller than 5%.

The denominator in Equation (5) reflects the fraction of predicted binaries that could have been detected at a distance d_i and magnitude difference ΔG; it accounts for the fact that at large distances and large ΔG, only binaries with large s can be detected. We set s_max = 10^5.25 au (i.e., ∼0.8 pc), omitting the widest separation bin owing to a high contamination rate beyond 10^5.25 au. The choice of s_min has no effect on our results because the integrand in Equation (5) goes to 0 at small separations. We set s_min = 10⁻² au.

If either member star of a binary is too faint, the binary probably cannot be detected. Therefore, whether a binary can be observed or not also depends on the apparent magnitude of both stars. However, this has no effect on the inferred separation distribution as long as the undetected binaries have the same intrinsic separation distribution.

We use emcee (Foreman-Mackey et al. 2013) to sample the posterior distribution of the broken power-law model of binary separations (Equation (3)), for the three subsamples of disk-like, intermediate, and halo-like binaries. We use flat priors for γ₁, γ₂, and $\mathrm{log}({s}_{{\rm{b}}}/\mathrm{au})$ and use an exponential prior on Λ, taking a prior of exp(−3Λ) for Λ. The top panel of Figure 7 displays the contours and marginalized probabilities of the model parameter posteriors ${\boldsymbol{m}}$, while the bottom panel illustrates the inferred intrinsic separation distributions for the three subsamples. The median values of ${\boldsymbol{m}}=({\gamma }_{1},{\gamma }_{2},\mathrm{log}({s}_{{\rm{b}}}/\mathrm{au}),{\rm{\Lambda }})$ are

• 1.  
  $(-{1.51}_{-0.03}^{+0.03},-{2.07}_{-0.19}^{+0.69},{3.97}_{-0.32}^{+0.67},{0.51}_{-0.39}^{+0.59})$ for the disk-like binaries,
• 2.  
  $(-{1.56}_{-0.04}^{+0.03},-{2.84}_{-0.38}^{+0.77},{4.39}_{-0.21}^{+0.34},{0.57}_{-0.41}^{+0.29})$ for the intermediate binaries,
• 3.  
  $(-{1.55}_{-0.06}^{+0.05},-{3.33}_{-0.68}^{+0.72},{4.59}_{-0.29}^{+0.24},{0.67}_{-0.34}^{+0.26})$ for the halo-like binaries.

Zoom In Zoom Out Reset image size

These results, which are illustrated in Figure 7, can be summarized as follows:

• 1.  
  At s < s_b, the three subsamples have the same separation distributions, i.e., a power law with an index of γ₁ ∼ −1.54 shown as the black dashed lines in the bottom panels of Figure 7.
• 2.  
  At s > s_b, the separation distributions become significantly steeper for all three subsamples, excluding single power-law models over the entire separation range (see the bottom panel of Figure 7). This profile steepening (Δγ ≡ γ₂ − γ₂) increases from 0.56 ± 0.50 for the disk-like binaries, through 1.28 ± 0.60 for the intermediate binaries, to 1.78 ± 0.70 for the halo-like binaries. There is thus a strong hint of a steeper fall-off at ultrawide separations (s > 10^4.5 au) for the halo-like sample than the disk sample, although the marginalized constraints on γ₂ are quite broad owing to degeneracies with other parameters.
• 3.  
  The Δγ for the disk-like sample is significant, but small enough to explain why ER18 found comparably good fits for MS/MS binaries with a single power law with a slope −1.6 over 500 < s/au < 50,000. The Δγ = 1.78 ± 0.70 for the halo-like sample, however, indicates a quite substantive change in the separation distribution slope at s ∼ 10^4.3 au (i.e., ∼0.1 pc).
• 4.  
  The smoothing parameter Λ differs among the subsamples: the halo-like subsample has the strongest break (highest Δγ) yet the smoothest power-law transition, i.e., the highest parameter Λ.
• 5.  
  The posterior distributions of the free parameters (except γ₁) gradually present a double-peak feature from the disk-like to the halo-like subsamples. It indicates that the intrinsic separation distributions favor two breaks for all three subsamples, in particular for the old population. One break takes place at s ∼ 10^3.8 au; the corresponding smoothing parameter is Λ ∼ 0.05. The other break is at s ∼ 10^4.5 au, and its Λ ∼ 0.7.
• 6.  
  The three free parameters, i.e., γ₂, $\mathrm{log}({s}_{{\rm{b}}}/\mathrm{au})$, and Λ, are degenerate and correlated with each other. Hence, we will discuss the two breaks with two fixed Λ, i.e., Λ = 0.05 (Figure 8) and Λ = 0.7 (Figure 9), respectively, in the Appendix.

At smaller separations, the observed separation distribution of binaries presumably reflects a combination of birth separations and subsequent evolution (e.g., El-Badry et al. 2019; Liu 2019; Moe & Di Stefano 2017, and references therein). For the very wide binaries (beyond a few thousand au) of primary interest here, the formation process leading to the "initial" separation distribution at the time of birth cluster dispersal and the subsequent evolutionary processes are not firmly established. Wide binaries presumably did not form within a common disk, but resulted from turbulent fragmentation (at closer separations) and stochastic pairing during cluster dissolution (at wider separations). After formation, wide binaries could retain their original orbital parameters, except (1) that they may be affected or become disrupted by gravitational encounters with molecular clouds or other massive objects (e.g., Chanamé & Gould 2004); (2) that either one or both binary members evolve off the MS, inducing rapid mass loss (e.g., ER18); or (3) that dynamical evolution leads to widening within unstable triples (e.g., Reipurth & Mikkola 2012).

This present work was originally initiated with the goal of detecting or constraining signs of binary disruption by massive objects. In such an analysis one must presume to have confident prior knowledge about the initial separation distribution to ∼1 pc, e.g., that it was a perfect single power law at large s. We first pursue this approach of modeling, asking what population of massive scattering objects could lead to the breaks seen. This, however, turns out to have physically implausible implications, at variance with other constraints. In light of this, we then ask what separation distribution we may have expected at birth, and we find it to be consistent with the observations.

7.1. An Initial Single Power Law, Broken by a Population of Scatters?

We start out by noting that the outer slope of the disk-like binary is, if anything, flatter than that of the halo-like populations, showing a steepening by only ≲0.6: this argues against, or at least does not provide evidence for, molecular clouds being a dominant mechanism in setting the outer slope of wide disk-like binaries by means of tidal disruption that would induce a cutoff.

We now use the approach of previous studies, assuming that the separation distribution for wide binaries follows a single (initial) power law. By implication, any break in such a power law is then attributed to external perturbations, e.g., binary disruption as a consequence of encountering MACHOs. Disk-like binaries are more likely to encounter molecular clouds—acting as perturbers—during much of their orbits, and any such encounters will be at lower velocities. It has therefore been inferred (e.g., Chanamé & Gould 2004) that such samples should not be used to probe the properties of MACHOs, as their perturbative effect would be subdominant. However, halo binaries spend only a small fraction of their lifetime in the disk and cross it at high speed, so MACHOs (if they exist) could be the dominant perturbers, with the observed separation distribution providing constraints on their density and mass.

Such external perturbations can be usefully divided into two regimes: a diffusive regime in which perturbations are described by a multitude of individually weak gravitational encounters, and a catastrophic regime in which perturbations are dominated by a single closest encounter. As described in Section 6, we detect steepening at ∼0.1 pc in the separation distribution, which is particularly distinct for the halo-like subsample.

For those two disruption regimes (Binney & Tremaine 2008, Equation 8.65(a) and (b)), the timescales to disrupt a solar-mass binary at a separation of ∼0.1 pc can be expressed as

$$\begin{eqnarray}&&{t}_{d,\mathrm{diff}}\approx 3\,\mathrm{Gyr}\left[\displaystyle \frac{v}{200\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right]\,\left[\displaystyle \frac{30{M}_{\odot }}{M}\right]\,\left[\displaystyle \frac{{\rho }_{H}}{\rho }\right]\,\left[\displaystyle \frac{0.1\,\mathrm{pc}}{s}\right]\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{t}_{d,\mathrm{cat}}\approx 3\,\mathrm{Gyr}\left[\displaystyle \frac{{\rho }_{H}}{\rho }\right]\,{\left[\displaystyle \frac{0.1\mathrm{pc}}{s}\right]}^{3/2}.\end{eqnarray} \tag{ 7 }$$

Here v is the relative velocity of a MACHO perturber passing by one of the binary components, M is the individual mass of the perturbers, and ρ is their mean mass density near the Sun. The local halo mass density is denoted by ρ_H, and we adopted ρ_H = 0.01 M_⊙ pc⁻³. In the catastrophic regime, the disruption timescale at a given separation is independent of the mass and velocity of the perturber, set merely by ρ. If the perturber mass

$$\begin{eqnarray}&&M\gt 10\,{M}_{\odot },\end{eqnarray} \tag{ 8 }$$

t_d,diff would be shorter than the typical time, i.e., 10 Gyr. The condition to see a distinct break due to disruption in the diffusive regime is $\tfrac{30{M}_{\odot }}{M}\cdot \tfrac{{\rho }_{H}}{\rho }\,\lt \,1$. Therefore, M_crit ∼ 30 M_⊙ marks the transition between the two regimes. Previous authors expected to see a signature of binary separation function at s ∼ 0.1 pc in several studies (CG04; Y04; Quinn et al. 2009) but failed to detect the break signature owing to the paucity of observational wide binaries at large separations.

However, the implied lower limits on the MACHO density or MACHO mass implied by the above analysis (critically hinging on the assumption of an initial single power law for the binary separations) would be in direct conflict with the stringent upper limits from the survival of compact ultrafaint dwarf galaxies and the star cluster in Eridanus II. With the latter systems, Brandt (2016) closed the window of allowed MACHO mass in ∼20–100 M_⊙ and thereby ruled out the entire window of MACHO mass in >10⁻⁷ M_⊙ by combining his results with the existing constraints (Y04; Tisserand et al. 2007; Wyrzykowski et al. 2008; Quinn et al. 2009; Green 2016). We also note that the steepening in the observed separation distributions is relatively gradual and is not manifest as a sharp truncation or break. The arguably most obvious way to resolve this tension is to acknowledge that the assumption of an initial single power law was too restrictive and not realized in nature.

7.2. A Simple Model for the Primordial Wide Binary Separation Distribution

We now lay out a very simple model for the separation distribution of wide binaries (wider than a few thousand au) that might be expected as a consequence of their formation. Given their dynamical fragility, we presume that most could have formed only as the cluster (or association) was dispersing. A pair of stars in a dispersing star cluster at separation s will be and presumably remain bound provided that ${v}_{\mathrm{tot}}\lt {v}_{\max }(s,m)\,=\sqrt{2{{Gm}}_{\mathrm{tot}}/s}$, where m_tot  is the total mass of the pair and v_tot is their total relative velocity. For a Gaussian velocity distribution σ (the velocity dispersion of the birth cluster), and for "wide binaries" (v_tot ≪ σ), this will occur with probability

$$\begin{eqnarray}&&p(s\,| \,{m}_{\mathrm{tot}},\sigma )=\displaystyle \frac{4\pi }{3}{\left(\displaystyle \frac{{v}_{\mathrm{tot}}}{\sqrt{8\pi }\sigma }\right)}^{3}\simeq 0.1{G}^{3/2}{\sigma }^{-3}\times {m}_{\mathrm{tot}}^{3/2}{s}^{-3/2},\end{eqnarray} \tag{ 9 }$$

for separations smaller than the cluster's tidal radius, R_tidal. Here we have taken into account the fact that while phase space available to yield s scales ∝s²ds, the s² in this expression is canceled for a typical cluster profile ρ ∝ r⁻², within R_tidal, with no pairs larger than this radius.

Thus, for an ensemble of dissolving clusters with characteristic tidal radii ${R}_{\mathrm{tidal}}\sim { \mathcal O }({R}_{\mathrm{typ}})$, the separation distribution of subsequently bound wide pairs for s ≪ R_typ will always be the same (i.e., s^γ, γ = − 3/2); the sum of such distribution will have the same power law. However, as these clusters will not all have the same R_tidal, the individual-cluster distributions will have breaks at s ∼ R_typ. Hence, for s ≪ R_typ, the joint distribution will be characterized by a true power law, while for s ≳ R_typ it will fall off more rapidly in a manner that presumably could be approximated as a power law, depending on the statistics of R_typ.

This simple model has two attractive features in this context. First, it provides an explanation of why wide binaries in the range ∼1000−10,000 au have a separation distribution of γ ≂ −1.5 across all populations. And it implies that the birth cluster size has an imprint on the (maximal) binary separation distribution even long after the cluster has been dispersed. If stars born in earlier epochs of the Milky Way's evolution (viz., the halo-like population) formed in clusters that were more compact (smaller R_tidal), or with a narrower distribution of R_tidal among the clusters, then their separation distribution should experience a more distinct departure from the γ ≃ −1.5 power law. On the other hand, halo binaries are subject to fewer gravitational perturbations owing to their higher velocities. These two effects work in opposite directions, and it is not immediately clear which will dominate. The fact that the observed separation distributions across the three samples are quite similar even at wide separations suggests that they partially cancel. If the clusters that gave birth at a later epoch to the disk-like binaries had a wider distribution of R_tidal, e.g., because a larger fraction of them were born in barely bound associations, then we would expect a less distinct break in the separation power law. In this model, the present-day separation distribution is still mostly a reflection of the initial distribution. Of course, some tidal disruption may have played a role as s approaches 1 pc. Note also that this model can only infer—not predict—the break radius, s_b, which is larger for the halo-like binaries than for the disk-like ones.

We have compiled an extensive yet pure catalog of candidate wide binaries (a < 1 pc) selected from Gaia DR2 in the solar neighborhood with distances d < 4 kpc, following a procedure similar to that of ER18. This initial candidate catalog consists of 807,611 possible binaries. Its contamination rates are lower than 10% at a < 20,000 au; however, the contamination rates quickly increase beyond 20,000 au, until up to 100% at the largest separation bin, i.e., a ∼ 1.0 pc. To address this, we subsequently applied additional selection criteria, tailored toward three kinematically selected, presumably pure subsamples: disk-like, intermediate, and halo-like binaries. The raw catalog is described in Table 2, which will be released online and available via the PaperData Service of China-VO.

Table 2.  Catalog Description

Column	units	Description
source_id		Gaia source id (int64); star 1
source_id2		Gaia source id (int64); star 2
astrometric_chi2_al		astrometric goodness of fit (χ2) in the along-scan direction; star 1
astrometric_chi2_al_2		astrometric goodness of fit (χ2) in the along-scan direction; star 2
astrometric_n_good_obs_al		number of good CCD transits; star 1
astrometric_n_good_obs_al2		number of good CCD transits; star 2
dec	deg	decl.; star 1
dec2	deg	decl.; star 2
ra	deg	R.A.; star 1
ra2	deg	R.A.; star 2
pairdistance	deg	angular separation between star 1 and star 2
parallax	mas	parallax; star 1
parallax2	mas	parallax; star 2
parallax_over_error		parallax divided by its error; star 1
parallax_over_error2		parallax divided by its error; star 2
phot_bp_mean_flux_over_error		integrated BP mean flux divided by its error; star 1
phot_bp_mean_flux_over_error2		integrated BP mean flux divided by its error; star 2
phot_bp_mean_mag	mag	integrated BP mean magnitude; star 1
phot_bp_mean_mag2	mag	integrated BP mean magnitude; star 2
phot_bp_rp_excess_factor		ratio of total integrated BP and RP flux to G-band flux; star 1
phot_bp_rp_excess_factor2		ratio of total integrated BP and RP flux to G-band flux; star 2
phot_g_mean_flux_over_error		integrated G-band mean flux divided by its error; star 1
phot_g_mean_flux_over_error2		integrated G-band mean flux divided by its error; star 2
phot_g_mean_mag	mag	G-band mean magnitude (Vega scale); star 1
phot_g_mean_mag2	mag	G-band mean magnitude (Vega scale); star 2
phot_rp_mean_flux_over_error		integrated RP mean flux divided by its error; star 1
phot_rp_mean_flux_over_error2		integrated RP mean flux divided by its error; star 2
phot_rp_mean_mag	mag	integrated RP mean magnitude; star 1
phot_rp_mean_mag2	mag	integrated RP mean magnitude; star 2
pmdec	mas yr−1	proper motion in the decl. direction; star 1
pmdec2	mas yr−1	proper motion in the decl. direction; star 2
pmdec_error	mas yr−1	standard error of proper motion in the decl. direction; star 1
pmdec_error2	mas yr−1	standard error of proper motion in the decl. direction; star 2
pmra	mas yr−1	proper motion in R.A. direction, i.e., ${\mu }_{\alpha }^{* }={\mu }_{\alpha }\cos \delta ;$ star 1
pmra2	mas yr−1	proper motion in R.A. direction, i.e., ${\mu }_{\alpha }^{* }={\mu }_{\alpha }\cos \delta ;$ star 2
pmra_error	mas yr−1	standard error of proper motion in R.A. direction; star 1
pmra_error2	mas yr−1	standard error of proper motion in R.A. direction; star 2
radial_velocity	km s−1	spectroscopic barycentric RV; star 1
radial_velocity2	km s−1	spectroscopic barycentric RV; star 1
radial_velocity_error	km s−1	standard error of spectroscopic barycentric RV; star 1
radial_velocity_error2	km s−1	standard error of spectroscopic barycentric RV; star 2
rv_nb_transits		number of transits used to compute RV; star 1
rv_nb_transits2		number of transits used to compute RV; star 2
s_AU	au	projected physical separation between two stars
binary_type		type of wide binary; different types are assigned different integer identifiers: MS–WD, WD–WD, and MS–MS are assigned with 1, 2, and 3, respectively.
num		number of pairs in neighborhood for each binary
U	km s−1	RV in the Cartesian coordinates (positive to the Galactic center); star 1
V	km s−1	azimuthal velocity in the Cartesian coordinates; star 1
W	km s−1	vertical velocity in the Cartesian coordinates; star 1
U2	km s−1	RV in the Cartesian coordinates; star 2
V2	km s−1	azimuthal velocity in the Cartesian coordinates; star 2
W2	km s−1	vertical velocity in the Cartesian coordinates; star 2
U_error	km s−1	error of RV in the Cartesian coordinates; star 1
V_error	km s−1	error of azimuthal velocity in the Cartesian coordinates; star 1
W_error	km s−1	error of vertical velocity in the Cartesian coordinates; star 1
U_error2	km s−1	error of RV in the Cartesian coordinates; star 2
V_error2	km s−1	error of azimuthal velocity in the Cartesian coordinates; star 2
W_error2	km s−1	error of vertical velocity in the Cartesian coordinates; star 2

Note. Each row in the catalog corresponds to a single binary; "star 1" and "star 2" designations in each binary are arbitrary. Full descriptions of Gaia variables can be found at https://gea.esac.esa.int/archive/documentation/GDR2/.

Download table as:  ASCIITypeset image

This raw catalog of candidate wide binaries can provide a starting point for several scientific applications beyond the scope of this paper. For instance, one could select a subsample of wide binaries containing a WD member from the raw catalog to determine stellar ages. The age of a WD can be easily constrained from its cooling age and mass given an initial−final mass relation and an initial mass−age relation, while the age of its companion MS star could not be measured precisely. The two members in a binary usually can be regarded as coeval, so the MS star's age can be simply obtained from its companion WD star (Fouesneau et al. 2019).

To reach low contamination to separations of ∼1 pc, we defined three subsamples of MS–MS binaries from the raw catalog via the systems' projected velocities with respect to the Sun, v_⊥,tot: 4361 disk-like binaries with v_⊥,tot ≤ 40 km s⁻¹, 10,090 kinematically intermediate binaries with 40 km s⁻¹ ≤ v_⊥,tot ≤ 85 km s⁻¹, and 4351 binaries with halo-like kinematics (v_⊥,tot ≥ 85 km s⁻¹). We presume that these velocity cuts represent a rough ordering in the binaries' age. Using three catalogs, i.e., the mock Gaia DR2 catalog (Rybizki et al. 2018), a shifted Gaia DR2 catalog, and around 7.2 million Gaia DR2 stars with RVs, we quantify the contamination rates of the three subsamples after we have applied additional cuts in the different separation bins, and we find that the contamination rates in the three subsamples are negligible in the range of separation s < 1.0 pc.

We then proceeded to model the separation distributions of these subsamples as a smoothly broken power law, with four free parameters ${\boldsymbol{m}}=({\gamma }_{1},{\gamma }_{2},\mathrm{log}({s}_{{\rm{b}}}/\mathrm{au}),{\rm{\Lambda }})$. Fitting these models for the three subsamples for s < 10^5.3 au and accounting for the selection function, several important findings have been discovered: (1) We confirm that the slope for separations from 10^2.5 to 10^4.0 au is p(s) ∝ s^γ = s^−1.54, as found in previous studies; this slope is essentially the same for all three subsamples. (2) We show for the first time at high significance how the slope steepens beyond 10⁴ au ≤ s ≤ 10⁵au, and we find first tentative evidence that this slope steepening differs among these subpopulations: by only Δγ ∼ 0.5 for stars of disk-like kinematics, but by Δγ = 1–3 for stars with halo-like kinematics. The actual values of Δγ would be even larger than these estimates, if we had unrecognized contamination at the largest separations. (3) We also find some, albeit tentative, evidence that s_b increases from 5000 au for the disk-like subsample to 20,000 au for the halo-like subsample.

We have offered interpretations for these observational findings in the two limiting cases: we start by presuming that the birth separation was a single power law to ∼1 pc, altered only by binary disruption effected by some possible MACHO populations. The break in ∼0.1 pc in the separation distribution of the halo-like subsample would imply M > 10 M_⊙ at a very high mean MACHO density. However, this window of MACHO candidates has been closed by other studies, e.g., Brandt (2016), and we must conclude that—in light of complex, but not truncated, separation distributions—MACHO limits or detections would be hard to derive, as no clearly defined null hypothesis for the undisturbed separation distribution exists.

We therefore build a conceptually different simple model asking what we should expect for the initial (birth) separation distribution. We find that a simple model, where binaries form by stars remaining stochastically bound as the cluster disperses, works remarkably well in two respects: it generically predicts a distribution power law of ${s}^{-\tfrac{3}{2}}$, as observed in all populations, and it implies that the separation distribution will steepen beyond separations in excess of the initial cluster size; therefore, (expected) structural differences between the birth clusters or associations of the different populations may now be reflected in the separation distribution of ultrawide binaries.

The authors thank Bing Yue, Qiang Yuan, Chao Liu, and Ling Zhu for helpful discussions. H.-J.T. acknowledges the National Natural Science Foundation of China (NSFC) under grants 11873034, U1731108, and U1731124. K.E.-B. was supported in part by an NSF graduate research fellowship and by SFB 881. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

In order to quantify and clearly illustrate the smaller-separation break in p(s), we sample the posterior distribution of the broken power-law model with a fixed value of Λ = 0.05, setting s_max = 10^4.6 au for each binary subsample. The resulting constraints on the intrinsic separation distributions are displayed in Figure 8, similar to Figure 7. The best-fit values of ${\boldsymbol{m}}=({\gamma }_{1},{\gamma }_{2},\mathrm{log}({s}_{{\rm{b}}}/\mathrm{au}))$ are

• 1.  
  $(-{1.52}_{-0.03}^{+0.03},-{1.74}_{-0.06}^{+0.07},{3.51}_{-0.19}^{+0.16})$ for the disk-like binaries,
• 2.  
  $(-{1.58}_{-0.03}^{+0.02},-{1.90}_{-0.14}^{+0.21},{3.88}_{-0.33}^{+0.18})$ for the intermediate binaries,
• 3.  
  $(-{1.56}_{-0.03}^{+0.03},-{1.90}_{-0.08}^{+0.11},{3.81}_{-0.13}^{+0.10})$ for the halo-like binaries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 is similar to Figure 8, but with Λ = 0.7 and s_max = 10^5.3 au. The best-fit values of ${\boldsymbol{m}}=({\gamma }_{1},{\gamma }_{2},\mathrm{log}({s}_{{\rm{b}}}/\mathrm{au}))$ are

• 1.  
  $(-{1.52}_{-0.02}^{+0.03},-{2.17}_{-0.20}^{+0.41},{4.13}_{-0.31}^{+0.34})$ for the disk-like binaries,
• 2.  
  $(-{1.54}_{-0.03}^{+0.02},-{2.77}_{-0.49}^{+0.45},{4.33}_{-0.33}^{+0.18})$ for the intermediate binaries,
• 3.  
  $(-{1.53}_{-0.04}^{+0.03},-{2.95}_{-0.54}^{+0.61},{4.45}_{-0.30}^{+0.21})$ for the halo-like binaries.

In the two extreme cases, we further confirm that (1) the slope for separations from 10^2.5 to 10^4.0 au is p(s) ≈ s^−1.5; (2) at s > 10^4.0 au, the separation distributions obviously deviate from the single power law; and (3) at s > 10^4.5 au, the slope steepens by a slightly larger amount for the old population than for the young population.

It is worth mentioning that the intermediate subsample presents bimodal parameter constraints, in particular for the parameters of γ₂ and $\mathrm{log}({s}_{{\rm{b}}}/\mathrm{au})$. The value of one peak is consistent with the value in the disk-like population; the value of the other one matches that of the halo-like subsample (see Figure 8, top panel).