Close binary fractions in accreted and in situ halo stars Open Access

ABSTRACT

The study of binary stars in the Galactic halo provides crucial insights into the dynamical history and formation processes of the Milky Way. In this work, we aim to investigate the binary fraction in a sample of accreted and in situ halo stars, focusing on short-period binaries. Utilizing data from Gaia Data Release 3 (DR3), we analysed the radial velocity uncertainty |$\sigma _{\mathrm{RV}}$| distribution of a sample of main-sequence stars. We used a novel Bayesian framework to model the dependence in |$\sigma _{\mathrm{RV}}$| of single and binary systems allowing us to estimate binary fractions F in a sample of bright (⁠|$G_{\mathrm{RVS}}$| < 12) Gaia sources. We selected the samples of in situ and accreted halo stars based on estimating the 6D phase space information and affiliating the stars to the different samples on an action angle versus energy (⁠|$L_{\mathrm{z}}{\!-\!}E$|⁠) diagram. Our results indicate a higher, though not significant, binary fraction in accreted stars compared to the in situ halo sample. We further explore binary fractions using cuts in E and |$L_z$|⁠, and find a higher binary fraction in both high-energy and prograde orbits that might be explained by differences in metallicity. By cross-matching our Gaia sample with APOGEE DR17 catalogue, we confirm the results of previous studies on higher binary fractions in metal-poor stars and find the fractions of accreted and in situ halo stars consistent with this trend. Our finding provides new insights into binary stars’ formation processes and dynamical evolution in the primordial Milky Way Galaxy and its accreted dwarf Galaxies.

1 INTRODUCTION

The Galactic halo, a vast and sparsely populated region enveloping the Milky Way, presents a unique environment for the study of stellar populations (Helmi 2008, 2020; Deason & Belokurov 2024). It comprises predominantly old stars, remnants of both the primordial galactic formation era and the subsequent accretion events. This dichotomy in the halo’s star formation history is epitomized in the two distinct groups of stars: ‘accreted’ halo stars, originating from smaller, consumed galaxies and the ‘in situ’ halo stars, formed within the nascent Milky Way (Deason & Belokurov2024).

The accreted component of the halo in the solar neighbourhood is exemplified by the Gaia-Sausage-Enceladus (GS/E; Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018), a prominent substructure discerned through analyses of Gaia data (Gaia Collaboration 2016, 2018) indicating the remnant of the Milky Way’s most consequential recent merger, an event estimated to have concluded around |$8{\!-\!}10$| billion years ago. The GS/E appears to dominate the halo around the Sun, but a wide variety of smaller halo substructures have also been detected (Helmi et al. 1999; Myeong et al. 2018a, b; Koppelman et al. 2019; Matsuno, Aoki & Suda 2019; Myeong et al. 2019; Naidu et al. 2020; O’Hare et al. 2020; Malhan 2022; Dodd et al. 2023).

Conversely, the in situ halo population includes the ‘splashed’ high-|$\alpha$| disc (Bonaca et al. 2017; Di Matteo et al. 2019; Gallart et al. 2019; Belokurov et al. 2020), and the Aurora stars (Belokurov & Kravtsov 2022; Conroy et al. 2022; Chandra et al. 2024; Zhang, Ardern-Arentsen & Belokurov 2024), both formed in the Milky Way proper before the arrival of the GS/E progenitor. Splash stars, with their relative metal-rich composition ([Fe/H] |$\gtrsim -1$|⁠) and highly eccentric orbits, are thought to have been displaced from the Galactic disc during the cataclysmic merger event that formed the GS/E structure. Having been formed before the merger, Splash stars are on average as old as the stars in the GS/E (see e.g. Belokurov et al. 2020). Aurora is characterized by its metal-poor nature ([Fe/H] |$\le$| −1) and a wide range of azimuthal velocities (Belokurov & Kravtsov 2022; Chandra et al. 2024; Zhang et al. 2024). Identified through chemical tagging, Aurora represents the oldest in situ stellar population of the Milky Way, formed before the Galaxy had developed a coherently rotating disc (see discussion in Dillamore et al. 2024; Semenov et al. 2024). Although first spotted in the vicinity of the Sun, Aurora likely follows a steep radial density, with most of its mass near the Galactic Centre (see Rix et al. 2022; Belokurov & Kravtsov 2023).

Alongside the accreted stars, these components of the in situ halo, Splash and Aurora, underscore the Milky Way’s dynamic history, from its earliest star formation episodes to the significant reshaping triggered by galactic mergers. This layered approach facilitated by the Gaia space observatory data provides a nuanced understanding of how different stellar populations within the Galactic halo offer insight into the chronological sequence of events that shaped the Milky Way.

In the context of the Galactic halo, the study of binary stars acquires additional significance. Binary stars are a fundamental component of our Galaxy’s architecture. They serve as critical probes for understanding stellar formation and evolution (Carney & Latham 1987; Latham et al. 2002; Tokovinin et al. 2006; Raghavan et al. 2010; Moe, Kratter & Badenes 2019; Mazzola et al. 2020; Price-Whelan et al. 2020; Bashi, Mazeh & Faigler 2023; Gaia Collaboration 2023b; Chen, Liu & Han 2024). The frequency and characteristics of binary systems in this region provide insights into the conditions of star formation and the dynamical processes in the early Galaxy and the smaller, assimilated dwarf galaxies. For a recent review of results in the study of binary star research, particularly those facilitated by Gaia, we refer the reader to El-Badry (2024).

