Wide-binary Stars Formed in the Turbulent Interstellar Medium

Gaia (Gaia Collaboration et al. 2016, 2019) has revealed a large number of wide binaries with semimajor axes a ≳ 10³ au (e.g., El-Badry & Rix 2018; Igoshev & Perets 2019; Hartman & Lépine 2020; Hwang et al. 2020; Tian et al. 2020; El-Badry et al. 2021). Due to their sensitivity to gravitational perturbations, wide binaries have been used to probe unseen Galactic disk material (Bahcall et al. 1985), dark matter substructure (Peñarrubia et al. 2016), and MAssive Compact Halo Objects (Chanamé & Gould 2004; Quinn et al. 2009; Monroy-Rodríguez & Allen 2014), and constrain the dynamical history of the Galaxy (Allen et al. 2007; Hwang et al. 2022b).

Despite their important astrophysical implications, the origin of wide binaries remains a mystery. Due to their large a's (comparable to the typical size ∼0.1 pc ≈2 × 10⁴ au of molecular cores; Ward-Thompson et al. 2007) and the dynamical disruption in star clusters, it is believed that very wide binaries with a ≳ 0.1 pc can hardly form and survive in dense star-forming regions (Deacon & Kraus 2020). Various wide-binary formation mechanisms have been proposed, including cluster dissolution (Kouwenhoven et al. 2010; Moeckel & Clarke 2011), dynamical unfolding of compact triple systems (Reipurth & Mikkola 2012), formation in adjacent cores with a small relative velocity (Tokovinin 2017), dynamical capture (Rozner & Perets 2023), gas-assisted capture (Rozner et al. 2023), during the formation of ultrafaint dwarf galaxies (Livernois et al. 2023), and in the tidal tails of stellar clusters (Peñarrubia et al. 2016). Observations of the large fraction of wide pairs of young stars in low-density star-forming regions (Tokovinin 2017), the chemical homogeneity between the components of wide binaries (Hawkins et al. 2020), the metallicity dependence of the wide-binary fraction (Hwang et al. 2021), and N-body simulations (e.g., Kroupa & Burkert 2001) suggest that wide-binary components are likely to form in the same star-forming region (Andrews et al. 2019).

In addition, Gaia's high-precision astrometric measurements reveal superthermal eccentricity distribution for wide binaries with separations ≳10³ au (Tokovinin 2020; Hwang et al. 2022c). Hamilton (2022) and Modak & Hamilton (2023) suggest that a superthermal eccentricity distribution cannot be produced dynamically, implying that the initial distribution must itself be superthermal.

The interstellar medium (ISM) and star-forming regions are turbulent, with the turbulent energy mainly coming from supernova explosions (Padoan et al. 2016). The observationally measured power-law spectrum of turbulence spans many orders of magnitude in length scales, from ∼100 pc down to ∼10⁻¹¹ pc in the warm ionized medium (e.g., Armstrong et al. 1995; Chepurnov & Lazarian 2010; Xu & Zhang 2017, 2020) and down to ∼0.01 pc in molecular clouds (MCs; Lazarian 2009; Hennebelle & Falgarone 2012; Yuen et al. 2022). The ubiquitous turbulence plays a fundamental role in the modern star formation paradigm (Elmegreen & Scalo 2004; Mac Low & Klessen 2004; McKee & Ostriker 2007). Our understanding of star formation has significantly improved thanks to the recent development in theories of magnetohydrodynamic (MHD) turbulence (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999; Cho & Lazarian 2002), MHD simulations (e.g., Stone et al. 1998; Federrath et al. 2011; Padoan et al. 2016; Kritsuk et al. 2017), and techniques for measuring interstellar turbulence (e.g., Lazarian & Pogosyan 2000; Heyer & Brunt 2004; Lazarian & Pogosyan 2006; Lazarian et al. 2018; Burkhart 2021; Xu & Hu 2021).

