Wide-binary eccentricity distribution in young star clusters: dependence on the binary separation and mass

We study the wide-binary eccentricity ($e$) distribution in young star clusters and the role of turbulence in setting the form of the $e$ distribution using magnetohydrodynamical (MHD) simulations of star cluster formation. The simulations incorporate gravity, turbulence, magnetic fields, protostellar heating, and jets/outflows. We find that (1) simulations that employ purely compressive turbulence driving produce binaries with a superthermal $e$ distribution ($\alpha>1$ in $p(e) \propto e^\alpha$), while simulations with purely solenoidal driving or natural mixture of driving modes produce subthermal/thermal distributions ($\alpha \leq$ 1), (2) the $e$ distribution over the full range of binary separations in our simulations is set at the early stages of the star cluster formation process, (3) while binaries (separation of $r_{\mathrm{pair}} \leq 1000\, \mathrm{AU}$) have subthermal to thermal $e$ distributions ($\alpha \sim 0.8$), wide binaries ($r_{\mathrm{pair}}>1000\, \mathrm{AU}$) have a superthermal distribution ($\alpha \sim 1.8$), and (4) low-mass binary systems (system masses of $M_{\mathrm{sys}} \leq 0.8\, \mathrm{M_\odot}$) have a highly superthermal distribution ($\alpha \sim 2.4$), whereas high-mass systems ($M_{\mathrm{sys}}>0.8\, \mathrm{M_\odot}$) exhibit a subthermal/thermal distribution ($\alpha \sim 0.8$). The binary eccentricity distribution is often modelled as a thermal distribution. However, our results suggest that the $e$ distribution depends on the range of separation of the sampled binaries, which agrees with the findings from recent Gaia observations. We conclude that the dependence of the $e$ distribution on the binary separation and mass is linked to the binary formation mechanism governed by the turbulent properties of the parent cloud.


INTRODUCTION
The prevalent existence of binary stars places them at the forefront of many fundamental astrophysical problems (e.g., Bahcall et al. 1985;Monroy-Rodríguez & Allen 2014;Peñarrubia et al. 2016;Liu et al. 2019;Andrews et al. 2022;Tauris & van den Heuvel 2023;Offner et al. 2023;Chen et al. 2024).As a key parameter of binary orbital dynamics, the eccentricity  is often modelled to follow a thermal 1 eccentricity distribution () = / = 2 (Jeans 1919).Recent Gaia observations with high-quality parallaxes and proper motions (Gaia Collaboration et al. 2023) have significantly expanded the sample of binaries (El-Badry et al. 2021), including wide binaries with semimajor axes  ≳ 10 3 AU (El-Badry & Rix 2018; Hwang et al. 2022a), which enables detailed statistical studies on their eccentricity distribution.The Gaia observations indicate a non-universal eccentricity distribution among binary stars in the solar neighbourhood, with a transition from a uniform distribution at separations ∼ 10 2 AU to a superthermal 1 distribution at separations greater than 10 3 AU (Hwang et al. 2022a).This shift challenges our traditional understanding of binary eccentricity distributions which are usually modelled as a thermal distribution.The superthermal ec-★ E-mail: sajay.mathew@anu.edu.aucentricity distribution of wide binaries cannot be accounted for by their evolution in the Galaxy under the effects of Galactic tides and encounters with passing stars (Modak & Hamilton 2023;Hamilton & Modak 2023).Hamilton & Modak (2023) suggest that the observed superthermal distribution is more likely a result of an even more superthermal initial distribution at the time of their formation.
The formation mechanisms and statistical properties of binaries, such as eccentricity, in the regime of  ≲ 100 AU, have been extensively studied with hydrodynamic and magnetohydrodynamic simulations in the literature (e.g., Ostriker 1999;Bate 2000;Offner et al. 2010;Ryu et al. 2017;Lee et al. 2019;Rozner et al. 2023;Kuruwita & Haugbølle 2023).For instance, Bate (2014) find that the  distribution is almost uniform in their star cluster formation simulations comprising of a large population of M-dwarfs and later type stars (see also Guszejnov et al. 2023).The median binary separation calculated in the simulations of Bate (2014) is ∼ 20 AU, which agrees with the peak of the separation distribution obtained for M-dwarf systems in observational surveys (Janson et al. 2012;Winters et al. 2019).At similar binary separation ranges, observations also find a uniform eccentricity distribution (Tokovinin & Kiyaeva 2016;Hwang et al. 2022a).With the new finding that the  distribution has a strong dependence on the binary separation (Tokovinin 2020;Hwang et al. 2022a), which in turn, is dependent on the binary mass (Raghavan et al. 2010;Janson et al. 2012;Winters et al. 2019) and the evolutionary state (Tobin et al. 2016(Tobin et al. , 2022;;Kuruwita & Haugbølle 2023) (see Fig. 4), it remains to be extensively studied what mechanisms play dominant roles in the origin of binary eccentricity on different mass and separation scales.In particular, the formation mechanisms and statistical properties of wide binaries with  ≫ 100 AU need further numerical investigation.
Recent Gaia observations have new evidence of the important role of turbulence in star and binary formation.Studies by Ha et al. (2021Ha et al. ( , 2022) ) have observationally demonstrated that young stars born in turbulent interstellar medium inherit turbulent velocities from their parent cloud.They found that the relative velocities of pairs of young stars are well coupled with the turbulent velocities of the surrounding gas, exhibiting a power-law dependence on their separations across scales within 1 − 100 pc.Xu et al. (2023) propose that wide binaries naturally form in turbulent molecular clouds, with the scaling relation between the initial relative velocities and initial separations of binary stars regulated by turbulence.They found this initial condition governed by turbulence naturally results in a superthermal distribution of eccentricities at birth.Using star cluster formation simulations, in this paper, we study the eccentricity distribution of binary star systems and their dependence on the turbulent properties of the parent cloud.
Our primary goal is to numerically probe the possible origin of the observed non-universal () of binaries with  ≳ 100 AU and its dependence on binary separation, system mass, and mass ratio.Previous numerical works focused on the role played by dynamical interaction and gas friction/gas accretion in the formation of closer binaries ( < 100 AU) (e.g., Ostriker 1999;Bate 2000;Ryu et al. 2017;Lee et al. 2019;Rozner et al. 2023).Since the importance of these two formation mechanisms in the regime of close binaries has been already extensively studied in the literature (e.g., Bate 2009Bate , 2012Bate , 2019)), we focus on the third new mechanism of scale-dependent turbulent velocities of gas inherited by wide binaries (≫ 100 AU), which has not been studied with numerical simulations before.We will test the analytical predictions by Xu et al. (2023) regarding the formation of wide binary stars in turbulent environments and the origin of their superthermal ().Specifically, we will also for the first time investigate the role of turbulence with different driving conditions in shaping ().The paper is organized as follows.In Section 2, we detail the simulation data utilised in our study and describe the methodology for calculating the orbital eccentricity of binary systems.In Section 3, we perform statistical analysis on () and its dependence on turbulence driving mechanism, binary separation, system mass, and mass ratio.In Section 4, we delve into the implications of our findings.Section 5 presents the limitations of this study, and in Section 6, we summarise the key results.