Previous studies have predominantly focused on binary systems within the disc and bulge of the Milky Way or have treated the halo as a homogeneous entity. One exception might be the recent work of Hwang et al. (2022), which focused on wide-binary stars fractions.

In recent years, the field of astrophysics has seen an influx of data from wide-field spectroscopic surveys such as LAMOST (The Large Sky Area Multi-Object Fiber Spectroscopic Telescope; Cui et al. 2012), GALAH (Galactic Archaeology with HERMES; De Silva et al. 2015), and APOGEE (Apache Point Observatory Galactic Evolution Experiment; Majewski et al. 2017). These surveys have been instrumental in measuring the radial velocities (RVs) and chemical abundances of hundreds of thousands of stars. A key aspect of some of these surveys is the collection of multi-epoch data for individual stars. This approach has enabled a statistical analysis of RV variability as a function of various stellar properties, including metallicity (Tian et al. 2018; Moe et al. 2019; Price-Whelan et al. 2020) and |$\alpha$| elements (Badenes et al. 2018; Mazzola et al. 2020; Niu et al. 2021).

For this study, we meticulously selected our sample from Data Release 3 (DR3; Gaia Collaboration 2023a), focussing on specific criteria to ensure the relevance and accuracy of our analysis. While Gaia DR3 does not include each source’s individual epoch RV measurements (Katz et al. 2023), we can exploit the available RV variability information such as the standard deviation of the epoch RV measurements of the individual source |$\sigma _{\mathrm{RV}}$|⁠, or the |$\Delta \mathrm{RV_{max}}$|⁠, the maximum detected shift in the RVs to assess binary frequency in a statistical framework (Maoz, Badenes & Bickerton 2012; Badenes et al. 2018; Moe et al. 2019; Mazzola et al. 2020; Andrew et al. 2022).

While it is generally hypothesized that accreted and in situ halo stars might exhibit different binary fractions due to their distinct formation histories, empirical evidence supporting this notion is sparse. This research aims to fill this gap by employing a novel methodological approach to estimate the close binary fractions among these stars, thereby shedding light on the diverse evolutionary pathways that have shaped the Galactic halo.

The paper is structured as follows. In Section 2, we present our sample selection of bright Gaia main-sequence (MS) sources. Section 3 discusses the dependence between |$\sigma _{\mathrm{RV}}$| and the stellar magnitude seen in the Gaia sample. In Section 4, we demonstrate our selection procedure of accreted and in situ halo stars. Section 5 presents our methodological approach using Bayesian modelling of the binary fractions. Sections 6 and 7 give our main results on the binary fractions in the accreted and in situ, as well as their dependence on metallicity. Section 8 summarizes our results and discusses their possible implications.

2 SELECTING A SAMPLE OF BRIGHT MAIN-SEQUENCE SOURCES

We targeted stars within 3 kpc of the solar vicinity with measurable RVs (Katz et al. 2023). The query used to obtain data from the Gaia DR3 was designed to filter and isolate a specific subset of MS and sub-Giant objects, focussing on those relevant to our study of binary fractions in the Galactic halo. Consequently, we limited our sample to stars with magnitudes estimated from the integration of the radial velocity spectrometer (RVS) |$G_{\mathrm{RVS}}$| brighter than 12 mag using rv_method_used = 1. |$G_{\mathrm{RVS}}$| is estimated using the flux in the RVS spectra (Sartoretti et al. 2023). This magnitude constraint was applied to ensure the reliability and quality of the RV data, as the epoch RVs of stars beyond this threshold were not considered reliable enough to derive the combined RVs as the median of the epoch RVs (Katz et al. 2023).

We further selected stars with well-constrained parallax measurements (parallax_over_error > 5 and parallax > 0), ensuring high-quality distance estimates. To minimize contamination from disc stars, we selected stars with a galactic latitude (b) satisfying the condition |$|b|\gt 10 ^{\circ }$|⁠. This criterion was essential to focus our study on halo stars and reduce the inclusion of stars from the galactic disc with different dynamical properties.

We used Gaia extinction values to accurately plot the colour–magnitude diagram (CMD) of our sample. The inclusion of stars with non-null extinction values ensures that we have reliable photometric extinction-corrected data, crucial for accurate luminosity and temperature determination. Specifically, we corrected the absolute G-band magnitudes (Abs. G) and BP–RP colours using the ag_gspphot and ebpminrp_gspphot values, respectively, to accurately represent the absolute magnitudes on the Y-axis and the BP–RP colours on the X-axis of the CMD plot.

Within the CMD, we employed a polygonal selection method to isolate mainly dwarf stars. Our decision to focus on MS stars stemmed from our objective of studying stellar variability caused by motion in a binary system. Using this selection, we avoided the added complexity in cases of evolved stars with their variable and evolutionary transient stages as well as the impact of stellar evolution on binary fraction.

To further constrain our sample to select stars that were observed with good coverage in epoch phase space, we selected sources with rv_visibility_periods_used larger than eight, where a visibility period is a group of transits separated from other groups by a gap of at least 4 d (Katz et al. 2023; Gaia Collaboration 2023a).