Turbulence not only regulates the dynamics of MCs and the formation of density structures, e.g., filaments, clumps, and cores, where star formation takes place, but also affects the kinematics of molecular gas and dust (Hennebelle & Falgarone 2012), density structures, and young stars over a broad range of length scales. Dense cores in MCs (Qian et al. 2018; Xu 2020) and young stars (Ha et al. 2021; Krolikowski et al. 2021; Zhou et al. 2022; Ha et al. 2022) inherit their velocities from the surrounding turbulent gas, and thus the statistical properties of their velocities are similar to that of the turbulent gas.

In this Letter, we will investigate how the initial turbulent velocities between pairs of stars affect the formation of wide binaries and their eccentricities. We first describe the initial turbulent velocities of pairs of stars in Section 2. In Section 3, we focus on the formation of wide binaries and their eccentricity distribution. More discussion is provided in Section 4. Our conclusions follow in Section 5.

Stars form in high-density structures, i.e., filaments and cores in turbulent MCs, which are generated by the compressions and shocks in highly supersonic turbulence (Federrath et al. 2009; Inoue et al. 2018; Mocz & Burkhart 2018; Xu et al. 2019). Dense cores arise at the collision interfaces of converging turbulent flows, with turbulent motions existing on subcore scales (Volgenau 2004). Naturally, the turbulent velocities of gas are imprinted in those of newly formed stars. As confirmed by recent Gaia observations (e.g., Ha et al. 2021; Zhou et al. 2022; Ha et al. 2022), the velocity differences and spatial separations of young stars statistically follow the power-law velocity scaling of interstellar turbulence.

We consider that in a star-forming region, the initial velocity differences, and spatial separations of stars at birth statistically satisfy the averaged turbulent velocity scaling,

$$\begin{eqnarray}&&{v}_{\mathrm{tur}}=\langle v\rangle ={V}_{L}{\left(\displaystyle \frac{r}{L}\right)}^{\alpha },\end{eqnarray} \tag{ 1 }$$

where r is the initial separation between a pair of stars, and v is the absolute value of their initial velocity difference. The injected turbulent speed V_L and the injection scale L have typical values as ∼10 km s⁻¹ and ∼100 pc for interstellar turbulence (Chamandy & Shukurov 2020). The power-law index α reflects the properties of turbulence, which is typically 1/3 for solenoidal turbulent motions with the Kolmogorov scaling, and close to 1/2 for highly compressive turbulent motions dominated by shocks (Federrath et al. 2009; Kowal & Lazarian 2010). A steeper turbulent velocity scaling (i.e., a larger α) corresponds to a more efficient energy dissipation in shock-dominated turbulence. The Kolmogorov scaling applies to most of the volume of an MC, while the steeper scaling is preferentially seen in small-scale high-density regions that result from shock compressions (Lazarian 2009; Xu 2020; Xu & Hu 2021; Rani et al. 2022; Yuen et al. 2022). For a driven supersonic turbulence, the transition from α = 1/2 to α = 1/3 is seen at the sonic scale, where the transition from supersonic to subsonic turbulent motions takes place (Federrath et al. 2021).

The Gaussianity of turbulent speed distribution f(v) has been extensively studied in the literature (e.g., Townsend & Taylor 1947; Dubinski et al. 1995; Wilczek 2011), which approximately follows (such that ∫f(v)dv = 1)

$$\begin{eqnarray}&&f(v)\approx \sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{1}{{\sigma }_{v}^{3}}\exp \left(-\displaystyle \frac{{v}^{2}}{2{\sigma }_{v}^{2}}\right){v}^{2},\end{eqnarray} \tag{ 2 }$$

where the parameter σ_v is comparable to v_tur (Equation (1)). The initial velocity distribution of stellar pairs and the turbulent velocity distribution of gas are related through the complex process of star formation. They are not necessarily to be identical, but here we will assume that they are equivalent and the initial speed distribution of stellar pairs at a given initial separation has the generic Gaussian form as given in Equation (2).