METHODOLOGY
To study the properties of young binary stars, we use the star cluster formation simulations from Mathew & Federrath (2021) and Mathew et al. (2023).These simulations are carried out in a 3D computational domain of 2 × 2 × 2 pc 3 dimensions with triply periodic boundaries.An adaptive mesh refinement (AMR) grid structure is used and the maximum effective resolution is 4096 3 cells with a minimum cell size of approximately 100 AU.The star cluster formation is modelled by solving the MHD equations with gravity using a modified version of the flash (version 4) code (Fryxell et al. 2000;Dubey et al. 2008).
The simulations incorporate a vast array of physical mechanisms including gravity, magnetic fields, turbulence, and protostellar feedback in the form of radiative heating and mechanical outflows.Sink particles are used to model protostellar objects.Whenever a spherical region of radius  sink = 250 AU centered on a cell satisfies the criterion for gravitational collapse (Federrath et al. 2010b), the gas within the region is replaced by a sink particle (see discussion on limitations in §5).The position, linear momentum, and angular momentum of the sink particle are determined by the respective quantities of the enclosed gas within the spherical volume of radius  sink (which is set to be the accretion radius of the sink particle).Sub-grid models are used for modelling the protostellar heating (Mathew & Federrath 2020) and outflow feedback (Federrath et al. 2014;Mathew & Federrath 2021).We encourage the reader to refer to Mathew & Federrath (2021) and Mathew et al. (2023) for more detailed description of the simulation model used here.

Turbulence driving
We drive turbulence within the gas in the computational domain by introducing a forcing term in the momentum equation of MHD.The forcing (acceleration) field is modelled using a stochastic Ornstein-Uhlenbeck (OU) process (Eswaran & Pope 1988), which enables us to drive turbulence continuously, with the field varying smoothly in space and time.The turbulence driving module injects kinetic energy on the largest scales, which naturally cascades down to smaller scales, producing a velocity power spectrum ∼  −2 or equivalently a velocity dispersion -size relation of   ∝ ℓ 1/2 (Larson 1981;Ossenkopf & Mac Low 2002;Heyer & Brunt 2004;Roman-Duval et al. 2011;Federrath 2013;Federrath et al. 2021).Our forcing approach allows us to control the fraction of compressive-to-solenoidal modes by performing a Helmholtz decomposition.When the total power in the forcing is contributed by compressive modes, we obtain a purely compressive driving (curl-free), and when the total power is associated with solenoidal modes, we get a purely solenoidal turbulence driving (divergence-free).When no decomposition is carried out, a natural mixture of driving modes emerges (see Federrath et al. 2008Federrath et al. , 2010a, for more details on the turbulence driving method adopted here).A public version of the turbulence driving module is available (Federrath et al. 2022).

Initial conditions
The gas density in the computational domain is initially uniform with  0 = 6.56 × 10 −21 g cm −3 .Since the size of the box  = 2 pc, the total cloud mass  cl =  0 ×  3 ≈ 775 M ⊙ , yielding a mean free-fall time of  ff ≈ 0.82 Myr.Once the turbulence driving module (Federrath et al. 2010a(Federrath et al. , 2022) ) is turned on, over-densities in the form of clumps and filaments are created due to the stirring of the gas and the presence of turbulent shocks.The magnetic field, which is uniform initially with  = 10 −5 G along the z-axis of the computational domain, is modified as a result of the twisting, mixing and compression (Cho & Lazarian 2002), and elongation of magnetic field lines by the turbulent motions (Seta & Federrath 2021), generating a field structure similar to that observed in real molecular clouds (Federrath 2016).
The turbulence driving is first performed without self-gravity for two turbulent crossing times, 2 turb = /(M s ) = 2 Myr (Federrath et al. 2010a), where M is the steady state sonic Mach number and  s = 0.2 km/s is the isothermal sound speed for solar metallicity, molecular gas at 10 K.The velocity dispersion on the turbulence driving scale is set as   =  s M = 1.0 km s −1 such that M = 5 and Alfvén Mach number M A = 2.9.The initial virial parameter is  vir = 2 kin / grav = 0.5, which is in the range of observed values (Falgarone et al. 1992;Kauffmann et al. 2013;Hernandez & Tan 2015).A steady state of turbulence will be reached as the simulation reaches two turbulent crossing times (Federrath et al. 2010a;Price & Federrath 2010), which is when self-gravity is turned on ( = 0).The turbulence driving is maintained during the further evolution.The simulation is then allowed to evolve till a star formation efficiency (SFE) of 5% is reached, i.e., when 5% of the total cloud mass has formed stars.During this evolution phase, the high-density regions (analogous to dense cores) produced by the turbulent shocks become gravitationally unstable and form star clusters (see Mathew & Federrath 2021;Mathew et al. 2023, for further description of the numerical setup).

Multiplicity identification algorithm
To identify the binary pairs in the simulations, we follow the iterative technique used in Bate (2009).In each iteration, we find the closest bound pair in the list of sink particles.This pair is replaced in the list by a binary object with a mass that is equal to the sum of the respective pair.The position and velocity of the binary object are given by the centre-of-mass position and velocity of the pair, respectively.In the next iteration, using the new list, we again search for the closest bound pair and replace them with an object of higher order, i.e., if the pair consists of two single objects, they are replaced by a binary object or if the pair comprises of a single object and binary object, they are replaced by a triple.The process of finding and replacing the closest bound pair is iterated till none of the objects in the set are bound anymore or the only possible combination is a quintuple.Systems of order higher than quadruples are rejected since high-order systems are highly unstable and they can easily decay dynamically.The output of the above algorithm is a list of multiple systems, i.e., a list consisting of singles, binaries, triples, and quadruples.From this list, we can calculate different multiplicity properties of the systems such as the eccentricity and semi-major axis (Offner et al. 2023).This method is already used in Mathew & Federrath (2021) and Mathew et al. (2023) to determine multiplicity properties.The eccentricity  of the binary orbit is given by, (1) Here , ℎ, and  are the specific energy (kinetic+gravitational), specific angular momentum, and the gravitational constant, respectively. 1 and  2 correspond to the mass of the two binary components.Using Eq. 1, we calculate the eccentricity associated with the binaries in our simulations and analyse the eccentricity distribution and its dependence on the turbulent driving mechanism, binary separation, mass and mass ratio, which is presented next.

RESULTS
To study the binary statistics, we utilise a total of 28 simulations of star cluster formation that form 1362 sink particles (stellar objects) overall.The sink particles themselves represent star+disc systems.However, to make discussions lucid, we will refer to sink particles as stars here although they might represent an under-resolved binary since disc fragmentation is not fully resolved (see discussion on limitations in §5).
Using the method described in §2.3, we categorise the stars formed in each of the simulations into multiple systems.The simulations have the same numerical setup and initial conditions as described in §2.2, except that they are characterised by different modes of turbulence driving.In the suite of simulations used here, there are 7 simulations (forming a total of 468 stars) with a purely compressive mode of turbulence driving (curl-free driving), 11 simulations (forming 445 stars) with purely solenoidal turbulence driving mode (divergence-free), and 10 simulations (forming 449 stars) with a natural mixture (equal power in compressive and solenoidal modes).In each of the three sets, all simulations have the same setup but have different turbulence realisations (random seed to generate numerical turbulence).For example, a different turbulence realisation is used in each of the 7 simulations with purely compressive driving, providing a statistically conclusive sample of stars.The number of stars formed depends on the type of turbulence driving (Mathew et al. 2023).Therefore, in order to obtain a similar number of stars (i.e., similar statistics) from each of the three simulation models, the number of simulations performed for each of the models is different.For example, the model with purely compressive turbulence driving produces the highest number of stars in a single simulation, and hence requires the least number of simulations to reach a total of ∼ 450 stars.