In the left panel of Fig. 1, we show a CMD of our sample of |$N=1,153,051$| MS and subgiants sources (out of a total stellar sample of |$N=2174\,836$| sources) delineated by green polygonal borders.¹ The colour scale represents the number of sources falling into each bin. The right panel of Fig. 1 displays a zoomed selection of the polygonal area, where we have superimposed the mesa Isochrones and Stellar Tracks (MIST) models (Paxton et al. 2011; Choi et al. 2016; Dotter 2016), which include isochrones interpolation for single stars (i.e. binary mass ratio, |$q=0$|⁠) and twin binaries (⁠|$q=1$|⁠). For single stars, the tracks span a range of metallicities, accurately tracing the progression of metal-poor to metal-rich stars within the MS. For the [Fe/H] = −1 and [Fe/H] = −2, a 12 Gyr isochrone was selected, while for the [Fe/H] = 0 we selected a 4.5 Gyr isochrone to better characterize the current population of stars in these metallicity bins. Additionally, for the case of twin binaries, the corresponding track is plotted 0.75 mag above the single star line (Hurley & Tout 1998) representing the photometric appearance of an equally bright companion in the CMD, effectively mimicking the luminosity of a normal metal-rich dwarf.

As can be seen by Fig. 1, the polygon borders are meticulously chosen to include mainly dwarf stars and capture the sequence of metal-poor stars, typically found at the lower extremity of the MS. Notably, the binary sequence of metal-rich stars (at BP_RP > 1) is clearly discernible within our selection and corresponds nicely to the twin stars [Fe/H]  = 0 isochrone. While most stars in this region are not halo stars, their inclusion was deemed essential for completeness. As a final remark, we note that in the case of a metal-poor star with an equally bright companion (twin), the photometric position on the CMD appears as a normal metal-rich dwarf.

3 STELLAR VARIABILITY: RADIAL VELOCITY UNCERTAINTY

To get the standard deviation of the epoch RV measurements of the individual bright source |$\sigma _{\mathrm{RV}}$|⁠, we used the following expression:

where |$\zeta _{\mathrm{RV}}$| is the uncertainty on the median of the epoch RVs (radial_velocity_error) to which a constant shift of 0.11 km s⁻¹ was added to take into account a contribution from the calibration floor (Andrew et al. 2022; Katz et al. 2023), and |$N_{\mathrm{RV}}$| is the number of transits used to derive the median RV (rv_nb_transits).

Fig. 2 shows the distribution of |$N=1153\,051$| MS sources in |$\sigma _{\mathrm{RV}}$| as function of |$G_{\mathrm{RVS}}$|⁠. As might be expected, the typical standard deviation of the epoch RV measurements increases towards fainter sources. Above this observable curve, we find a wider spread of sources of suspected binary systems with high |$\sigma _{\mathrm{RV}}$| values.

To demonstrate this expected trend, we show in Fig. 3 the |$\sigma _{\mathrm{RV}}$|–|$G_{\mathrm{RVS}}$| scatter of a sample of three types of systems: (i) fraction of planet-host stars listed in the NASA Exoplanet Archive² with a flag suggesting detection by RV and can be approximated as a well-constrained sample of single stars; (ii) a sample of single-lined spectroscopic binaries (SB1s) listed in the Gaia DR3 NSS catalogue (Gaia Collaboration 2023b) with a clean-score above 0.578 based on Bashi et al. (2022) metric; and (iii) eclipsing binaries (EBs) listed in the first Gaia catalogue of EB candidates (Mowlavi et al. 2023). As expected, the single-star sample sits nicely on the exponential curve and have the lowest |$\sigma _{\mathrm{RV}}$| for a given |$G_{\mathrm{RVS}}$|⁠, while the SB1s are found above this curve, and the EBs with typically shorter periods and larger semi-amplitudes have the highest |$\sigma _{\mathrm{RV}}$| values.

In Section 5, we will further elaborate on this dependence in developing our model to estimate binary fractions and the typical (mean) |$\sigma _{\mathrm{RV}}$| value as a function of |$G_{\mathrm{RVS}}$|⁠.

4 ACCRETED AND IN SITU SAMPLES

In addition to the criteria mentioned above, we implemented further cuts to refine our selection of halo stars. Specifically, we applied a velocity threshold, selecting stars with a total velocity (⁠|$V_{\mathrm{tot}} = \sqrt{U^2 + V^2 + W^2}$|⁠) greater than 220 km s⁻¹. This threshold is based on the understanding that the halo shows close to net zero rotation, and therefore in the Solar neighbourhood, the velocity of its stars is on average offset from that of the disc’s stars by the local value of the circular velocity (Nissen & Schuster 2010; Belokurov et al. 2018; Helmi et al. 2018; Bashi & Zucker 2022; Lane, Bovy & Mackereth 2022).

To accurately estimate the dynamical properties of our selected stars, we used the galpy package (Bovy 2015). Using galpy, we calculated the orbital energy (E) and the z-component of angular momentum (⁠|$L_z$|⁠) for each star in our sample. These calculations were based on the assumption of the Milky Way’s potential as described in the MWPotential2014 model, which is a well-established model that approximates the gravitational potential of the Milky Way, considering the contributions from the bulge, disc, and halo. We plot all energies offset by the potential of MWPotential2014 at infinity so that |$E\lt 0$| are bound orbits (e.g. Lane et al.2022).

In the |$E{\!-\!}L_z$| plane, the samples are confined by a bounding parabola defined by the escape velocity as a function of the change in energy at a fixed Galactocentric radius.