We note that due to different effects in star-forming regions, e.g., gravitational compression, gravity-driven turbulence (Guerrero-Gamboa & Vázquez-Semadeni 2020; Xu & Lazarian 2020), and protostellar outflows and winds (Mathew & Federrath 2021; Hu et al. 2022), turbulent velocities may not always follow a power-law scaling. The presence of magnetic fields also causes the scale-dependent anisotropy of turbulent velocity distribution (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999). In addition, a detailed modeling on the velocity distribution of wide binaries at birth also requires a better understanding on, e.g., mass accretion, redistribution of angular momentum, which are long-standing problems of (binary) star formation (Larson 2003; Fisher 2004). On the other hand, we are concerned here only with the eccentricity distribution of wide binaries, and in certain limits this distribution will turn out to depend only weakly upon the exact turbulent velocity scaling.

3.1. Formation of Wide Binaries

For a pair of stars to be gravitationally bound at birth, the relative velocity v and r must satisfy

$$\begin{eqnarray}&&v\lt {v}_{\mathrm{bon}}(r)\equiv \sqrt{\displaystyle \frac{2{GM}}{r}},\end{eqnarray} \tag{ 3 }$$

where M is the total mass of the binary. As an illustration, in Figure 1 we present v versus r for pairs of stars formed in turbulent gas over the range 0.01 pc < r < 1 pc (i.e., 2 × 10³ au < r < 2 × 10⁵ au). The colored region corresponds to v < v_bon for M = 2M_⊙. The color scale corresponds to f(v) given by Equation (2). We fix the value of M in this illustration, but note that there can be dependence of M on r (Bate & Bonnell 1997; Lund & Bonnell 2018).

Zoom In Zoom Out Reset image size

The intersection between v_bon and v_tur occurs at

$$\begin{eqnarray}&&{r}_{\mathrm{int}}={\left(\displaystyle \frac{2{{GML}}^{2\alpha }}{{V}_{L}^{2}}\right)}^{\displaystyle \frac{1}{2\alpha +1}},\end{eqnarray} \tag{ 4 }$$

which increases with increasing M and α. For instance, at α = 1/3 we have

$$\begin{eqnarray}&&{r}_{\mathrm{int}}\approx 0.023\,\mathrm{pc}{\left(\displaystyle \frac{M}{{M}_{\odot }}\right)}^{0.6}{\left(\displaystyle \frac{L}{100\,\mathrm{pc}}\right)}^{0.4}{\left(\displaystyle \frac{{V}_{L}}{10\,\mathrm{km}\ {{\rm{s}}}^{-1}}\right)}^{-1.2},\end{eqnarray} \tag{ 5 }$$

and at α = 1/2 we have

$$\begin{eqnarray}&&{r}_{\mathrm{int}}\approx 0.093\,\mathrm{pc}{\left(\displaystyle \frac{M}{{M}_{\odot }}\right)}^{0.5}{\left(\displaystyle \frac{L}{100\,\mathrm{pc}}\right)}^{0.5}{\left(\displaystyle \frac{{V}_{L}}{10\,\mathrm{km}\ {{\rm{s}}}^{-1}}\right)}^{-1}.\end{eqnarray} \tag{ 6 }$$

Obviously, at r < r_int most pairs of stars can be gravitationally bound and form wide binaries. With the increase of v_tur with r, at r > r_int, only a small fraction of pairs can form binaries.

We note that Equation (4) only provides a very rough estimate of the separation scale above which most stellar pairs are unbound. In detailed modeling of wide-binary formation and their population (e.g., Kroupa 1995), the values of the parameters depend on the star formation model and local turbulence properties.

3.2. Eccentricity Distribution

In this section we calculate the eccentricity distribution of wide binaries that results from random pairing of stars.

For a bound two-body system, the specific (per unit reduced mass) energy and angular momentum are defined as

$$\begin{eqnarray}&&E=\displaystyle \frac{1}{2}{v}^{2}-\displaystyle \frac{{GM}}{r},\,\,\,\,\,\,J={rv}\sin \theta ,\end{eqnarray} \tag{ 7 }$$

where θ is the angle between the separation vector r and the relative velocity v , and E and J are related to the eccentricity e by

$$\begin{eqnarray}&&{J}^{2}=(1-{e}^{2})\displaystyle \frac{{G}^{2}{M}^{2}}{(-2E)}.\end{eqnarray} \tag{ 8 }$$

Eliminating E and J from the above expressions leads to