Effect of the mode of turbulence driving
The mode of turbulence driving has a significant influence on the density distribution of a gas cloud and in turn on the stellar properties like the star formation rate (Federrath et al. 2008(Federrath et al. , 2010a;;Padoan & Nordlund 2011;Hennebelle & Chabrier 2011) and the initial mass function (IMF) (Schmidt et al. 2010;Lomax et al. 2015;Mathew et al. 2023).To determine which mode of turbulence driving dominates the driving in a cloud, the turbulence driving parameter  can be estimated. is the constant of proportionality in the relation between the standard deviation of gas density and Mach number,  / 0 = M (see Federrath et al. 2008Federrath et al. , 2010a, for a detailed description of ).
Purely compressive driving has  ∼ 1 and purely solenoidal driving has  ∼ 1/3 (Federrath et al. 2008).The  value of real clouds is between 1/3 and 1 (Federrath et al. 2016;Menon et al. 2021;Sharda et al. 2022;Dhawalikar et al. 2022;Gerrard et al. 2023), where a value of  ∼ 0.4 represents a natural mixture of the two driving modes (Federrath et al. 2008(Federrath et al. , 2010a)).Fig. 1 presents the density structure in three simulations that have the same numerical setup and turbulence realisation but have different turbulence driving modes -from left to right: purely compressive (COMP), natural mixture (MIX) and purely solenoidal (SOL).The simulations are compared at a star formation efficiency SFE = 5% (i.e., the ratio of stellar mass to initial cloud mass).The morphology of the star-forming regions in the three simulation models are different and the number of stars formed also varies between them, as discussed in detail in Mathew et al. (2023).Here we primarily focus on the eccentricities of the binary stars.Fig. 2 compares the eccentricity distributions between the three simulation models.The histograms are normalised such that the total area under the histogram is unity.In each of the eccentricity distributions, only the binary pairs that contain the primary star (highest mass star) in the corresponding multiple system are included.For example, in a triple system where a binary subsystem (secondary+tertiary) orbits a primary, only the bound pair that involves the primary is included in the distributions (the secondary+tertiary subsystem is excluded).With such a selection criterion, the COMP model yields a total of 140 pairs, the SOL model gives 132 pairs, and the MIX model yields 117 pairs.The power-law fits to the eccentricity distributions (() ∝   ) are obtained by Monte-Carlo sampling where the red solid line represents the 50 th percentile and the shaded region is bracketed by the 16 th and 84 th .The dashed line corresponds to a thermal distribution with () = 2.The mean eccentricities for the COMP, MIX, and SOL models are  = 0.71 ± 0.02, 0.61 ± 0.02, and 0.63 ± 0.03, respectively (see Tab. 1).The  distribution of the COMP model is superthermal ( ∼ 1.6) while the  distributions for the MIX and SOL models are in the subthermal-thermal range ( ∼ 0.6 and  ∼ 0.8, respectively), and are similar to each other ( values are within the 1-sigma variations of each other).We performed Kolmogorov-Smirnov (K-S) tests to quantify the differences between the distributions.The p-values obtained when comparing the distributions for COMP and MIX is 0.034, and that for COMP and SOL is 0.11.These p-values suggest that () for COMP is likely different from that of MIX or SOL.The p-value of 0.59 obtained from comparing MIX and SOL implies that the  distributions for MIX and SOL are likely similar.This is expected because the turbulence driving parameter  for MIX and SOL is similar (0.4 and 0.3, respectively), i.e., MIX and SOL yield statistically very similar distributions.The variations caused by randomness are likely comparable to that caused by the small change in  between MIX and SOL.Overall, it can be inferred that the  distribution has a considerable dependence on the turbulence driving mode (SOL or MIX on one hand vs. COMP on the other).When compressive turbulence driving dominates, the  distribution tends to be superthermal, while for solenoidal or mixed driving, it corresponds to a subthermal/thermal distribution.
The effect of the turbulent properties of the cloud on the eccentricity statistics will be investigated further in a follow-up work.Here we will mainly focus on the effects of the binary parameters including the binary separation and mass.For the purpose of increasing our sample size, from here on, all analyses will be carried out by compiling together data from the three simulation models, i.e, using the binary pairs from the COMP, MIX, and SOL simulations, yielding a total of 389 binaries.We justify this choice of adding all the models together, in the next paragraph.
The data compilation we obtain will be analogous to binary data obtained from star-forming clouds in the Milky Way since the initial conditions employed in the simulations are typical of low-mass starforming regions in the Milky Way.Similar to the different modes of turbulence driving obtained in the simulations ( ∼ 1,  ∼ 0.4, and  ∼ 0.3), the extent to which a particular driving mode dominates varies between the clouds in the Milky Way (i.e., a range of  values can be associated with the clouds).The binaries in the solar neighbourhood collected by observational surveys possibly originate

SFE = 5%
Figure 3. Eccentricity distribution at SFE = 5% for the three simulation models compiled together.We see that the power-law fit (solid) agrees with a thermal distribution (dashed).
from clouds with different turbulence properties.Since these surveys generally do not distinguish the binary systems based on the initial turbulence environment, combining our simulation models is reasonable.Combining all the simulations also ensures that we have more statistics than the individual cases.
Fig. 3 presents the measured () when the data from the three simulation models are considered as a whole (i.e., plotting the data in the three panels of Fig. 2 together).The resulting collection provides 389 binary pairs.The mean eccentricity is  = 0.65 ± 0.01, denoting the presence of a relatively high fraction of highly eccentric orbits.It is evident from the figure that the fit to the distribution compares well with a thermal distribution of () = 2.The thermal distribution here is simply a result of combining superthermal (COMP) and subthermal/thermal (MIX, SOL) distributions together.Fig. 2 clearly shows that the physics of the turbulence gives rise to different distributions.This suggests that observations might be encouraged to look for non-thermal distributions, which might have been the results of different turbulence in different environments.The distri-bution corresponds to the simulation time at which SFE = 5%.Since the star formation rate differs greatly between the three simulation models, the time at which the SFE = 5% is reached varies between the simulations here.However, as compared to the field stars, the stars here are still very young, with even the oldest stars having a lifetime (time between star formation and end of simulation) of < 1 Myr.We examined the () at SFE = 2, 3, 5%, combining the different simulation models together.The value of  measured for the three groups can be found in Tab. 1.The p-value returned while comparing the  distribution corresponding to SFE = 2% and 3% is 0.23.Comparisons between () at SFE = 2% and SFE = 5%, and () at SFE = 3% and SFE = 5% gives p-values of 0.69 and 0.93, respectively.The distributions are similar and approximately thermal for all SFE, i.e., at different points in time and star formation efficiency, implying that the overall form of the  distribution is decided at the early stages of the cluster formation process.
Observational studies of the eccentricity in binaries generally find mean eccentricity values of < 0.5 (Duquennoy & Mayor 1991;Abt 2006;Raghavan et al. 2010), which is comparatively lower than the value calculated from our simulations ( ∼ 0.65).However, the sample of binaries used in these studies have orbital periods of less than a few 10 5 days, i.e., a semi-major axis of < ∼ 100 AU.Recently, it has been shown that the  distribution is dependent on the separation between pairs or the semi-major axis (e.g., Tokovinin 2020; Hwang et al. 2022a).Hence, the mean estimation will be sensitive to the selection scheme in the survey, i.e., whether the sample has a higher population of wide or close pairs.The binary pairs in our simulations represent a sample with a wide range of semi-major axes, varying from around 10 AU to 10 5 AU, although we are underestimating the number of close binaries.Bate (2014) find a mean eccentricity of 0.33 ± 0.02 with a standard deviation of 0.28 in their simulations, which have comparatively much higher maximum spatial resolution (see also Bate 2009Bate , 2012Bate , 2019;;Guszejnov et al. 2023).The mean eccentricity in Bate (2014) is almost a factor of 2 less than our estimate, because a major fraction of their binary population has separations of around 10−20 AU, while most of the binaries from this work have separations of > 200 AU.The binary statistics, including the eccentricity, are significantly different for different separation or semi-major axis ranges, which is discussed next.

Semi-major axis distribution
Fig. 4 shows the binary semi-major axis formed in the three simulation models combined together, in comparison with the  distributions collected from the radiation hydrodynamical simulations of Bate (2012) (dash-dotted histogram in Fig. 4) and some observational surveys (Raghavan et al. 2010;Winters et al. 2019;Tobin et al. 2022).The peak of the distribution occurs at around 200 AU, although there exists a considerable number of pairs at separations greater than 1000 AU.The number of closer pairs (< 100 AU) is only a lower limit since fragmentation on small scales (at < ∼ 200 AU) is not resolved (see §5 for more details).The criterion for sink particle formation in our simulations does not allow new sink particles to form within the accretion radius of an existing sink particle (250 AU).Therefore, in our simulations, the pairs at separations less than 250 AU are probably the ones that formed at larger separations which later migrated to closer separations due to dynamical interactions (Reipurth & Clarke 2001;Bate et al. 2002b;Goodwin & Kroupa 2005;Offner et al. 2010).Although new sink particles cannot form within the accretion radius of an existing sink particle (in that case the gas is accreted onto the existing sink particle), there is no resolution-dependent restriction on how close two sink particles can come together as the trajectory of sink particles are decided by Lagrangian dynamics, i.e., the gravitational interaction between each other and with the gas.However, the gravitational potential is softened at distances of around 100 AU to prevent artificial collision between sink particles (see §5).Raghavan et al. (2010) find that, in the solar neighbourhood, the distribution of projected separation of binaries with solar-type primaries is Gaussian (on a logarithmic scale) with a peak at 50 AU (see dotted curve in Fig. 4) (see also Duquennoy & Mayor 1991;Mathieu 1994).For M-dwarf primaries in the solar neighbourhood, the distribution is Gaussian as well, but with a slightly lower peak of 20 AU (Janson et al. 2012;Winters et al. 2019).The standard deviation of the distribution in Winters et al. (2019) is  log(/AU ) ≈ 1.16 (see dash-double-dotted curve in Fig. 4).The separation distribution of pre-main sequence (PMS) stars is found to be different from the distribution of the field stars with median separations estimated at a few hundred AU (Mathieu 1994;Duchêne et al. 2004;Connelley et al. 2008;Chen et al. 2013;Tobin et al. 2016).
The VANDAM survey (Tobin et al. 2016) looked at a sample of 94 protostars comprising of Class 0 and Class I objects within Perseus (Enoch et al. 2009) and obtained a bimodal separation distribution with a peak at ∼ 75 AU and another peak at ∼ 1000 AU.Tobin et al. (2022) find that the separation distribution is still bimodal with a peak at ∼ 100 AU and a second peak at ∼ 1000 AU after combining the samples from the Perseus (Enoch et al. 2009;Tobin et al. 2016) and Orion molecular clouds (Fischer et al. 2010;Stutz et al. 2013;Furlan et al. 2016), although the minimum projected separation that can be measured is ∼ 20 AU (see dashed histogram in Fig. 4).Tobin et al. (2022) suggest that the bimodal behaviour is mainly driven by the Class 0 protostars since the peak at ∼ 1000 AU becomes insignificant for Class I and more evolved protostars (see also Kuruwita & Haugbølle 2023).The age of the stars in our simulations lie in the range 10 3 − 10 5 yr which resembles the Class 0 and Class I lifetimes of approximately 10 4 -10 5 yr and a few 10 5 yr, respectively (Froebrich et al. 2006;Dunham et al. 2014).Therefore, the separation distribution in our simulations is more relatable to the distribution for young star-forming regions than that for field stars.However, the comparison of the position of the peak in the distribution with the observationally derived peaks should be done with caution because of the inherent limitations in the numerical resolution.This paper is focused on the relative differences in the binary statistics, particularly the variation in the eccentricity distribution for different separation and mass ranges, i.e., the relative change in the value of .

Dependence on the binary separation
The top row in Fig. 5 depicts the difference in () (at SFE = 5%) between binary pairs with a separation of  pair ≤ 1000 AU (close binaries) and those with  pair > 1000 AU (wide binaries).We see that wide binaries have a higher fraction of highly eccentric orbits.The mean eccentricities for the close and wide binaries are  = 0.63±0.02and  = 0.73 ± 0.02, respectively.It is clear from the fits that the distribution is in the subthermal to thermal range for close binaries ( ∼ 0.8) and superthermal ( ∼ 1.8) for wide binaries, which is consistent with the observational findings by Tokovinin (2020) and Hwang et al. (2022a).The p-value derived from the K-S test on the distributions is 0.013, which means that the distributions are most probably different.Xu et al. (2023) propose that the superthermal nature of wide binaries is an outcome of star formation in the turbulent interstellar medium, where the scaling relations between the velocity differences and the initial separation of stars is regulated by turbulence.Their analytical derivation for () (see Eq. 12 in Xu et al. 2023) agrees comparatively well with our  distribution for wide binaries (see top-right and bottom-right panels in Fig. 5).
The bottom row presents the  distribution when the distinction is based not on the separation but on the semi-major axis .The p-value obtained from the K-S test is 10 −3 , which means that we can reject the hypothesis that the distributions are derived from the same underlying distribution.We see that the trend is similar to that in the top row with  = 0.61 ± 0.02 and  ∼ 0.6 for pairs with  ≤ 1000 AU and  = 0.74 ± 0.02 and  ∼ 1.9 for binaries with  > 1000 AU.It is possible that in a highly eccentric binary that has a Notes.Column 1: MHD simulation model from which binaries are obtained ('ALL' means data from all simulation models are included, i.e., COMP, MIX, and SOL combined).Column 2: star formation efficiency (SFE) at which the calculations are made.Columns 3-5: selected range of separation ( pair ), system mass ( sys ), and mass ratio (), with 'FULL' denoting the full range.Column 6: number of binaries that satisfy the constraints in Columns 1-5.Columns 7 and 8: mean eccentricity () and power-law exponent () in the eccentricity distribution, determined by fitting  () ∝   .
large semi-major axis, the secondary (lower-mass star) may be close to the periastron (closest point) when the separation is measured and hence the separation between the components may be small at that instant of time, i.e., large  but small  pair .However, since the velocity is low close to the apastron (furthest point) and high close to the periastron, at any time instance, it is more probable for the components to be at large separations in the case of wide pairs (large ).Hence the corresponding distributions in the top and bottom rows are statistically similar, i.e., the separation and semi-major axis are analogous when analysing the  distribution (see also Dupuy & Liu 2011;Hwang et al. 2022b).
The mean eccentricity is found to be increasing on increasing the separation range by an order of magnitude with  = 0.62 ± 0.02 in the range 1.5 ≤ log( pair /AU) < 2.5,  = 0.67 ± 0.02 in the range 2.5 ≤ log( pair /AU) < 3.5, and  = 0.76 ± 0.03 in the range 3.5 ≤ log( pair /AU) < 4.5.The distributions also transition from subthermal ( ∼ 0.7) in 1.5 ≤ log( pair /AU) < 2.5 to a thermal distribution ( ∼ 1.2) in 2.5 ≤ log( pair /AU) < 3.5, and then to a superthermal distribution ( ∼ 2.5) in 3.5 ≤ log( pair /AU) < 4.5 (see Tab. 1).This shows that the type of distribution can be different even when the separation ranges of interest differ by just an order of magnitude.
Tokovinin (2020) finds that the  distribution at 200 − 1000 AU separations is nearly thermal and at < 200 AU, the distribution has a clear shortage of highly eccentric pairs.They find that, at > 10 3 AU separations, the distribution is slightly superthermal.Hwang et al. (2022a) find that the () is uniform at ∼ 100 AU, thermal at ∼ 10 2.5 − 10 3 AU and superthermal at > 10 3 AU.Our findings are consistent with these observations.We remind the reader that the actual lower error limits in the mean calculations and  estimates might be even lower.This is because we are probably underestimating the number of low eccentricities due to the inherent numerical limitations (see §5).The underestimation is not taken into account while calculating error limits because an extensive resolution study is needed to quantify the extent of underestimation, which would require a separate study altogether.The underestimation becomes particularly significant for close binaries, e.g., in the above discussed binary range 1.5 < log( pair /AU) < 2.5.
Fig. 6 compares the value of  obtained in the three separation ranges (1.5 ≤ log( pair /AU) < 2.5, 2.5 ≤ log( pair /AU) < 3.5, and 3.5 ≤ log( pair /AU) < 4.5) with the measurements of  as a function of the binary separation by Hwang et al. (2022a) using the Gaia data.The  estimates in the ranges 1.5 ≤ log( pair /AU) < 2.5 and 2.5 ≤ log( pair /AU) < 3.5 agree well with the measurements of Hwang et al. (2022a).However, our simulations produce a higher value of  in the range 3.5 ≤ log( pair /AU) < 4.5 as compared to The  distribution at SFE = 5% where the semi-major axis of the bound pairs is  ≤ 1000 AU (left) and  > 1000 AU (right).For the small separations, the distribution is sub-thermal, while the distribution is super-thermal for the large separations.This is consistent with observational findings by Tokovinin (2020) and Hwang et al. (2022a).The dash-dotted line represents the analytically derived  () for wide binaries by Xu et al. (2023).that for the wide binaries with separations of > 10 3.5 AU in Hwang et al. (2022a).The excess of highly eccentric orbits is most likely because the stars in our simulations represent a very young cluster as compared to the ones probed by the Gaia observations.The value of  for the wide binaries is likely to be reduced with further evolution in the Galactic field as a result of physical processes like dynamical scatterings between binaries and passing stars and molecular clouds (Modak & Hamilton 2023;Hamilton & Modak 2023).

Dependence on the binary mass and mass ratio
Fig. 7 shows () for low-mass binary systems ( sys ≤ 0.8 M ⊙ , left panel) and high-mass systems ( sys > 0.8 M ⊙ , right panel), where  sys =  1 +  2 .A large fraction of the low-mass binaries have highly eccentric orbits ( > 0.9) with  = 0.70 ± 0.03.The p-value obtained from the K-S test on the distributions is 0.022, and therefore we can say that the two distributions are different with 98% confidence.Further, the low-mass systems have a highly superthermal distribution ( ∼ 2.4) while the distribution for the high-mass systems is almost thermal ( ∼ 0.8,  = 0.64 ± 0.01), i.e., () is dependent on the total mass of the binary.This suggests that the form of the  distribution obtained in observations would vary depending on the mass range of binaries selected in the survey (see also Dupuy & Liu 2011).
We compared () for binaries with a mass ratio ( =  2 / 1 , where  2 <  1 ) of  ≤ 0.5 and those with  > 0.5 and find that the distributions for both groups are similar (p-value = 0.96) The  distribution at SFE = 5% for binary systems with  pair ≤ 500 AU and  value given by  ≤ 0.5 (left) and  > 0.5 (right).Bottom: The  distribution at SFE = 5% for binaries with  pair > 500 AU and  value given by  ≤ 0.5 (left) and  > 0.5 (right).and represent thermal distributions, suggesting that the form of the  distribution is generally independent of the  value (see Tab. 1).However, when looked at on different separation scales, the  value does seem to influence ().The top row in Fig. 8 compares the () for binaries with  ≤ 0.5 (left panel) and  > 0.5 (right panel) in the separation range  pair ≤ 500 AU.The distribution is thermal ( ∼ 0.9,  = 0.64 ± 0.02) in the former and subthermal in the latter ( ∼ 0.6,  = 0.61 ± 0.02).However, a K-S test returns a p-value of 0.39, which means it is possible that the distributions are similar.The opposite trend can be seen in the bottom panel, where the same comparison is made, but for binaries with  pair > 500 AU.For  ≤ 0.5 (left panel), we have a thermal distribution ( ∼ 1.1,  = 0.67 ± 0.03), while for  > 0.5 (right panel), the distribution is highly superthermal ( ∼ 2.7,  = 0.75±0.03).The p-value obtained from the K-S test is 0.097, i.e., we can say with 90% confidence that the distributions are different.The features seen in Fig. 8 suggest again that the properties related to the eccentricity distribution are dependent on the separation range.The superthermal distribution obtained for the wide binaries with high mass ratios agrees with the observational finding in Hwang et al. (2022b).Fig. 9 shows a scatter plot of the system mass  sys as a function of the semi-major axis, .The markers are colour-coded based on the time of formation (in units of the total simulation time) of the primary star in the pair.It is clear that the pairs with high system masses form in the early stages of cluster formation.These pairs also fall under a wide range of semi-major axis values ranging from around 10 AU to a few 10 5 AU.The low-mass systems on the other hand form towards the last stages of the cluster formation and they lie mostly in the 100 − 1000 AU range.This, along with the finding that low-mass systems have an  distribution that is different to that of high-mass systems (see Fig. 7), suggests that different mechanisms may be responsible for binary formation on different mass scales, which is discussed next.

Eccentricity distribution and binary formation mechanism
In the following, we discuss several binary formation mechanisms, ranging from large-scale (cloud) to small-scale (disc) fragmentation.

Fragmentation on cloud scales (fragmentation by turbulence-induced density fluctuations)
Star formation is typically initiated by the turbulent shocks prevalent in molecular clouds.These shocks induce overdensities in the gas cloud, which, when they become gravitationally unstable, lead to the birth of stars (Robertson & Goldreich 2018;Mocz & Burkhart 2018).
When a pair of stars (originating from two individually collapsing overdensities) are gravitationally bound at birth, they form as a binary.Due to the log-normal distribution of gas density in turbulent environments (e.g., Vazquez-Semadeni 1994;Padoan et al. 1997;Kowal et al. 2007;Kritsuk et al. 2007;Federrath et al. 2008;Federrath 2013;Hopkins 2013b;Federrath & Banerjee 2015;Kuffmeier et al. 2019;Seta & Federrath 2022), these overdensities vary significantly in size and mass.This variation gives rise to diverse stellar formations, including binary systems with different separations and system masses.Since the overdensities are produced by supersonic turbulent shocks, the binaries formed by pairing-up of individually collapsing overdensities will generally have separations comparable to the supersonic and transonic scales, i.e., ≳ 0.1 pc (Federrath et al. 2021).
Different mechanisms for the formation of wide binaries have been proposed in the literature (e.g., Kouwenhoven et al. 2010;Moeckel & Clarke 2011;Reipurth & Mikkola 2012;Peñarrubia et al. 2016;Tokovinin 2017;Livernois et al. 2023;Rozner & Perets 2023;Rozner et al. 2023).Xu et al. (2023) clearly addressed the role of turbulence and analytically derived that wide binaries (∼ 0.01 − 1 pc) formed from turbulence-induced density fluctuations have superthermal  distributions by considering that the initial velocity difference and separation of the binary components follow the turbulent gas velocity scaling.The superthermal () obtained for the high  pair binaries in our simulations confirms their theoretical prediction (see Fig. 5).We also find that the form of the  distribution obtained for the simulations that employ purely compressive turbulence driving is different from that of the simulations with purely solenoidal driving or a natural mixture of turbulence driving modes.This indicates that the turbulent properties of the parent cloud play a critical role in the binary eccentricity statistics.
In addition, protostellar outflow feedback promotes fragmentation by driving core-scale turbulence (Federrath et al. 2014;Mathew & Federrath 2021;Guszejnov et al. 2021;Hu et al. 2022).The fragmentation on core scales would produce binaries with separations in the range of ∼ 10 3 AU to 0.1 pc (Offner et al. 2023).On smaller scales, the turbulence transitions from transonic to subsonic (Federrath et al. 2021), where gravity, thermal effects, magnetic fields, and rotation begin to dominate over turbulence.The formation of systems with very low masses and with close separations (≲ 500 AU) at the time of their formation by turbulent fragmentation is relatively uncommon, as it is difficult to generate sufficiently strong overdensities through turbulent shocks on core scales, as those scales are intrinsically transto subsonic, at best very midly supersonic (Elmegreen & Elmegreen 1978;Lubow & Pringle 1993;Clarke 1999;Bonnell et al. 2008;Offner et al. 2023).
The low-mass systems formed in the present simulations, which represent a relatively small population, have highly superthermal  distributions (see Fig. 7).These systems are found to have formed during the late stages of the cluster formation and have separations predominantly in the range 100 − 1000 AU (see Fig. 9).These characteristics are unlike those of higher-mass systems, which form relatively early (likely via turbulence-induced density fluctuations in shocks; see above), and span wider separations of up to a few 10 5 AU (see Fig. 9).The variation in the characteristics of () can be explained by the significant changes in the environment in which the binaries form, at early vs. late stages or on large vs. small scales.As stars continue to form in the shocked regions of gas, the local density increases and the local virial parameter decreases (more unstable gas) in cluster-forming regions due to the concentration of gas there.At late stages, the increase in the density and the increased influence of gravity allow the small collapsing overdensities to develop substructures with low Jeans length and low Jeans mass (Bonnell et al. 2008).The binaries that form by the gravitational fragmentation in these collapsing overdensities typically posses lower system masses and shorter separations compared to their counterparts that formed earlier through turbulent fragmentation.Numerical works (Bonnell et al. 2008;Mathew & Federrath 2021), including the simulations presented here (see Fig. 10 in Mathew & Federrath 2021), have shown the existence of a negative correlation between the velocity at the time of formation and the final stellar mass, i.e., the low-mass (late-forming) stars have relatively high velocities as compared to the high-mass (early-forming) stars.This agrees with the argument that the low-mass stars that form in the late stages fragment out of the infalling gas (Bonnell et al. 2008).Such a negative correlation was also obtained in the recent observations of Wei et al. (2023) (see Fig. 10 & 11 in Wei et al. 2023).The superthermal nature of the binaries formed from gravitational fragmentation is possibly due to the preferential stretching of the binary orbit by the gravitational pull in the direction of the nearby sub-cluster.However, a detailed study would be required to test this scenario.

Fragmentation on disc scales
Fragmentation in extended discs (Bate et al. 2002a;Goodwin & Whitworth 2007;Stamatellos et al. 2007Stamatellos et al. , 2011;;Rogers & Wadsley 2012;Thies et al. 2015) is likely to happen in our simulations.However, disc fragmentation on typical scales of ≲ 200 AU (see reviews by Kratter & Lodato 2016;Lee et al. 2020, and the references therein) does not occur here, as these scales are not sufficiently resolved in our simulations (see §5).We expect disc fragmentation (e.g., Adams et al. 1989;Shu et al. 1990;Bonnell & Bate 1994;Stamatellos & Whitworth 2009;Stamatellos et al. 2011;Kratter & Lodato 2016) to be an important binary formation mechanism (see also Hennebelle et al. 2019) on these scales, favouring a uniform  distribution (Raghavan et al. 2010;Duchêne & Kraus 2013;Moe & Di Stefano 2017;Hwang et al. 2022a;Ceppi et al. 2024), where very close binaries become circularised due to tidal effects.Thus, with higher numerical resolution, the  value in the short separation ranges is expected to be somewhat lower than suggested based on our current simulation results.However, since these simulations incorporate a large array of physical mechanisms and we carry out multiple such simulations to produce statistically significant data, achieving higher resolutions is not feasible at the moment, due to the associated computational cost, but remains a high priority for future work.

Dynamical interactions and gas friction/gas accretion
Stars generally form in highly clustered environments (Lada & Lada 2003), and therefore dynamical interactions are a natural outcome of the star formation process.A considerable number of binaries undergo dynamical decay due to interactions with the gas or with other systems and migrate to closer separations in ∼ 100 kyr (Reipurth & Clarke 2001;Bate et al. 2002b;Goodwin & Kroupa 2005;Offner et al. 2016;Lee et al. 2019).This why our simulations have a significant number of binaries with separations in the range ∼ 10 − 200 AU, even though fragmentation on these scales is not resolved in our simulations.A recent study by Kuruwita & Haugbølle (2023) finds that core fragmentation and dynamical capture can produce a considerable number of low-mass close binaries (with disc-scale separations) via efficient in-spiral.During the migration, the eccentricity can change to a different value due to the dynamic nature of the orbital decay.For instance, the hydrodynamic drag forces cause eccentric orbits to become more circular (Szölgyén et al. 2022).Similarly, initially close binaries can be widened due to interactions with other multiple systems and some of the initially wide binaries can even become unbound, since they are generally more weakly bound than close binaries.

Future prospects
An initial thermal eccentricity distribution (Jeans 1919) is frequently adopted in theoretical and numerical modelling of binary populations and star clusters (Kroupa 1995;Dabringhausen et al. 2021).However, adhering to a purely thermal eccentricity distribution can lead to inaccuracies in predicting the evolution of binary populations and their merger rates (Geller et al. 2019).Our study provides physically justified and numerically tested initial eccentricity distributions of wide binaries for future studies.These distributions account for the complexities and variations introduced by different MHD turbulence properties in the surrounding medium.The superthermal eccentricity distribution of wide outer binaries in triple stellar systems also has important implications for the formation channel of black hole binary mergers (Su et al. 2024).Su et al. (2024) find that the outer eccentricity distribution can remain significantly superthermal for modestly hierarchical black hole triples.
The turbulence in our current simulations is weakly magnetized with M A ≈ 3. Observations suggest variations of M A values in different molecular clouds (Hu et al. 2019).Earlier studies show that strong magnetic fields can play an important role in the formation of binaries, especially for the formation of massive close binaries (e.g., Lund & Bonnell 2018;Harada et al. 2021).We will investigate the effect of magnetic fields on shaping the eccentricity distributions of binaries in our future work.

Numerical resolution
Due to the limitations in the numerical grid resolution, the accretion discs are not fully resolved in our simulations.Further, to prevent artificial fragmentation due to resolution limitations, the sink particle formation criteria (Federrath et al. 2010b) in our simulations does not allow the formation of a new sink particle within the accretion radius of an existing one (250 AU).Hence, the presence of close binaries (≲ 200 AU) in our simulations cannot be directly linked to fragmentation on those scales, e.g., disc fragmentation.These binaries are the ones that formed at relatively large separations initially and then decayed to small separations (Reipurth & Clarke 2001;Bate et al. 2002b;Goodwin & Kroupa 2005;Offner et al. 2010).
The resulting effect of the resolution limitation is that, although stars can form in extended discs in our simulations (Bate et al. 2002a;Goodwin & Whitworth 2007;Stamatellos et al. 2007Stamatellos et al. , 2011;;Rogers & Wadsley 2012;Thies et al. 2015), binary formation due to fragmentation on typical disc scales (∼ 100 AU), particularly the formation of spectroscopic binaries (period of a few days), is underestimated.Spectroscopic binaries are found to have eccentricities close to zero (Abt 2006) as a result of the tidal circularisation (Zahn 1977;Tassoul & Tassoul 1992).Hence the fraction of pairs with low eccentricities is underestimated in our simulations.

Gravitational softening
When the separation between two star particles becomes very small, the gravitational acceleration becomes exceptionally high.As the distance approaches zero, the acceleration will go to infinity and the simulation timestep will tend to zero, which will stall the simulation.To overcome such a problem and also because the sink particle accretion radius is limited to the gas resolution, gravitational softening is introduced for distances shorter than a given softening radius (here equal to the sink particle accretion radius of 250 AU).We utilise the version of spline softening (Price & Monaghan 2007) used in Federrath et al. (2010b) to soften the gravitational interaction.The interaction is unaffected at distances greater than the softening radius, while at shorter distances, the acceleration smoothly approaches zero.Such a softening scheme may result in some of the orbits having artificial eccentricities if the binary spends time at separations shorter than the softening radius.The underestimation in the acceleration will be around a factor of 2 at a distance of half of the softening radius (see Fig. 1 in Federrath et al. 2010b).
The limitations in numerical resolution and the application of gravitational softening imply that the actual values of  are affected.However, the main objective of this study is to understand the relative change in , for example, when considering binaries in different separation and mass ranges.Therefore, at least the relative inferences from this paper still hold.

CONCLUSIONS
We use the star cluster formation simulations of Mathew & Federrath (2021) and Mathew et al. (2023) to investigate binary statistics, particularly the eccentricity distribution and its dependence on the binary separation, mass, and mass ratio.The simulations employ gravity, turbulence, magnetic fields, and protostellar feedback in the form of radiative heating and jets/outflows, and therefore include most of the main relevant physical ingredients for star and binary formation.The dataset comprises a total of 28 simulations, which together produce 1362 stars, of which 389 are binaries, allowing us to make statistically conclusive remarks.We find that (1) the eccentricity distribution () is dependent on the mode of turbulence driving (see Fig. 2).The simulations with purely compressive driving produce  distributions with the form () ∝   with  > 1, while simulations with purely solenoidal and a natural mixture of driving modes have  < 1. Kolmogorov-Smirnov tests performed on the three simulation sets suggest that (with ≳ 90% confidence), the  distribution obtained from the purely compressive driving simulations is different from that of simulations with purely solenoidal or a natural mixture of driving modes.This suggests that the turbulent properties of the cloud could play a significant role in shaping the  distribution.We conclude that this is because the turbulence regulates the cloud dynamics including the density and velocity statistics of the gas from which the binaries emerge, and therefore influences their orbital parameters as well.
(2) () is thermal ( ∼ 1) when data from all the simulations are compiled together (see Fig. 3).Recent observational surveys also find that the  distribution is thermal for broad separation ranges comparable to that in our simulations (10 − 10 5 AU).It is likely that these surveys include binary stars that formed in different turbulent conditions and environments, including different turbulence driving modes.It is also worth noting that the binary systems included in the surveys generally have undergone significant dynamical evolution in comparison to the young systems in our simulations, which tends to drive the eccentricity distribution towards thermal (e.g., Hamilton & Modak 2023).
(3) () is similar with  ∼ 1 at different points in time in our simulations where the cluster formation process is followed up to around 1 Myr (see Tab. 1).This shows that the overall form of the  distribution is imprinted in the early stages of the formation process, and small-scale dynamical interactions (within a cluster) do not significantly influence the  value.
(4) () is dependent on the binary separation or the semi-major axis, where wide binaries ( pair > 1000 AU or  > 1000 AU) have superthermal distributions with  ∼ 2 as compared to binaries with  pair ≤ 1000 AU or  ≤ 1000 AU, which have subthermal distributions with  < 1 (see Fig. 5).K-S tests suggest that the hypothesis that the two distributions are the same can be ruled out.The value of  transitions from subthermal to thermal, to superthermal, with increasing binary separation (see Tab. 1), which concurs with the recent observational findings of e.g., Tokovinin (2020), Hwang et al. (2022a) from Gaia (see Fig. 6).We also find that the birth eccentricity distribution of wide binaries is more superthermal than the observed one in the Galactic field.
(5) () also depends on the system mass of the binary with lowmass systems ( sys ≤ 0.8  ⊙ ) producing highly superthermal distributions ( ∼ 2.4), and higher-mass systems ( sys > 0.8  ⊙ ) producing subthermal to thermal  distributions ( ∼ 0.8) (see Fig. 7).Based on a K-S test, we can say with 98% confidence that the two distributions are different.
(6) () is independent of the binary mass ratio () when looked at in the whole separation range (see Tab. 1).However, for relatively close binaries (here  pair ≤ 500 AU), the distribution is thermal ( ∼ 0.9) when  ≤ 0.5 and subthermal ( ∼ 0.6) when  > 0.5 (see Fig. 8).However, a K-S test suggests that there is a 40% chance that the two distributions are similar.In the separation range  pair > 500 AU, the distribution is thermal ( ∼ 1.1) when  ≤ 0.5 and highly superthermal ( ∼ 2.7) when  > 0.5, i.e., an opposite trend in the change of  as compared to that for close binaries ( pair ≤ 500 AU).We find that the above distributions are different with 90% confidence based on a K-S test.Our finding is consistent with the observational finding on the highly superthermal eccentricity distribution of wide twin binaries (Hwang et al. 2022b).
Our study suggests that the often adopted thermal eccentricity distribution (Jeans 1919) may not always be valid or appropriate, and properties of the cloud in which stars are formed have a direct influence on the form of the eccentricity distribution (see Tab. 1).The dependence of the eccentricity distribution on the binary separation and mass indicates that different mechanisms are responsible for binary formation on different scales (see also Bate 2009Bate , 2012Bate , 2014Bate , 2019;;Guszejnov et al. 2023).Our results suggest that wide binaries ( pair ≳ 1000 AU) predominantly result from turbulenceinduced density fluctuations (Xu et al. 2023) and turbulent fragmentation in cores (Padoan & Nordlund 2002;Krumholz & McKee 2005;Hennebelle & Chabrier 2008, 2009;Offner et al. 2010;Hopkins 2012;Bate 2012;Guszejnov & Hopkins 2015;Kuffmeier et al. 2019;Offner et al. 2023).Our findings are consistent with the theoretical predictions of Xu et al. (2023) that the superthermal nature of the eccentricity distribution of wide binaries is an outcome of the turbulent characteristics of the cloud in which they are born.In the intermediate separation range (100 <  pair (AU) < 1000), turbulent fragmentation, gravitational fragmentation (Bonnell et al. (2008); see also Lee & Hennebelle (2018)) and fragmentation in extended discs (Bate et al. 2002a;Goodwin & Whitworth 2007;Stamatellos et al. 2007Stamatellos et al. , 2011;;Rogers & Wadsley 2012;Thies et al. 2015) may have a greater influence on the eccentricity distribution.
It is to be noted that, due to the highly dynamic nature of the star formation process, close binaries can form from initially wider binaries.The decrease of their separations can be caused by e.g., gas dynamical friction (Ostriker 1999;Lee et al. 2019;Rozner et al. 2023;Rozner & Perets 2024), gas accretion (Bate 2000), and dynamical interactions in unstable multiple systems (Reipurth & Clarke 2001;Bate et al. 2002b;Goodwin & Kroupa 2005;Offner et al. 2010).Similarly, dynamical interactions can also widen the binary orbit (Hwang et al. 2022b).On scales less than a few hundred AU, disc fragmentation (e.g., Adams et al. 1989;Shu et al. 1990;Bonnell & Bate 1994;Stamatellos & Whitworth 2009;Stamatellos et al. 2011;Kratter & Lodato 2016) may be a dominant mechanism leading to a uniform eccentricity distribution (Raghavan et al. 2010;Duchêne & Kraus 2013;Moe & Di Stefano 2017;Hwang et al. 2022a;Ceppi et al. 2024).However, disc fragmentation is not resolved in our simulations and hence on those scales, only binaries that formed from turbulent fragmentation, which then dynamically decay to short separations exist.
We conclude that the properties of binaries, in particular their eccentricity, separation, total mass, and mass ratio, depend on the turbulence properties of the parent cloud from which they form.The statistical studies of binaries in young clusters provide valuable constraints on binary formation mechanisms and physically justified initial conditions for binary population synthesis.

Figure 1 .
Figure 1.Column density maps (mass-weighted) of the purely compressive driving (COMP; left panel), mixed driving (MIX; middle panel), and solenoidal driving (SOL; right panel) simulations, at a star formation efficiency (SFE) of 5%.The circular markers correspond to the sink particle (star+disc system) positions, and the colour bar on the right represents the mass of the sink particles.The size of the markers is also scaled by the mass of the sink particles.

−−−Figure 2 .N
Figure2.Eccentricity distribution of the bound pairs at SFE = 5% for each of the three simulation models (from left to right: COMP, MIX, SOL).The solid curves represent power-law fits (  () ∝   ) to the distributions, and the dashed curves correspond to a thermal eccentricity distribution  () ∝  1 .The mean eccentricity is also indicated in the legends, which is higher for the COMP case than for MIX and SOL cases, which have similar values.The COMP case shows superthermal ( > 1) eccentricity distribution whereas the MIX and SOL cases show subthermal to thermal ( ≤ 1) distributions.This demonstrates that the eccentricity distribution of binary stars depends on the turbulent driving mode of the parent cloud in which they are formed.

Figure 4 .
Figure 4.The distribution of the semi-major axis  of the bound pairs from our simulations (histogram with solid edges).The solid line corresponds to a Gaussian fit with a peak at  = 206 AU and a standard deviation  log(/AU) = 0.98.The dash-dotted histogram (adapted from Bate 2012) represents the  distribution (including the pairs in binary, triple, and quadruple systems) in the radiation hydrodynamical simulations ofBate (2012).The dotted curve represents  () for solar-type binaries obtained byRaghavan et al. (2010).The dash-double-dotted curve represents  () for M-dwarf systems derived in the survey byWinters et al. (2019).The dashed histogram (adapted fromTobin et al. 2022) corresponds to the bimodal  distribution obtained for a binary sample in the Orion and Perseus molecular clouds combined together (∼ 400 protostars), comprising mainly of Class 0 and Class I objects.

Figure 5 .
Figure 5. Top: The  distribution at SFE = 5% where the separation between the bound pairs is  pair ≤ 1000 AU (left) and  pair > 1000 AU (right).Bottom:The  distribution at SFE = 5% where the semi-major axis of the bound pairs is  ≤ 1000 AU (left) and  > 1000 AU (right).For the small separations, the distribution is sub-thermal, while the distribution is super-thermal for the large separations.This is consistent with observational findings byTokovinin (2020) andHwang et al. (2022a).The dash-dotted line represents the analytically derived  () for wide binaries byXu et al. (2023).

Figure 6 .
Figure 6.The power-law index  of the  distribution at SFE = 5% as a function of the separation.The circular markers within the rectangular boxes represent the 50 th percentile estimate of  in our simulations in the separation range denoted by the width of the boxes.The height of the boxes represents the 16 th to 84 th percentile range of .The triangular markers correspond to the  values measured in Hwang et al. (2022a) at different binary separations, and the solid curve represents the fit to their data points.

Figure 7 .N−N−N−Figure 8 .
Figure 7.The  distribution at SFE = 5% for binaries with system masses  sys ≤ 0.8 M ⊙ , i.e., binaries with mainly M-dwarfs and later type stars (left) and  sys > 0.8 M ⊙ , i.e., binaries with mainly solar and earlier type stars (right).The low-mass systems ( sys ≤ 0.8 M ⊙ ) have highly superthermal  distribution while the relatively higher-mass systems have subthermal/thermal distribution.

Figure 9 .
Figure 9. System mass ( sys vs. semi-major axis () at SFE = 5%.The dashed line represents  sys = 0.8 M ⊙ .The marker colour is scaled to the time of formation of the binary ( form ) in units of the simulation end time ( end ).

Table 1 .
Summary of eccentricity distribution results for different cases.