In Fig. 4, we show |$N=1929$| stars on the |$E{\!-\!}L_z$| plane that were left after applying the above cuts and a background grey contours marking the position of the overall MS sample that mainly consisted of disc stars. We mark the Sun’s position at |${\rm E}_{\odot }=-1.04 \times 10^5~\mathrm{km^2\,s^{-2}}$| and use it to gauge the boundary energy we select to separate between the accreted and in situ samples. Follow Belokurov & Kravtsov (2023), and to reduce contamination, we selected conservative cuts in the E and affiliated accreted stars as sources with |$E\, \gt\, {\rm E}_{\odot } + 0.2$| and |$L_z\, \lt \, 2$| and in situ stars as |$E \,\lt \, {\rm E}_{\odot }-0.2$| and |$L_z \lt 0$|⁠.

By incorporating these velocity and dynamic cuts, we were able to identify stars that are characteristic of the Milky Way’s halo more precisely and were left with |$N=181$|accreted and |$N=291$|in situ halo stars. This refined selection ensures that our analysis of binary fractions is focused on the intended stellar population, thereby enhancing the reliability of our binary fractions estimates.

5 METHODOLOGY: MODELLING BINARY FRACTIONS USING THE RV UNCERTAINTY DISTRIBUTION

Single stars typically exhibit minimal RV variation, primarily due to factors such as stellar activity and instrumental noise, resulting in low |$\sigma _{\mathrm{RV}}$|values that follow a trend similar to the cyan points of Fig. 3 with fainter stars having a larger variation. On the other hand, in binary systems, the orbital motion of each star causes periodic changes in their RVs. For a given binary semi-amplitude, K, in the most simplified scenario of a circular orbit, assuming that the measurements are uniformly distributed across the orbital phase, the added RV uncertainty from the binary nature of the system would be |$\sigma _{\mathrm{RV}}^b=K/\sqrt{2}$|⁠. Thus, the standard deviation of the value of the epoch RV measurements is a direct indicator of the binary nature of the system (Maoz et al. 2012; Badenes et al. 2018; Moe et al. 2019; Mazzola et al. 2020; Andrew et al. 2022), where a larger |$\sigma _{\mathrm{RV}}$| typically suggests a closer or more massive companion as |$K \propto P^{-1/3} q (1+q)^{-2/3}$|⁠, where P is the binary orbital period and q is the binary mass ratio.

By analysing the |$\sigma _{\mathrm{RV}}$|–|$G_{\mathrm{RVS}}$| distribution for a given sample of Gaia stars, we can statistically estimate the binary fraction. Given a typical |$\sigma _{\mathrm{RV}}$| of |$\sim 5$| km s⁻¹ at |$G_{\mathrm{RVS}}$| = 12, our method is sensitive to close binary systems with orbital period |$P \lesssim 1000$| d for a solar mass twin binary system. None the less, it important to note that detecting twin sources with mass ratios close to one (⁠|$q \simeq 1$|⁠) using Gaia RVS is more challenging because the Doppler shifts of the two similar stars partially cancel each other out in the combined spectrum. This cancellation reduces the observable RV variability, and consequently the |$\sigma _{\mathrm{RV}}$|value.

A data set consisting predominantly of single stars would exhibit a narrow distribution centred around and following the trend of dense sources shown in Fig. 2, thereby reflecting minimal RV variation. In contrast, a data set with a high binary fraction would show a broader distribution, extending towards higher |$\sigma _{\mathrm{RV}}$|values for a given |$G_{\mathrm{RVS}}$|⁠, indicative of the significant RV changes caused by binary companions. A mixed population of single and binary stars would likely display a bimodal distribution, with two distinct peaks corresponding to the two populations.

In most cases, our method treats higher order multiple-star systems, such as triples, as binary stars. In this scenario, we might be sensitive to the close binary substructure where the ratio of the outer period |$P_{\mathrm{out}}$| to the inner period |$P_{\mathrm{in}}$|⁠, adjusted for eccentricity |$e_{\mathrm{out}}$|⁠, is expected to have |$P_{\mathrm{out}}\left(1-e_{\mathrm{out}} \right)^3/P_{\mathrm{in}} \gt 5$| (Tokovinin 2004) while we may miss the RV variation of the outer companion.

To model the population of sources in log |$\sigma _{\mathrm{RV}}$|⁠, we define a density function, which is a sum of two Gaussian distributions to model the populations of single and binary stars. To properly model the fractions of single and binary stars, we multiply each Gaussian by the expected fractions, where F marks the binary fraction. The final density function will then be

where |$\theta = {\left(F, a, b, G_{\mathrm{min}}, d, \sigma _s, \sigma _b\right)}$| defines our parameter space and |$\mathcal {N}_s\left(x| \mu _s(G_{\mathrm{RVS}}), \sigma _s^2\right)$|⁠, |$\mathcal {N}_b\left(x| \mu _b(G_{\mathrm{RVS}}), \sigma _b^2\right)$| are the Gaussian density functions of the single and binary populations with means |$\mu _s(G_{\mathrm{RVS}})$|⁠, |$\mu _b(G_{\mathrm{RVS}})$| and standard deviations |$\sigma _s$|⁠, |$\sigma _b$|⁠, respectively.

To model the dependence between |$\sigma _{\mathrm{RV}}$| and |$G_{\mathrm{RVS}}$|⁠, we used an exponent function and defined |$\mu _s(G_{\mathrm{RVS}})$| as

where |$a,~ b,~ G$| are free parameters which control the shape of the exponential function. Similarly, we defined the binary mean, |$\mu _b(G_{\mathrm{RVS}})$| in quadrature as

where d is a free parameter that accounts for the extra variability in |$\sigma _{\mathrm{RV}}$| as a result of the system’s binary nature.