$$\begin{eqnarray}&&\sin \theta =\displaystyle \frac{\sqrt{1-{e}^{2}}}{2u\sqrt{1-{u}^{2}}},\,\mathrm{with}\,u\equiv v{/v}_{\mathrm{bon}}(r).\end{eqnarray} \tag{ 9 }$$

Note that u ∈ (0, 1) for bound binaries. The condition $\sin \theta \in (0,1)$ further constrains the range of u as

$$\begin{eqnarray}&&\sqrt{\displaystyle \frac{1-e}{2}}\lt u\lt \sqrt{\displaystyle \frac{1+e}{2}}\end{eqnarray} \tag{ 10 }$$

at a given e.

Now we imagine that our random pairing process forms an ensemble of randomly oriented binaries of fixed M at a given initial separation r. The resulting number density of these binaries in (u, e) space satisfies (Equation (9))

$$\begin{eqnarray}\begin{array}{rcl}{dN} & \propto & f(v){dvd}\cos \theta \\ & \propto & \displaystyle \frac{f(u){du}}{u\sqrt{1-{u}^{2}}\sqrt{{e}^{2}-{\left(2{u}^{2}-1\right)}^{2}}}{ede},\end{array}\end{eqnarray} \tag{ 11 }$$

where f(u) is the distribution function of u. Therefore, we find the distribution function of e at a given r as

$$\begin{eqnarray}&&p(e)=C\,e{\int }_{\sqrt{(1-e)/2}}^{\sqrt{(1+e)/2}}\displaystyle \frac{f(u){du}}{u\sqrt{1-{u}^{2}}\sqrt{{e}^{2}-{\left(2{u}^{2}-1\right)}^{2}}},\end{eqnarray} \tag{ 12 }$$

where C is a normalization constant, and the integral bounds are given by Equation (10).

As an illustration, Figure 2 shows p(e) (Equation (12)) by assuming that f(u) takes the form of Equation (2) for various values of σ_u = σ_v /v_bon. When σ_u is much less than unity, p(e) has an excess (compared to the thermal distribution p(e) = 2e) at large e's and a deficiency at small e's (see blue line in Figure 2). The deficiency at small e's is caused by the small f(u) near $u=\sqrt{1/2}$ (corresponding to the circular orbit speed).

Zoom In Zoom Out Reset image size

When σ_u is larger than unity, we approximately have f(u) ∝ u² over 0 < u < 1. In this case, Equation (12) simplifies to

$$\begin{eqnarray}&&p(e)\approx \displaystyle \frac{1.06\,e}{\sqrt{1+e}}K\left(\sqrt{\displaystyle \frac{2e}{1+e}}\right),\end{eqnarray} \tag{ 13 }$$

where K is the complete elliptic integral of the first kind (Gradshteyn & Ryzhik 1994). As shown in Figure 2, the above expression corresponds to a superthermal p(e) and agrees well with the cases of σ_u ≳ 1 (purple and green lines).

As an example, in Figure 3 we present p(e) (Equation (12)) by assuming that f(u) takes the form of Equation (2) for M = 2M_⊙, L = 100 pc, V_L = 10 km s⁻¹, and α = 1/3. At r = 0.01 pc with r < r_int (see Equations (5) and (6), Figure 1), it falls in the regime where σ_u is much less than unity. Therefore, we see a deficiency at small e's and a significant excess at large e's. At r = 0.1 pc with r > r_int and σ_u larger than unity, the corresponding p(e) is well described by Equation (13) and is superthermal. We see that irrespective of the value of σ_u , a superthermal p(e) is generally expected. This result does not depend on the exact scaling relation between σ_v and r. In particular, as long as σ_u ≳ 1 and thus f(u) ∝ u², p(e) given by Equation (13) always stands.

Zoom In Zoom Out Reset image size