In our analysis, we sought to account for potential variations in |$\sigma _{\mathrm{RV}}$| attributed to line-broadening where early-type stars at a fixed |$G_{\mathrm{RVS}}$| have larger errors (Andrew et al. 2022). To explore this dependence, we extend our model by introducing an additional factor to account for the relation between |$\sigma _{\mathrm{RV}}$| and BP_RP. Despite this, the outcomes from our ‘simplified’ model were found to align with those from a model that explicitly incorporates the spectral type. This agreement suggests that the parameters |$\sigma _s$| and |$\sigma _b$| within our model effectively encapsulate the influence of spectral type, thereby moderating the dispersion of data points around the fitted single- and binary-star models.

Finally, since stars whose combined RV formal uncertainty was greater than 40 km s⁻¹ were flagged as invalid and removed from the final Gaia DR3 catalogue (Katz et al. 2023), we used a truncated normal distribution to cope with this abrupt cut in the observed distributions.

The final likelihood expression would then be

To test our algorithm, we applied it on different sets of simulated data sets and used a Markov chain Monte Carlo routine with 40 walkers and |$10^4$| steps, using uninformative priors on the seven-free parameters of our model as listed in Table 1. We used the pythonemcee package (Foreman-Mackey et al. 2013) to find the parameter values (and their uncertainties) that maximize the sample likelihood.

In the first test set, we attempted to find the best-fitting model parameters |${\left(F, a, b, G_{\mathrm{min}}, d, \sigma _s, \sigma _b\right)}$| from a simulated sample based on the defined density function of equation (2).

In the second test set, we simulated a population of single versus short-period binary systems using a random selection of binary separations, mass ratios, eccentricity, distribution, inclination, and phase (e.g. Mazzola et al. 2020) and attempted to get the best-fitting parameters of |${\left(F, a, b, G_{\mathrm{min}}\right)}$|⁠.

In the third test set, we used a mixed sample of known single stars and Gaia SB1s (Bashi et al. 2022; Gaia Collaboration 2023b) based on the subsamples used in Fig. 3 and attempted to get the best-fitting parameters of |${\left(F\right)}$|⁠.

We found good agreement between the posterior distributions and the actual simulated parameters in all test cases, thereby supporting our Bayesian model and encouraging us to proceed and estimate binary fractions on real data sets.

6 BINARY FRACTIONS IN HALO STARS

6.1 The accreted and in situ subsamples

Given our selection of |$N=181$| and |$N=291$|accreted and in situ sources, respectively, we list in Table 2 the best-fitting posterior values of our Bayesian model (median and 16 per cent, and84 per cent quantiles uncertainties). We present in Fig. 5 the posterior distributions of the accreted (orange) and in situ (green) binary fractions. We find a |$\sim 2\sigma$| difference in favour of a higher binary fraction in accreted stars, |$F_{\mathrm{accreted}}=0.238\pm 0.038$|⁠, compared to in situ stars, |$F_{\mathrm{in-situ}}=0.182\pm 0.032$|⁠.

In Fig. 6, we show a scatter plot of |$\sigma _{\mathrm{RV}}$| versus |$G_{\mathrm{RVS}}$| for the accreted and in situ samples. The red dotted lines mark our best-fitting posterior distributions of the mean values using equation (3) that describes the relation between |$\sigma _{\mathrm{RV}}$| and |$G_{\mathrm{RVS}}$| in the case of a single source solution. In general, the best-fitting posteriors (a, b, |$G_{\mathrm{min}}$|⁠), which determine the shape of the curves of the accreted and in situ samples, are consistent within |$1\sigma$| as can be seen in Table 2.

To explore the distributions of |$\sigma _{\mathrm{RV}}$| after removing the effect of stellar |$G_{\mathrm{RVS}}$|⁠, we show in Fig. 7 the ‘flattened’ normalized distributions after subtracting the observed |$\sigma _{\mathrm{RV}}$| values with the best-fitting model |$\mu _s(G_{\mathrm{RVS}}; a, b, G_{\mathrm{min}})$|⁠. Comparing these distributions might also shed some light on the nature of these binary systems. As expected, both histograms have clear peaks around zero, marking the distributions of single sources with standard deviation |$\sigma _s$|⁠. The extended tail of binaries is also evident where we find the accreted sample having more sources at high |$\sigma _{\mathrm{RV}}-\mu _s(G_{\mathrm{RVS}}; a, b, G_{\mathrm{min}})$| values compared to the in situ subsample. This observed trend is also supported by our posterior model results where the accreted best-fitting model on d and |$\sigma _b$|⁠, which are responsible to mark the position and spread of the binary sequence, is larger than in the in situ case.

A significant limitation encountered in using the Gaia sample relates to the obscured insights into the metallicity content of the two halo samples. Metallicity is known to anticorrelate with binary fraction (Badenes et al. 2018; Moe et al. 2019; Mazzola et al. 2020; Price-Whelan et al. 2020). Consequently, any disparity in the metallicity content between the accreted and in situ samples can lead to variations in binary fractions that are not inherently related to their dynamical properties. To further explore this, we examine in Section 7 a subset of our samples with reported metallicity estimates using the APOGEE DR17 catalogue (Abdurro’uf et al. 2022).

6.2 Cuts in E and |$L_z$|

We next attempted to estimate binary fractions using different cuts in E and |$L_z$| that might be associated with different Galactic components (Helmi & White 1999; Belokurov et al. 2018; Matsuno et al. 2019; Myeong et al. 2019; Naidu et al. 2020; Belokurov & Kravtsov 2022, 2023). While accreted stars and especially the GS/E are more common at higher E values in the solar neighbourhood (Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Naidu et al. 2020), other substructures also have their expected places on the |$E{\!-\!}L_z$| diagram. For example, at high |$L_z$|⁠, we are evident to the extended high-|$\alpha$| sequence of disc stars as well as Splash stars (Bonaca et al. 2017; Di Matteo et al. 2019; Gallart et al. 2019; Belokurov et al. 2020). This sample might be contaminated by a higher fraction of metal-rich stars. Similarly, at low-energy values of retrograde orbits, we should expect to have a higher population of Aurora, in situ stars (Belokurov & Kravtsov 2022; Conroy et al. 2022; Chandra et al. 2024; Zhang et al. 2024).

In the upper panel of Fig. 8, we show the binary fractions using different cuts in energy, |$E^{\mathrm{cut}}$| that separate our sample to high- and low-energy sources (see red and violet points, respectively). Our analysis demonstrates that below a critical energy threshold (⁠|$E \lt E^{\mathrm{cut}} \approx {\rm E}_{\odot }$|⁠), binary fractions remain relatively similar. This consistency can be explained by the high contamination rate of in situ stars which dominates the sample of high-energy stars, making the two subsamples of low- and high-energy similar in stellar properties. Conversely, above this energy limit (⁠|$E \gt E^{\mathrm{cut}} \approx {\rm E}_{\odot }$|⁠), we observe a noticeable increase in binary fractions. This trend intensifies at higher energy levels, where the sample likely comprises a purer selection of accreted stars. These results align with the findings of our previous section, where we estimated binary fractions in a conservative subsample of accreted and in situ stars. The increase in binary fraction at higher energy bins suggests a more significant presence of dynamically heated, accreted structures especially of the GS/E, with lower metallicities.

The variation of binary fractions with angular momentum also presents an interesting trend. In the bottom panel of Fig. 8, we show the binary fractions as function of the angular momentum cut-off (⁠|$L_z^{\mathrm{cut}}$|⁠) that separates our sample to retrograde-prograde and slow-fast rotation sources (see red and violet points, respectively). We find consistent higher binary fractions in sources with |$L_z \lt L_z^{\mathrm{cut}}$| compared to the complement subgroup. In stark contrast, the most prograde sources, predominantly associated with the high-|$\alpha$| disc and Splash (Bonaca et al. 2017; Di Matteo et al. 2019; Gallart et al. 2019; Belokurov et al. 2020; Dodd et al. 2023), exhibit significantly lower binary fractions than the overall halo sample. Possible explanations for this difference might include the higher metal content associated with this subsample compared to the slow and retrograde sources. Furthermore, as recently suggested by Mazzola et al. (2020), high abundances of |$\alpha$| elements also suppress multiplicity at most values of [Fe/H]. Consequently, the combination of higher metallicity and elevated |$\alpha$| element abundances likely contributes to the lower binary fractions observed.

7 METALLICITY DEPENDENCE

In this section, we delve into the relationship between metallicity and binary star fractions while considering samples of accreted and in situ halo stars. We used the APOGEE DR 17 catalogue (Abdurro’uf et al. 2022) to obtain estimates of stellar metallicity and cross-matched our Gaia MS sample with it. To minimize effects by stellar mass on binary fraction and to select a representative sample similar to our halo sample, we limited our selection to sources with APOGEE effective temperature (TEFF) between |$4600\,{\text{and}}\,6200$| K and stellar mass, approximated using the Torres, Andersen & Giménez (2010) relation, within the |$0.5{\!-\!}1.2{\rm M}_{\odot }$| range.

On the one hand, while there are larger catalogues with reported metallicity estimates (e.g. Andrae, Rix & Chandra 2023), it is extremely important when dealing with metal-poor binary stars to use a high spectroscopic derivation of the metal content to avoid overestimation of the derived metallicities. On the other hand, performing our binary search analysis solely using APOGEE stars would not be applicable since most targets had only a single exposure with a small fraction of sources having 2–4 RV epoch observations (Mazzola et al. 2020; Abdurro’uf et al. 2022). This limited number of epochs is insufficient to get good estimates of |$\sigma _{\mathrm{RV}}$| compared to Gaia’s DR3 typical |$20{\!-\!}30$| epochs (Katz et al. 2023).

As part of our initial filtering, we removed stars with measurements that are believed to be problematic, as indicated by the following general flags: VERY_BRIGHT_NEIGHBOR, LOW_SNR, as well as these element pipeline analysis flags: FE_H_FLAG, MG_FE_FLAG, AL_FE_FLAG. However, we found that even while using these cuts, there were a few cases where the extracted metallicity estimates listed in APOGEE were underestimated. In most of these cases, the outliers strongly correlated with high values of |$\sigma _{\mathrm{RV}}$|⁠, suggesting that the binary nature of these sources, especially in cases of high mass ratios, have biased the elemental estimation. These cases are common particularly in blended sources and multiple-star systems with more than one star that contribute significantly to the overall flux. A similar issue was also highlighted in Price-Whelan et al. (2020), where, in their analysis aimed at characterizing close binaries, it was assumed that their sample consisted exclusively of single-line sources.