In Figure 4, we compare our p(e) of wide binaries taken from Figure 3 and that of the wide binaries formed from dynamical unfolding of triple systems at 1 Myr taken from Reipurth & Mikkola (2012). In the latter case, p(e) for stable bound triples declines at large e's. A superthermal p(e) of the wide binaries at binary separations >10³ au in the solar neighborhood is indicated by Gaia observations (Hwang et al. 2022c; see Figure 4). Compared with the observations, we see a more significant excess at large e's that we derive for young wide binaries. This excess is likely to be reduced by other physical processes that may occur after the binary formation, e.g., dynamical interaction/scatterings between binaries and passing stars and MCs. The superthermal p(e) we find may therefore be important for understanding the observations, especially since those observations likely reflect the formation process of wide binaries rather than their subsequent dynamical interactions (Hamilton 2022; Modak & Hamilton 2023).

Zoom In Zoom Out Reset image size

Observationally, wide binaries at separations >10³ au in the solar neighborhood have a superthermal eccentricity distribution (Hwang et al. 2022c). Since the effect of Galactic tides cannot produce the superthermal eccentricity distribution on its own and diffusive scattering with passing stars and MCs would tend to make the distribution more thermal (Hamilton 2022; Modak & Hamilton 2023), the observed eccentricity distribution is likely a relic of the wide binary formation mechanism.

Simulations show that wide binaries formed from cluster dissolution would have a thermal eccentricity distribution (Kouwenhoven et al. 2010) and therefore cannot explain the observation. While wide tertiaries formed from dynamical unfolding of compact triples are predicted to be highly eccentric (Reipurth & Mikkola 2012), Hwang (2023) finds that the eccentricities of wide tertiaries in triples are similar to the wide binaries at the same separations, suggesting that the dynamical unfolding scenario plays a minor role. The remaining formation channels that form eccentric wide binaries are turbulent fragmentation (Bate et al. 1998) and random pairing during star formation phase (Tokovinin 2017). Simulations of turbulent fragmentation suggest that wide binaries at >10³ au are eccentric (e > 0.6) (Bate 2014), although the number of binaries in simulations is too low to have well-characterized eccentricity distributions.

In this paper, we have focused on the eccentricity distribution that results from the random pairing scenario under the consideration that the initial relative velocities of pairs of stars are likely drawn from a characteristic turbulent velocity distribution of the star-forming gas. Given the typical turbulence conditions in star-forming regions, wide binaries formed from random pairing are expected to be highly eccentric, and the eccentricity is left unchanged during stellar evolution with adiabatic mass losses (Silsbee & Tremaine 2017), which may be (partly) responsible for the observed superthermal eccentricity.

We note that while the random pairing scenario may explain the high-eccentricity wide binaries, it is not clear whether it dominates the wide binary formation at >10³ au. Wide binaries formed from random pairing are weakly bound and their further gravitational interactions with other stars may disrupt the binaries (Dabringhausen et al. 2022). However, even if random pairing only contributes a few per cent of wide binaries at >10³ au, they may still significantly change the eccentricity distribution to superthermal (e.g., Hwang et al. 2022a).

Turbulence plays a fundamental role in star formation. The turbulent motions in molecular gas are inherited by stars at their birth. Different from previous studies relying on the long-term dynamical evolution of stars and star clusters, we suggest that the formation of wide binaries with initial separations in the range 0.01 pc ≲ r ≲ 1 pc (i.e., 2 × 10³ au ≲ r ≲ 2 × 10⁵ au) is a natural consequence of star formation in the turbulent interstellar medium. With the velocity differences and separations of newly formed stars statistically following a turbulent velocity scaling as suggested by recent Gaia observations, a pair of stars with a sufficiently small velocity difference at a given separation is gravitationally bound and can form a wide binary.

For a turbulent velocity distribution of pairs of stars, we find that the resulting eccentricity distribution of the bound pairs is generally superthermal. This is true regardless of the exact turbulent velocity scaling and whether the average turbulent velocity v_tur is smaller or larger than the escape velocity v_bon of a gravitationally bound binary. The superthermal p(e) of wide binaries under this formation channel may be important for explaining the observed superthermal p(e) in the solar neighborhood. Comparisons with future measurements on p(e) in star-forming regions will provide further testing of our theory.

S.X. acknowledges the support for this work provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51473.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. H.-C.H. acknowledges the support from the Infosys Membership at the Institute for Advanced Study. This work was supported by a grant from the Simons Foundation (816048, C.H.).

Software: MATLAB (MATLAB 2021).