Consequently, to remove these outliers, we found the value on RV_CCFWHM, which marks the full width at half-maximum (FWHM) of the cross-correlation peak from the combined versus best-match synthetic spectrum, as a good indicator of ‘well-behaved’ solutions in cases with values below 150 km s⁻¹. The cross-correlation technique is employed to match observed stellar spectra against a library of synthetic spectra. The FWHM of the cross-correlation peak, denoted by RV_CCFWHM, quantifies the width of this peak at half of its maximum height, reflecting the degree of match between the observed and synthetic spectra. Thus, a peak that is too broad might indicate an imprecise match.

The cross-match of the Gaia-APOGEE sample, yielded |$N=48\,484$| sources (disc and halo) in common, with only |$N=9$| and |$N=27$|accreted and in situ sources, respectively, having using Section 4 subsample selection. In the left panel of Fig. 9, we show a scatter plot of the sources on the |$E{\!-\!}L_z$| plane similar to Fig. 4 while on the right panel we show the distribution on the [Fe/H]–[Al/Fe] diagram. As expected, our kinematic selected halo subsamples dominate the metal-poor fraction of stars, with accreted sources having lower [Al/Fe] values compared to in situ stars (Hasselquist et al. 2021; Belokurov & Kravtsov 2022, 2023; Deason & Belokurov 2024). In addition, we note on a few other in situ metal-poor prograde stars seen in this plot (grey points), probably affiliated with the Aurora population, which while assumed to be kinematically hot, has an approximately isotropic velocity ellipsoid and a modest net rotation (Belokurov & Kravtsov 2022; Zhang et al. 2024).

Using the Gaia-APOGEE sample, we can attempt to estimate binary fractions in bins of metallicity. The black crosses in Fig. 10 present a detailed analysis of the binary fraction in various metallicity bins. We focus on the median values and the 16 per cent and 84 per cent quantile values of the binary fractions posteriors. In addition, we plot by red and violet crosses the fraction of retrograde (⁠|$L_z \lt 0$|⁠; |$N=51$|⁠) and prograde (⁠|$L_z \ge 0$|⁠; |$N=83$|⁠) halo sources, respectively.

The low number of accreted (⁠|$N=9$|⁠) and in situ (⁠|$N=27$|⁠) halo stars in our Gaia-APOGEE cross-matched subsample is not sufficient to produce valuable estimates of binary fractions using our Bayesian model. Instead, we decided to estimate the median and 16 per cent and 84 per cent quantile values on the metallicity values in each halo subsample and place the binary fractions found in Section 6 as the expected values. Consequently, we mark by orange and green boxes the accreted and in situ binary fractions to stress that these are not a direct calculation of the binary fractions of corresponding samples shown in Fig. 9.

Examining the trend in field stars’ binary fraction as a function of [Fe/H], Fig. 10 corroborates previous research findings that have consistently shown higher binary fractions in metal-poor stars (Badenes et al. 2018; Moe et al. 2019; Mazzola et al. 2020; Price-Whelan et al. 2020). This trend is prominently visible across the different metallicity bins, reinforcing the significant role of metallicity in the formation of binary stars.

In our context, one further notable observation from Fig. 10 is the consistency of the accreted and in situ with the field stars trend, suggesting that the anticorrelation in close-binary fractions with metallicity continues into the domain of halo stars. Secondly, we find an interesting discrepancy in field binary fractions compared to the halo prograde stars. A similar discrepancy is also evident when comparing the trend of field binary fractions in Fig. 10 to the highest prograde bin (⁠|$L_z \ge L_z^{\mathrm{cut}}=0.5$|⁠) in the bottom panel of Fig. 8. Assuming that the prograde bin includes high-|$\alpha$| disc and Splash stars, we might expect it to be populated mainly by sources with [|$\alpha$|/Fe] |$\ge 0.2$|⁠. However, our bin estimate of a binary fraction of |$F\approx 0.10$|⁠, is lower than the expected value observed in the trend of Fig. 10 for similar metallicity. None the less, given the Gaia-APOGEE sample is dominated by thin disc stars, most of the sources in this diagram have lower |$\alpha$| values and have a higher stellar mass compared to thick disc stars. Both of these parameters inspire a higher binary fraction (Raghavan et al. 2010; Moe et al. 2019; Mazzola et al. 2020) that can mitigate this effect.

8 CONCLUSIONS

This study provides valuable insight into close binary star fractions within the Galactic accreted and in situ halo populations, which stands as a pivotal component in understanding the dynamical history and formation of the Milky Way Galaxy.

By leveraging the comprehensive data from Gaia DR3 (Gaia Collaboration 2023a), our methodological approach, utilized a novel Bayesian framework to analyse the RV uncertainty distribution of Gaia RVS sources (Katz et al. 2023; Sartoretti et al. 2023). Using this approach, we were able to model the populations of single and binary stars as a sum of two Gaussian distributions, allowing a nuanced understanding of close-binary star fractions in different Galactic environments using only Gaia data.

Our findings reveal a higher binary fraction in accreted stars, compared to in situ halo stars. On the one hand, this observation suggests that the processes that influence binary star formation and evolution might differ across these distinct populations. On the other hand, given the slightly lower metal content in accreted stars compared to in situ halo stars, the difference in binary fractions might be explained by the difference in metallicity. A similar conclusion was also reached for wide-binary fractions in Hwang et al. (2022), although their selection of accreted and in situ sources was different from the one we used. Furthermore, their selection of stars on the |$E{\!-\!}L_z$| plane was more contaminated, causing the binary fractions to be much more consistent in any case.

In our study, we also explored binary fractions across different galactic components, employing varied energy and angular momentum cuts to differentiate these groups. We found that, below a specific energy threshold (⁠|$E \lt E^{\mathrm{cut}} \approx E_{\odot }$|⁠), binary fractions remained relatively uniform due to a high contamination rate of in situ stars that dominate the high-energy sample (Naidu et al. 2020; Belokurov & Kravtsov 2023; Deason & Belokurov 2024; Zhang et al. 2024), resulting in similar stellar properties between low- and high-energy stars. Above this energy threshold, binary fractions increased, especially at higher energy levels where the sample likely consists more of accreted stars, of the GS/E, indicating a prevalence of dynamically heated metal-poor structures. Additionally, when examining the variation of binary fractions with angular momentum, we observed consistently higher binary fractions in retrograde sources (⁠|$L_z \lt L_z^{\mathrm{cut}}$|⁠) compared to their prograde counterparts. This disparity is stark, as the most prograde sources, which are typically metal-rich, exhibited significantly lower binary fractions than those in the overall halo sample. These findings suggest that the observed trends in binary fractions can provide insights into the dynamic and compositional differences within the Galaxy.

In general, our research corroborates previous studies, indicating higher binary fractions in metal-poor stars (Moe et al. 2019; Mazzola et al. 2020; Price-Whelan et al. 2020). Our analysis reveals that the binary fractions of accreted and in situ samples are consistent with the field stars trend, as evidenced in Fig. 10. This observation suggests that the anticorrelation between close-binary fractions and metallicity continues into the domain of halo stars. Particularly, it underscores the influence of metallicity during the early stages of the Milky Way’s formation and in the accreted dwarf galaxies. The higher prevalence of binary systems in metal-poor environments suggests a fundamental aspect of binary formation in the early universe, where lower metallicity was more common. Such subtleties in the data not only enhance our understanding of the binary star dynamics but also pose intriguing questions about the historical processes that have shaped these populations.

Building on this, the impact of primary mass on binary fraction adds another layer of complexity to our analysis. Since the binary fraction depends on the primary mass of the system (Raghavan et al. 2010; Moe et al. 2019), exploring this relationship could bias our results. To address this, we conducted a comparative analysis of the Abs. G distributions between the two halo samples, hypothesizing that it could serve as a proxy for primary mass. Our findings revealed that the distributions in Abs. G between the accreted and in situ samples were comparable, suggesting that despite the challenges posed by metallicity, the primary mass, as inferred from Abs. G, does not significantly differentiate the binary fractions of the two halo samples. Similarly, in the analysis of the Gaia-APOGEE metallicity bins, we intentionally restricted our sample to narrow region of effective temperature |$4600{\!-\!}6200$| K and stellar mass within the |$0.5{\!-\!}1.2\,{\rm M}_{\odot }$| range to decrease possible biases with primary mass. However, this in any case would not change dramatically, nor reverse, the observed trend of higher binary fraction towards metal-poor stars. In fact, as Moe et al. (2019) pointed out, as metal-poor stars tend to be of later type, the trend seen of higher binary fraction in metal-poor stars should become more significant.

The latest discoveries of dormant compact objects, such as Gaia BH3 (Gaia Collaboration 2024) and a list of Neutron star candidates (El-Badry et al. 2024), have highlighted the critical importance of exploring binary fractions in metal-poor regions, particularly among halo stars. These findings suggest that metal-poor environments, often associated with older galactic structures, may harbour a significant number of these elusive compact remnants. The presence of such objects in binary systems can provide valuable insights into the evolutionary paths that lead to their formation.

Future work could aim to expand the halo sample utilized in this study to obtain more accurate estimates of the mean values of binary fractions and their uncertainties. For instance, we limited our analysis to stars with |$G_{\mathrm{RVS}}$|<12 to ensure the reliability and quality of the RV data. None the less, the Gaia RVs are available to stars up-to |$G_{\mathrm{RVS}}$| = 14. Another potential strategy involves leveraging the RV scatter, incorporating data from various ground-based spectroscopic surveys such as LAMOST (Cui et al. 2012), GALAH (De Silva et al. 2015), and APOGEE (Majewski et al. 2017). This approach would allow us to extend our sample, including a broader range of stars suspected to be in binary systems.

In terms of extending our methodical approach, future works might also exploit our Bayesian model to estimate the likelihood of individual Gaia sources being affiliated as single stars or as binary stars. This refined focus could significantly enhance our understanding of the distinct characteristics and behaviours of specific systems allowing further study of their specific nature.

ACKNOWLEDGEMENTS

We would like to thank the reviewer for their helpful comments and suggestions. DB acknowledges the support of the Blavatnik family, the British council, and the British Friends of the Hebrew University (BFHU) as part of the Blavatnik Cambridge Fellowship. DB also extends his gratitude to Didier Queloz for his hospitality and encouragement. We wish to thank Shay Zucker for a helpful discussion on the Bayesian methodology and Jo Bovy for assisting with the definition of the Galpy’s galactic potential. This work has also made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Funding for the Sloan Digital Sky Survey IV (SDSS IV) has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss4.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics|Harvard & Smithsonian (CfA), the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA, Heidelberg), Max-Planck-Institut für Astrophysik (MPA, Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. This research also made use of TOPCAT (Taylor 2005), an interactive graphical viewer and editor for tabular data.

This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

This research has made use of the SIMBAD data base, CDS, Strasbourg Astronomical Observatory, France.

DATA AVAILABILITY

The research presented in this article predominantly relies on data that is publicly available and accessible online through the Gaia DR3 and APOGEE archives.

Footnotes

REFERENCES