© The Authors 2023

1. Introduction

The majority of Local Group (LG) dwarf galaxies (DGs) are pressure-supported (e.g., Wheeler et al. 2017) and devoid of HI gas (e.g., Putman et al. 2021), which means that their dynamical mass is inferred from measurements of the velocity dispersion of their stellar component. While some first measurements of the velocity dispersion on the plane of the sky exist for an handful of DGs satellites of the Milky Way (MW; Massari et al. 2018, 2020; del Pino et al. 2022), the vast majority of dynamical mass estimates of LG DGs rely on the velocity dispersion measured from line-of-sight (LoS) velocities, σ_los. Since the first such study by Aaronson (1983) for the MW satellite Draco, the LoS velocity dispersion has been measured for more than 80 LG DGs, with values ranging from ${\sigma }_{\mathrm{los}}=0.9_{-0.5}^{+0.6}$ km s⁻¹ for Tucana III to 35 ± 5 km s⁻¹ for NGC 205 for systems, where σ_los could be resolved (see recent review by Battaglia & Nipoti 2022, and references therein). A result that has been consolidated over the years states, contrary to globular clusters, the σ_los of LG DGs is significantly higher than the values that one expects from the contribution of baryons, with dynamical mass-to-light ratios (M_1/2/L_1/2)_dyn ≳ 10 in solar units¹ within the 3D half-light radius for the vast majority of these galaxies. In addition, the dynamical mass-to-light-ratio is inversely proportional to the DG luminosity (e.g., Mateo 1998; McConnachie 2012; Walker 2013; Strigari 2013; Battaglia & Nipoti 2022), with (M_1/2/L_1/2)_dyn ≳ 1000 for the so-called ultra-faint dwarfs (UFDs)². In standard Newtonian dynamics, these high mass-to-light-ratios are explained by the fact that dwarf galaxies are embedded in dark matter halos several order of magnitudes more massive than the luminous component. It is still possible (particularly for the faintest systems) that a non-negligible fraction of the velocity dispersion measured in DGs is a consequence of other effects, one of them being the presence of unresolved binary stars (e.g., see early studies by Hargreaves et al. 1996; Olszewski et al. 1996).

It has been shown that the effect of unresolved binary stars is minimal for galaxies with high intrinsic velocity dispersions (Hargreaves et al. 1996; Minor 2013), namely, in the regime of “classical” DGs. For example, a binary fraction equal to 1 (under the assumption of period, mass and eccentricity distributions as in the solar neighborhood) would inflate an intrinsic LoS velocity dispersion = 8 km s⁻¹ to 9 km s⁻¹ (Spencer et al. 2017). However for galaxies with low intrinsic dispersions, ≲2 km s⁻¹, the same study shows that the measured σ_los can be 1.5−4 times higher than the intrinsic dispersion already with a modest binary fraction of 0.3.

The value for which the measured σ_los, and, hence, the inferred dynamical mass-to-light ratio, inflated by binary stars, is strongly dependent on the fraction of stars located in binary systems, as has also been concluded by Pianta et al. (2022)³. In particular, if the fraction of binaries with periods lower than 10 yr exceeds 5%, there is a non-negligible (> 5%) chance that the values of σ_los from single-epoch observations might be completely due to binary systems for, for instance, Segue 1, Segue 2, Willman 1, Bootes II, Leo IV, Leo V, and Hercules (McConnachie & Côté 2010). Thus far, in studies of UFDs for which multi-epoch spectroscopic measurements have been acquired, removing (or accounting for) stars showing obvious velocity variations reduces the value of the measured σ_los (e.g., Martinez et al. 2011; Ji et al. 2016; Kirby et al. 2017; Minor et al. 2019; Buttry et al. 2022), even though it is not enough to bring the dynamical mass-to-light ratio in line with expectations from the baryonic component only. Nonetheless, it remains of interest to determine how strongly the inferred dynamical mass-to-light ratios of individual UFDs might be inflated by unresolved binary stars.

Unlike the Solar neighborhood, the properties and the fraction of binary stars is generally unknown in LG dwarf galaxies. In the Milky Way, Duquennoy & Mayor (1991) found an overall binary fraction of ∼2/3 for stars near the middle of the main sequence, with the binary fraction decreasing from f ≃ 0.5 − 0.6 for very metal-poor stars down to f ≃ 0.1 for the stars with super-solar metallicities (Raghavan et al. 2010; Badenes et al. 2018; Moe et al. 2019). Moreover, Hettinger et al. (2015) found that disk stars, mainly metal-rich, are 30% more likely to belong to a short-period (< 12 days) binary system than metal-poor halo stars.

Regarding dwarf galaxies, the binary fraction has been investigated in a few of the Milky Way satellite galaxies, mostly in the so-called “classicals” dSphs (e.g., Minor 2013; Spencer et al. 2017, 2018), but also in some UFD (e.g., Martinez et al. 2011; Minor et al. 2019). Typically, these studies have to make the assumption that properties of the population of binary stars in DGs, such as their period and mass ratio distributions, are the same as in the solar neighborhood. What makes it particularly challenging is that this type of explorations necessitate accurate and precise LoS velocity determinations from multi-epoch observations of a large number of stars, as pointed out by Martinez et al. (2011), for instance. Minor (2013) found (with the notable exception of Carina) that three “classical” MW dwarf galaxies (Fornax, Sculptor, and Sextans) have a uniform binary fraction of f ≃ 0.5, similar to the value found for the metal-poor stars in the Milky Way. However, these results were recently contested by Spencer et al. (2018, hereafter S18), who revisited the binary fraction in these galaxies plus three others (Ursa Minor, Draco and Leo II). They found that these seven dwarf galaxies cover a wide range of binary fractions; they rejected at a 99% level the universality of the binary fraction in dwarf galaxies (i.e., that they have the same binary fraction within a range of width 0.2), provided that they share similar period distributions; or conversely, that if they share similar binary fractions, then their period distributions need to be significantly different. Bonidie et al. (2022) examined the velocity variation of stars across different epochs in the Sagittarius DG and compared them with a sample of stars from the MW of similar stellar atmospheric parameters. These authors found that Sagittarius stars display larger velocity variations than the comparison Milky Way stars, so they ruled out all the possible explanations – leaving just the notion of a larger binary fraction for close binaries.

If the results by Spencer et al. (2018) were confirmed, this would exclude the possibility of using the binary stellar population of dSphs as templates for those in UFDs; however, constraining its properties is even more challenging due to the scarcity of target stars bright enough to be followed up with current facilities. In this paper, we revisit the determination of the binary fraction in these galaxies. As aptly highlighted in Minor (2013), knowing the properties of binary stars in MW DGs is also important in constraining star formation theories in galactic environments of low metallicity and covering a variety of star formation histories. Moreover, since DGs are more diffuse than globular clusters, they present a lower chance of stellar encounters, making DGs a more controlled environment in which to test stellar formation theories, especially since formation models predict a binary fraction ranging from f = 0.05 (Hurley et al. 2007) to f = 1.0 (Ivanova et al. 2005) for them.

In this work, we first focus on the Sculptor dwarf spheroidal (dSph), presenting results from a homogeneous spectroscopic dataset for 96 red giant branch (RGB) stars obtained at the Very Large Telescope (VLT) with the multi-object spectrograph Fiber Large Array Multi Element Spectrograph (FLAMES). This dataset was based on the same instrument set-up between 2003 and 2015, with 9 to 12 epochs per star. Given the results obtained for this galaxy, we propose some modifications to Bayesian methods existing in the literature (Spencer et al. 2018), so as to obtain more robust uncertainties. Using this new method, we re-analize the binary fraction of the six other DGs studied in S18 (Draco, Ursa Minor, Sextans, Carina, Leo II, and Fornax).

The paper is organized as follows. In Sect. 2, we introduce the datasets used in this work, including the new multi-epoch VLT/FLAMES dataset for Sculptor. In Sect. 3, we describe the methodology employed and present the changes with respect to S18. Our results are presented in Sect. 4, where we first determine the binary fraction in Sculptor (Sect. 4.1) and then those of the other six MW dSphs (Sect. 4.2). In Sect. 5, we present an outlook for what future observations with next-generations multi-object spectrographs such as WEAVE will contribute to studies of the population of binary stars in MW dSphs (Dalton et al. 2012; Jin et al. 2023). Finally, we draw our conclusions in Sect. 6. In Appendix A, we discuss whether there might be a common binary fraction that would cover all the DGs.

2. Data

Here, we present a new set of FLAMES/GIRAFFE observations for the Sculptor dSph (Sect. 2.1), followed by the literature datasets used for Draco, Ursa Minor, Leo II, Carina, Fornax, and Sextans (hereafter, Dra, UMi, Leo II, Car, Fnx, and Sxt), as explained in Sect. 2.2.

2.1. FLAMES/GIRAFFE Sculptor dataset

Our analysis of the binary star fraction in the Sculptor dSph relies on LoS velocities of individual RGB stars measured at different times over a large time baseline. The LoS velocities are obtained from spectroscopic data acquired at the VLT with the instrument FLAMES in GIRAFFE mode. FLAMES is perfectly suited for the objective of this work as it is known to be a very stable spectrograph (Pasquini et al. 2004). The section is organized as follows. In Sect. 2.1.1, we describe the main characteristics of the data and briefly comment on the data-reduction process and the determination of LoS velocities. In Sect. 2.1.2, we perform an analysis of the possible zero-point velocity offsets between the different observations of the same sources. In Sect. 2.1.3, we analyze the robustness of the velocity uncertainties.

2.1.1. Observations, data-reduction, and determination of velocities

Our targets are 96 stars with Gaia early Data Release 3 (eDR3) (Gaia Collaboration 2016, 2021) G-magnitudes, approximately between 16.5 and 19 and located within 25 arcmin from Sculptor center (see Fig. 1). These were selected among the likely Sculptor members from the datasets presented in Tolstoy et al. (2004), Battaglia et al. (2008, 2010)⁴ and followed up on as part of an ESO monitoring program over two years (ID: 593.D-0309; PI: Battaglia). Specifically, the monitoring consisted of eight repeated observations of 2700 s, between 2014 and 2015, with the same instrument set-up and fiber configuration, in order to reduce the source of errors on the repeated measurements of the LoS velocities. The grating used was LR8, covering the region of the near-infrared Ca II Triplet (CaT). The choice of placing the pointing at a location displaced from Sculptor center was due to the wish of probing a region with a mean metallicity more akin to that of ultra-faint dwarf galaxies. In the monitoring programme, we aimed for observations spaced approximately logarithmically over time. For the actual observation dates, see Table 1.

Fig. 1. Distribution of Gaia eDR3 sources (in gray) within a 50 arcmin radius around the central coordinates of the Sculptor dSph, with overlaid the 96 targets of the FLAMES-GIRAFFE multi-epoch observations (large dots). Left panel: distribution on the sky. Right panel: distribution on the color-magnitude diagram. The color indicates the number of repeated observations (see legend).

In the earlier programs the targets were drawn from, the spectra for these 96 stars were acquired over multiple pointings, between 2003 and 2007 – again, with the same instrument and grating, but with a different fiber configuration. In these nights not all the same 96 stars were observed, but each of them has at least one extra observation in the 2003–2007 period. Therefore, the sample consists of 96 stars, with at least nine repeated measurements, with a minimum time baseline of about 8 yr (2940 days). Figure 2 summarizes the number of stars for a given number of repeated measurements and time-baseline. Table 1 lists the dates of the observation and the number of stars proceeding from those exposures.

Fig. 2. Distribution of time baseline and number of repeated observations for the 96 stars that are part of our Sculptor dataset.

The approach adopted for the data reduction, extraction of the 1D spectra, sky-subtraction, and determination of the velocities is presented in Tolstoy et al. (2023), to which we refer for more details. Of particular relevance here is the set of LoS velocities and uncertainties adopted. In Tolstoy et al. two methods for determining LoS velocities were used: (1) cross-correlation with a template with zero continuum and three Gaussians mimicking the CaT lines in absorption, as in Battaglia et al. (2008), with the uncertainties being derived from the scatter of the velocities returned by each of the 3 Gaussians; (2) a maximum-likelihood-method where 3 parameters (v_los, Δv_los and best fit Gaussian width for blurring the delta function of the CaT template) are explored to obtain the most likely combination. In order to decide which set to adopt, we did our own error analysis for the 96 stars of the monitoring program, since it is a key parameter in our work. Here, we adopted the cross-correlation set, since it has more reliable errors.

As a confirmation that the targets of the monitoring program are members to Sculptor, we checked that their average LoS velocity ranges between 89 km s⁻¹ and 137 km s⁻¹, which is well within 3-σ from the systemic LoS velocity of Sculptor (the systemic velocity being v_los, sys = 110 ± 0.5 km s⁻¹ and the LoS velocity dispersion σ_los, sys = 10.1 ± 0.3 km s⁻¹, see Battaglia et al. 2008). These values are in agreement with what is found in Tolstoy et al. (2023) with an updated version of this dataset that features a total of 1604 probable member stars with velocity measurements, v_los, sys = 111.20 ± 0.25 km s⁻¹.

The monitoring program was designed to deliver spectra with intermediate-to-high signal-to-noise ratio (S/N) and, therefore, precise LoS velocities. The minimum S/N in the spectra observed between 2014 and 2015 is around 20; for the spectra from the previous nights, 2003–2007, the minimum S/N is around 5.6. These spectra produce a mean velocity error of ∼1 km s⁻¹. In the final dataset, namely, taking into account all observations, the individual spectra display a mean S/N of 56. In Fig. 3, we show three spectra of the same star at different dates, in which we can see velocity variation.

Fig. 3. Spectra of an example Sculptor member star, observed during three different nights, as indicated in the legend. The top panel shows the whole spectra. The bottom panels show a zoom-in view of the CaT lines, with the x-axis indicating the velocity shift with respect to the rest wavelength (indicated with dashed lines) of a given CaT line.

2.1.2. Correction for zero-point offsets

Even though the observations were carried out with the same instrument and set-up, systematic differences in the v_los between the exposures taken during different nights can still exist. Those systematic errors would increase the differences between the v_los for a given star, potentially inducing an inflated inferred binary fraction.

Here, we compare the LoS velocities of each star j obtained in each exposure i, v_i, j, with the LoS velocities measured for the same star during a reference night, v_ref, j. For simplicity, hereafter, we drop the subscript j. As reference night, we adopt that with the best S/N and G-magnitude relation, which was on 2014-08-19. An example of the distribution of the differences, v_ref − v_i, is shown in the upper left part of Fig. 4, where we present the case only for three nights in order to better visualize the results. The typical zero-point shifts, measured as the mean values of v_ref − v_i (listed in the second column of Table 1), are typically rather small, around 0.18 km s⁻¹ for the nights of the monitoring program; nonetheless, we subtracted those shifts from the velocities of the stars observed in a given night. Examples of the distributions corrected for these offsets are given in the upper right panel of Fig. 4: now all the distributions are centered on 0 km s⁻¹. Since we do not have enough statistics for the 2003-09-30 night (only 1 out of our 96 stars was observed that night), those velocities have not been corrected. The other case with low statistics is the night of 2007-09-16. For this one, we decided to correct by offsets. However, this is not a problem since we checked that removing any of those nights, does not impact the results at all.

Fig. 4. Distributions of velocity differences and α. Upper panels: histograms of the difference between velocities for the different nights with respect to the reference night, as indicated in the legend, before the zero-point offset correction (left) and after the correction (right). Lower panels: same as the upper panels but for α instead of the velocity difference. Only three nights have been shown for ease of visualization, but all nights have been corrected.

2.1.3. Robustness of the velocity uncertainties

In order to check whether the uncertainties on the v_los measurements are robust, we examined the distribution of the differences between the v_los of each star in a given night and the reference night, v_ref − v_i, weighted by the respective velocity uncertainties:

$$\begin{array}{c}\hfill \alpha =\frac{v_{\mathrm{ref}}-v_i}{\sqrt{{\left(\mathrm{\Delta }{v}_{\mathrm{ref}}\right)}^2+{\left(\mathrm{\Delta }{v}_i\right)}^2}},\end{array}$$(1)

where Δv_ref and Δv_i are the velocity uncertainties of a given star j in the nights ref and i. Table 1 shows the scaled median absolute deviation (sMAD)⁵ for the parameter α for each given night, after applying the offsets’ correction (Col. 4). We prefer the sMAD to the standard deviation as the former is more robust to the presence of outliers. In the lower part of Fig. 4, we can see examples of the distribution of α before and after the correction. For a set of measurements where the only source of errors are statistical ones, we expect the distribution of α to be a Gaussian centered on 0 and with sMAD of 1; values of sMAD lower or larger than 1 would indicate over- or under-estimated velocity uncertainties, respectively. In Table 1, we may notice that the sMAD of the distribution of α (Col. 4) for the nights between 2014 and 2015, where the same exact fiber configuration was used, is systematically lower than 1, which would point to the uncertainties being overestimated. Dividing them by a factor 1.15 leads to sMAD values centered around 1 (Table 1, Col. 5), so we adopted this correction factor for our uncertainties. However, after the correction, we still obtain a range of sMAD between 0.81 and 1.23. There are a couple of reasons as to why the sMAD might differ from 1 even after correcting the velocity uncertainties: (a) in this case, we are treating the dataset as if there were no binary stars and (b) we are studying a sample of 96 points as a statistical distribution.

Focusing on point b, we sampled 20 000 Gaussian distributions centered on 0 and with standard deviation 1 with 96 points and then computed their sMAD. The final distribution of the values of sMAD is a Gaussian centered on 1 and with a width of σ = 0.12. Therefore, all the values of sMAD we obtained after correcting the velocity uncertainties are within 2σ of the sMAD expected for a distribution sampled with 96 points. The results from the 2003–2007 nights have larger deviations from unity; however, in those cases, the number of stars in the exposure is lower than 96, hence it is natural to observe larger deviations⁶. Therefore, in general, our results are in line with the expectations under the conditions of the observations and we consider our velocity errors to be estimated as well, once they are corrected by a factor of 1.15.

The jitter due to the movement of material in the surface of the RGB stars is also a factor to take into account when it comes to uncertainties. This is because it can occur on short timescales, on the order of minutes (Yu et al. 2018) and, therefore, we are going to measure it as it was noise in the velocity measurements. The effect of jitter on v_los is a function of the stars’ surface gravity (e.g., Hekker et al. 2008). Therefore, we have cross-matched our dataset with APOGEE to obtain a subsample of stars with known values of surface gravity. The stars in our dataset have 2 > log(g)[cgs] > 0.2, which correspond to jitter velocities between 0.12 < v_jitter [km s⁻¹] < 1.40. On average, the jittering velocities of our sample are around 0.5 km s⁻¹, so it should already be accounted for in our errors, and in the correction factor that we have applied.

We also checked for possible trends in the velocity determination, that is, those produced by the wavelength calibration of the spectra, by searching for trends in α as a function of v_ref. In Cols. 6 and 7 from Table 1, we can see the slope m of the best linear fit of α as a function of v_ref and the dispersion coefficient, r². We find that both m and r² are close to 0 in all the cases, indicating that there is no increasing or decreasing trend of α with v_ref and no relation between both variables.

2.2. Data from the literature

In this work, we also re-analyze spectroscopic datasets used in the literature for studies of the binary stellar population of MW classical dSphs. Specifically, we use the datasets from Draco and Ursa Minor presented in S18; the one from Leo II presented in S17; and the ones from Sextans, Carina, and Fornax presented in Walker et al. (2009). All of these were used by Spencer et al. (2018), while Minor (2013) adopted only the datasets from Walker et al. (2009), available at the time.

3. Methodology

The methodology, we adopt for determining the binary fraction, f, of the dwarf galaxy stellar population largely follows that of S18, but with three main differences, as highlighted in the next sections. Broadly speaking, it consists of producing mock datasets with the same characteristics as the observational data, in terms of number of stars, time sampling, and velocity uncertainties, for different values of f, and then comparing them with the observational data.

In Sect. 3.1, we give the details of the simulation strategy, in Sect. 3.2, we recall the formula for the v_los, of a component of a binary system and the parameters it depends on, and we discuss the distributions we use for these parameters. In Sect. 3.3 we present the maximum a posteriori method (MAP) we follow to compare the simulated and the observational data. Finally, in Sect. 3.4, we highlight the differences with respect to similar methodologies used for this same objective by S17 and S18.

3.1. Simulation strategy

The following steps to produce the simulated datasets are the same as published in S18. We summarize them below.

(1) First, we chose a value for the binary fraction, f.

(2) For each observed star, we randomly labeled it as a binary or single star, according to the binary fraction in Eq. (1). In practice, this is done by randomly extracting a value from a uniform distribution between 0 and 1. If the number extracted is lower than f, the star is labelled as binary; otherwise, it is not. If it labeled as binary, we randomly assign to it a set of the seven parameters extracted from the distributions in Sect. 3.2.

(3) We computed the v_los of the star at the time of each observation. If the star is a binary, we use Eq. (3), taking into account the time variation of the true anomaly along all the observation epochs. If the star is a single star, we assign it a velocity of 0 km s⁻¹. Here we neglect the change in v_los, due to the intrinsic motion of stars within the galaxy, as it is negligible on our time baseline.

(4) Next, we computed the Gaussian deviations of the velocities with a standard deviation equal to the observational error of each observation. We added this value to the one obtained in the third step for both binaries and single stars. We would like to highlight that it is not necessary to add the systemic motion of the galaxy, since we are interested in the temporal velocity variations of each given star in the sample with respect to its velocity at a reference time.

(5) We repeated steps three and four for all the stars in the dataset.

(6) Then, for statistics, we repeated steps two to five 10 000 times. At this point, we have 10 000 different simulated mocks for each binary fraction f.

(7) Finally, we repeated steps one to six for all the binary fractions we are interested in. We sample binary fractions between 0 and 1, with a step of 0.01. As in S18, the parameter that serves as the basis for this study is:

$$\begin{array}{c}\hfill \beta =\frac{|v_i-v_j|}{\sqrt{{\left(\mathrm{\Delta }{v}_i\right)}^2+{\left(\mathrm{\Delta }{v}_j\right)}^2}},\end{array}$$(2)

where the subscripts i and j denote different observations of the same star, v is the measured velocity, and Δv the error of this measurement. The value of this parameter represents how large the velocity variation is between the different observations of the same star, weighted by the velocity uncertainties. We compute this parameter for both the observed and simulated data (β_obs and β_sim, respectively). For a star with n observations, there will be $\frac{n(n-1)}{2}$ different values of β. In our particular case, we have 96 stars observed between 9 and 12 times. This leads us to 3986 different values of β. In previous works, the observations produced 2015 β for Draco (Spencer et al. 2018), 1278 for Ursa Minor (Spencer et al. 2018), and 723 for Leo II (Spencer et al. 2017), far less than with our observations. So, we can see that we have not only significantly more data, but also the uncertainties in our velocity measurements are better, as we saw in Sect. 2.1.

3.2. Parameters distributions

Assuming that a binary system behaves as a two point mass distribution with only gravitational interaction, v_los, of one of the components, with respect to the system’s center of mass, can be derived analytically (e.g., see Green 1985):

$$\begin{array}{c}\hfill v_{\mathrm{los}}=\frac{qsin(i)}{\sqrt{1-e^2}}{\left(\frac{2\pi G{m}_1}{P{(1+q)}^2}\right)}^{1/3}(cos(\theta +\omega )+ecos(\omega )).\end{array}$$(3)

The equation contains four intrinsic (q, m₁, P, and e) and three extrinsic (i, ω, and θ) parameters. Among the intrinsic ones: m₁ is the mass of the primary star; $q=\frac{m_2}{m_1}$ is the ratio between the mass of the secondary star, m₂, and of the primary star m₁; e, the eccentricity of the orbit around the center of mass; and P the period of the orbit. Of the extrinsic parameters: ω is the argument of the periastron, the angle corresponding to the point of the orbit that is closest to the center of mass; i is the inclination, that is, the angle between the plane determined by the binary system’s orbit and the plane perpendicular to the LoS; and θ is the true anomaly, the only parameter that changes with time, as it indicates the position of the star along the orbit.

We adopted the same distributions used in previous studies of binary fractions of DGs (Spencer et al. 2017, 2018; Minor 2013), so that we can directly compare our results with theirs. We adopted the following distributions (see also Fig. 5). For the mass of the primary star, m₁, we used a fixed value of m₁ = 0.8 M_⊙, since this is the typical mass expected for the majority of our target old and metal-poor RGB stars (see e.g., Martig et al. 2016)⁷. For the mass ratio, q, we used a Gaussian distribution centered on μ_q = 0.23 and with dispersion σ_q = 0.42 (see upper left panel of Fig. 5), as found by Kroupa et al. (1990). It has been shown to fit correctly the distribution of mass ratios for G-dwarfs belonging to binary systems in the solar neighborhood (Duquennoy & Mayor 1991). However, we also tested a uniform distribution, as Kroupa et al. (1990) also showed that the distribution of mass ratios is flattened for binaries with short periods. This is assuming the uniform distribution yields very similar results to the Gaussian one. As for P, most studies about binary systems in the solar neighbourhood have shown that the observed period distribution is a log-normal function (Duquennoy & Mayor 1991; Fischer & Marcy 1992; Marks & Kroupa 2011):

$$\begin{array}{c}\hfill f(logP[\mathrm{d}])exp\left(\frac{-{(logP[\mathrm{d}]-\overline{logP[\mathrm{d}]})}^2}{2{\sigma }_{logP[\mathrm{d}]}^2}\right),\end{array}$$(4)

Fig. 5. Probability distribution functions for the parameters defining the properties of the binary star populations. The black lines show the theoretical distributions while the colored histograms show the results of our sampling. Upper panels: period distribution for the two limits in mass-ratios, as indicated by the legend, assuming a constant minimum orbital distance amin = 0.21 AU (left) or the taking into account the RLOF when computing amin (right). The limits for each mass-ratio are plotted with vertical lines. Middle left panel: distributions adopted for the mass ratio (green: GS distribution from Duquennoy & Mayor 1991; salmon: uniform distribution). Middle right panel: distribution of inclinations. Lower left panel: eccentricity distribution; the differences between the uniform theoretical distribution and the one sampled are due to the dependence of e on the semi-minor axis and therefore, the orbital period, as explained in the main text. Lower right panel: distribution of the argument of the periastron.

where $\overline{logP}$ is the mean of the distribution, with the period in units of days, and σ_{log P} is its width.

Even though previous works (Minor 2013; Spencer et al. 2017, 2018) have shown that the choice of the mean of the distribution and σ_{log P} is crucial for the determination of the binary fraction in dwarf galaxies, the available datasets do not yet allow for such parameters in dwarf galaxies to be inferred, and strong degeneracies are present between the values of these parameters and the binary fraction that best reproduces the observational data. Throughout the article, we adopt the typical baseline case of σ_{log P} = 2.3 and $\overline{logP[\mathrm{d}]}=4.8(\approx 180\mathrm{yr})$ as it has been done previously in the literature (Minor 2013; Spencer et al. 2017, 2018). These values have been obtained by fitting the period distribution of stars between the types F7 and G9 in the solar neighborhood (Duquennoy & Mayor 1991). In Sect. 4.1.2, we will use other realistic values for $\overline{logP}$ and σ_{log P}, to see whether some are preferred by the data. For σ_{log P}, we will use values between 1.7 and 2.9. For the mean of the distribution, we will use as lower limit a value of $\overline{logP[\mathrm{d}]}=3.5$, obtained from the analysis of K-type stars and M-dwarf stars in binary systems by Fischer & Marcy (1992) (corresponds to separations of 3 AU for 0.8 M_⊙ RGB stars). As an upper limit, we adopted $\overline{logP[\mathrm{d}]}=5.8$, which was determined by S18, they fit the period distribution of the evolved binaries population synthesis for dwarf irregular galaxies made by Marks & Kroupa (2011).

The period distribution will be limited at the lower and upper ends by the minimum and maximum semi-major axis distances, a, expected for the type of targets we are considering, RGB stars, and their galactic environment. Using the third law of Kepler, $P^2=\frac{4{\pi }^2}{G{m}_1(1+q)}{a}^3$; this implies a dependence also on q. Typically, the minimum value used for the semi-minor axis is a_min = 0.21 AU, that this is the typical radius of RGB stars; it can be obtained by assuming a mass of 0.8 M_⊙ and surface gravity of 10 cm s⁻². Considering the extreme cases of q = 0.1 and q = 1, the minimum period is then log(P_min[d]) = 1.57 for q = 0.1, and log(P_min[d]) = 1.44 for q = 1.

However, when considering the Roche lobe overflow (RLOF), the less dense star of the binary system (usually the RGB star) will lose material before both components are within 0.21 AU and the equilibrium of the system will be broken. At this moment, either the envelope of the giant star is lost and its core and the companion keep orbiting around the center of mass, or a supernovae type Ia is produced. Simulations of this type of interaction for a 0.88 M_⊙ star and a compact companion of 0.6 M_⊙ (similar conditions of what we typically expect to have) have shown that the typical time at which this happens is only of some tens of years (Reichardt et al. 2019). Therefore, we do not expect to observe binary stars in the overflow phase and therefore orbiting closer than the Roche limit.

In the approximation by Eggleton (1983), the Roche lobe, r_R, can be calculated as:

$$\begin{array}{c}\hfill \frac{r_{\mathrm{R}}}{a}=\frac{0.49q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})}·\end{array}$$(5)

The condition for not having overflow is that r_R > 0.21 AU; this translates into a minimum orbital distance of a_min > 1.02 AU for q = 0.1 and a_min > 0.56 AU for q = 1. Those values correspond to minimum periods of log(P_min[d]) = 2.60 for q = 0.1 and log(P_min[d]) = 2.07 for q = 1. We use the limits corresponding to the RLOF approximation in the main part of our analysis⁸.

For the maximum value of the semi-major axis, a_max, we require that the two components of the binary systems are gravitationally bound at all times. We will consider that the two stars are no longer bound when an interaction with another star is produced, that is to say, when a is larger than the mean free path of a star in the galaxy. The mean free path can be obtained as (πσ_lostλ)^−1/2, where σ_los is the galaxy’s LoS velocity dispersion, t is the average age of the stars, and λ is the volume number density of stars. In order to obtain the volume number density of stars, we assume the luminosity density listed in Table 2, together with the luminosity of the primary being L/L_⊙ = (M/M_⊙)⁴ (relation that works well in our range of masses; Reid 1987), and the fact that with m₁ = 0.8 M_⊙ and with the q distributions examined, each star has an average mass of around 0.4 M_⊙. Once we have the mean free path, that is, a_max (listed in Table 2, together with the parameters used to compute it), we can use again the third law of Kepler to compute the value of P_max, this time it will be different for each galaxy and will depend on q again. The period distributions are shown in the upper right panel of Fig. 5, and their lower and upper limits summarized in Table 3.

As mentioned previously, the period of binary systems with our target stars as a primary is unlikely to be shorter than log(P_min[d]) = 2.7, which corresponds to 500 days. For periods longer than this value, the eccentricity distribution is expected to be rather uniform (Duquennoy & Mayor 1991; Raghavan et al. 2010). However, in our sampling of a uniform distribution, we will take into account that certain combinations of P and q can make the stars collide for very elliptical orbits. This upper limit in the allowed eccentricity is, as a function of the semi-minor axis, e_max = 1 − (a_min/a). Therefore, the order we have followed to sample the e distribution (shown in the lower left panel of Fig. 5) is to simulate random values of q according to the assumed distribution of this parameter, then for P, and, finally, for e, taking the various dependencies into account.

The inclination (i) describes the angle between the normal to the plane of the orbit and the LoS to the observer. After some trigonometry, one finds that the distribution of i is f(i) = sin(i) (see middle left panel of Fig. 5), with a face-on orbit having i = 0 rad and an to edge-on orbit $i=\frac{\pi }{2}$ rad.

The argument of the periastron (ω) is the angle between the ascending line of nodes of the elliptical trajectory of the star and the periastron. Since this orientation is random, the distribution of this parameter is a uniform distribution. We can see the distribution in the lower panel of Fig. 5.

The true anomaly (θ) is the angle between the direction of periastron and the position of the star, from the main focus of the ellipse. For a star with multiple observations, we need to take into account that this parameter changes with time. In order to sample θ, we made use of an auxiliary parameter called mean anomaly (M), that is to say, the fraction of an elliptical orbit’s period that has elapsed in a time, t. So, M = M₀ + kt, where M₀ is the initial mean anomaly and k is a constant, due to the second law of Kepler, as a star that has an elliptical orbit sweeps out equal areas in equal times. Thus, the parameter k can be computed by using the third law of Kepler $k=\frac{2\pi }{P}$. We sample M₀ with a uniform distribution. Now, we can relate the mean anomaly of a star with its eccentric anomaly E with the equation M = E − e ⋅ sin(E), where e is the eccentricity. We solve the equation numerically for E using the function solve of the python package SciPy. Then, the true anomaly can finally be obtained with the relation $cos(\theta )=\frac{cos(E)-e}{1-e·cos(E)}$. All the distributions above have been sampled by implementing a Monte Carlo rejection method.

3.3. Maximum a posteriori method

Now, for each f, we are going to use the 10 000 simulations to build a likelihood function that describes how well the observational data can be described in terms of the properties of the binary stars population. Then, we evaluate this function using our observational data to obtain the posterior distribution function (PDF) of a given binary fraction.

Recalling Bayes theorem, we have that:

$$\begin{array}{c}\hfill P\left(f|D,M\right)=\frac{P\left(D|f,M\right)P\left(f|M\right)}{P\left(D|M\right)}·\end{array}$$(6)

In the above, P(f|M) is the prior probability of a binary fraction f using a model, M, which we assume is a uniform distribution between 0 and 1, since we do not have any prior knowledge on the binary fraction. P(D|f,M) is the likelihood function, it represents the probability that our observational data is well described by a fraction f using a theoretical model, M; this is the function that we have to build using the parameter β and the 10 000 simulated mocks for every fraction, f. The next term, P(D|M), is a normalization factor that we will choose so that the integral of the PDF (P(f|D,M)) is 1. Therefore, since P(f|M) is uniform, the term P(f|D,M) that represents the probability that a fraction, f, fit well our observational data D for a theoretical model, M, is directly proportional to the likelihood P(D|f,M). Then, we only need the likelihood function to obtain the PDF on f.

To build our likelihood function, we will be comparing the cumulative distribution function (CDF) of β produced by the observations and by the simulated datasets (see Fig. 6), which we refer to as N_obs(β ≤ ℬ) and N_sim(β ≤ ℬ), respectively:

$$\begin{array}{c}\hfill P\left(f|D,M\right)\propto {{\displaystyle (\sum _{i=1}^{10000}\int |N_{\mathrm{obs}}(\beta \le \mathcal{B})-N_{\mathrm{sim}}^{f,i}(\beta \le \mathcal{B})|\mathrm{d}\beta )}}^{-\gamma },\end{array}$$(7)

Fig. 6. Cumulative distribution function of β for Sculptor. The black lines show the CDF obtained for Sculptor from the FLAMES dataset; the CDFs from 750 examples of mock datasets for a binary fraction of 0, 0.55 and 1 are shown in green, blue, and salmon respectively. As there is some overlap between the cases of f = 0.55 and f = 1, the median cumulative curve is highlighted in each case with brighter colors and a thicker line-width. Inset: Example for just one simulation, where in gray we highlight the area between the CDF of the simulated case and that for the Sculptor observations, that is what we use when computing the likelihood function.

for all ℬ positive. Hereafter, for simplicity, we write N_sim(β) and N_obs(β) instead of N_sim(β ≤ ℬ) and N_obs(β ≤ ℬ).

In Eq. (7), we are considering that the smaller the area between the N_obs(β) and N_sim(βB) curves, the better the observations are reproduced by the datasets simulated for a given binary fraction, f. The sum is over the 10 000 simulations for each f. The integral is performed numerically by summing up the difference between the 3986 sorted values of β along the y-axis in Fig. 6.

In order to turn the outcome of Eq. (7) into an actual probability distribution, we made use of the definite-positive (free) parameter γ. In practice, we first found the value of f maximizing Eq. (7), with γ = 1; let’s call it f_max. We then simulated 100 mock datasets with a binary fraction equal to f_max, using the same simulation strategy explained in Sect. 3.1. Afterwards, we varied the parameter γ from 1 to 10 with 0.5 steps (range and precision that we have tested work well), analyzed them as done for the observations, and for each we normalized the PDF. Next, we computed which value of γ is needed so that 68.2% of the PDFs include f_max within 1σ, 95% within 2σ and 99% within 3σ. This way, we make sure that the uncertainties provided with these PDFs are robust. We note that since γ is only an exponent, it will change only the width of the PDF – not the value of f_max.

This procedure has to be repeated for every dataset we use (see Sects. 2.1 and 2.2), since the quality of the data affects the required γ. Once we have the correct value for γ, we only have to compute the γ power of our posterior. Finally, we normalized this posterior and we obtained a PDF that peaks at the maximum of the likelihood and that provides well-defined uncertainties. As we see in the next sections, given the lack of this re-normalization step through the parameter γ, the uncertainties on f can be seriously underestimated.

3.4. Changes with respect to S18 methodology

As discussed, our methodology relies closely on the approach by S18, except for three main changes, which we implemented for our analysis of the Sculptor FLAMES dataset, as well as for a re-analysis of literature data for Draco, Ursa Minor, Leo II, Carina, Sextans, and Fornax. We detail these changes here.

First, S18 bin the values of β and analyzes the distribution of β in those bins individually by fitting them with Poisson distributions; in our case, we are using the individual values of β and analyzing their CDF. Second, once S18 computes their PDF, the width of the distribution is not corrected and, therefore, the uncertainties quoted result into formal uncertainties. We note however that the authors provide results from mock simulations (see their Fig. 7), from which it can be appreciated that the formal precision on f has a significantly smaller value than what was suggested by the analysis of mock datasets, even though they do not discuss this implication. In our case, our PDF is corrected so that the quoted 1-, 2-, and 3-σ uncertainties encompass the 68%, 95%, and 99% of the mock simulations, respectively. Third, S18 uses a constant value of the minimum orbital distance, fixing it to the average radius of a RGB star: 0.21 AU. In our case, we have also considered the RLOF, so that our minimum orbital distance is the one that makes the Roche-lobe radius larger than 0.21 AU, so that there is no overflow between the components of the binary system.

4. Results

In Sect. 4.1, we analyze the results from the application of our method to the Sculptor data presented in Sect. 2.1 and test the method on mock datasets, while in Sect. 4.2, we revise our current understanding on the binary fraction of MW classical dSphs, based on literature datasets.

4.1. Results from the FLAMES/GIRAFFE Sculptor data

Figure 7 shows the PDF from the new VLT/FLAMES dataset for Sculptor, both for the case of a_min = 0.21 AU and when accounting for RLOF. We can see that the median f value returned by the PDF in these two cases is considerably different, even though the sizeable uncertainties bring them in agreement within; essentially, 1σ ($f=0.{34}_{-0.16}^{+0.18}$ in the first case and $f=0.{55}_{-0.19}^{+0.17}$ in the second case).

Fig. 7. PDFs obtained for Sculptor, using the VLT/FLAMES dataset. In salmon, we show the result assuming amin = 0.21 AU and in blue the one taking into account the RLOF, when computing amin. The median and confidence levels of 1σ and 2σ are indicated with horizontal lines.

As for distinguishing what hypothesis leads to a better representation of the data between a_min = 0.21 AU and when accounting for RLOF, we can gain information by looking in more detail into the comparison of N_obs(β) with the CDF of β produced by the most probable model. Figure 8 shows the comparison to the mean CDF produced by 10 000 simulations of the most likely model for each case (f = 0.55 for the RLOF a and f = 0.34 for 0.21 AU), separated in different ranges of β for ease of visualization. We can see that for low values of βs (left panel), the distributions produced under the two assumptions for a_min are very similar, where this is because those βs are mainly due to stars that do not belong to a binary system or that have velocities with large uncertainties. In the middle panel, we can see that now both distributions tend to the observational one, with the RLOF case displaying a slightly larger β than the 0.21 AU case at a given value on the y-axis, until around β = 4. This is the range where the long period binaries are acting. For β, from 4 to 5, we see a change in the behavior, with the red curve (the case for a_min = 0.21 AU), starting to show larger βs, with respect to the RLOF case at a given N_β. This is the regime where the short period binaries come into play. Finally, in the right panel, showing the regime of larger βs, the binaries with short periods are dominating, that is why we see a significant difference between the red and blue curves. Among the two, the curve produced with the RLOF approximation stays closer to the observational data than the curve obtained under the assumption of a_min = 0.21 AU, which departs significantly from the observed behavior. A χ² test using bins of size = 1 and using the CDF of β from the observations as the ‘expected’ distribution yields a χ² = 329 for a_min = 0.21 AU, four times larger than that for the RLOF case, χ² = 75. Based on this comparison, we conclude that the RLOF case is to be preferred with respect to fixing the minimal value of the orbital distance to just the radius of the red giant star. Therefore, from now on, all the results shown will be derived under this assumption⁹. The FLAMES data return a value of $f=0.{55}_{-0.19}^{+0.17}$ for the binary fraction of Sculptor.

Fig. 8. Cumulative distribution function of β for Sculptor. In black, we show the observational distribution, in red the one produced by averaging the result of the 10 000 simulations for the most likely model using amin = 0.21 AU and in blue for the RLOF when computing amin. The panels show different ranges for the axes in order to better visualize the results.

4.1.1. Testing the method and comparison to S18

If we concentrate on the results for the a_min = 0.21 AU case (shown in the previous section), which is the assumption adopted by S18, it would appear that our approach is less precise than used in the aforementioned study; however, as we will show here, this is just a consequence of our more robust values of the quoted uncertainties on f.

To tackle this aspect, we explore what results we would have obtained by applying the methodology by S18 to our FLAMES dataset for Sculptor, without the modifications we have introduced, and then check with mock datasets the robustness of the uncertainties on f.

As a test to check that our implementation of the S18 method is correct, we run it for Draco and Ursa Minor, using the same datasets as the authors and making the same assumptions for the binary parameter values and distributions. Figure 9 shows the resulting PDF of the binary fraction, f, and Table 4 summarizes the median of the PDF of the binary fraction together with the 1σ, 2σ, and 3σ confidence levels. We can see that our results for Draco and Ursa Minor are in excellent agreement with those obtained by S18. The relative differences for f is only on the order of 2%, which is well within the 1-σ error. Those small differences are expected since we changed the size and numbers of the bins used to derive the distribution of β¹⁰.

Fig. 9. Posterior distribution function for f obtained for Sculptor, Draco and Ursa Minor generated using the methodology by S18.

The application to the FLAMES data returns a much tighter confidence interval for f at a given confidence level (see Table 4), being 0.07 for a confidence level of 68.2%, to be compared to 0.32 for S18 and 0.40 for Minor (2013) for the same level. This is probably the consequence of various factors. First, the velocity uncertainties are larger in the previous studies, with a median velocity error of 2.1 km s⁻¹ versus 1.4 km s⁻¹ in our case. Second, even though their sample is larger than ours (190 stars vs. 96), our dataset has a larger number of repeated observations and, therefore, of βs (around 250 from Walker et al. 2009 vs. 3986). Finally, our maximum baseline is three times larger (12 yr vs. 4 yr).

Nonetheless, we note that the application of the S18 methodology to mock datasets suggests that the formal uncertainties on f are significantly underestimated. This is visible in Fig. 10, where we show the result for 100 mock datasets built mimicking the characteristics of the FLAMES Sculptor dataset and a binary fraction f = 0.6. We can see that the dispersion in the distribution of the mean of these 100 PDFs is larger than the formal uncertainties that the PDFs themselves provide. Table 5 lists the percentage of models in which the simulated binary fraction f = 0.6 is within the formal 1σ, 2σ, and 3σ uncertainties for Sculptor, Draco and Ursa Minor: in all cases, this is much lower than 68%, 95%, and 99% (for example for Sculptor it is 20%, 28%, and 42%). On the other hand, in the same table we list the percentages obtained through the application of our modified methodology and we can see that they are much closer to what expected for 1-, 2-, and 3-σ uncertainties.

Fig. 10. PDFs for 100 mock galaxies simulated reproducing the characteristics of the FLAMES dataset of Sculptor, in terms of number of stars, observations of each star, and velocity uncertainties of each observation, for a binary fraction f = 0.6. Lower panel: PDFs of the 100 mock galaxies overlaid. Upper panel: distribution of the maxima of each PDF with a Gaussian fit overlaid. In both panels, the green and orange vertical lines indicate the median of the distribution of maxima and the range that encloses the 68% of the simulations (1σ confidence level).

4.1.2. Period distribution-binary fraction relation

Now, we study how the distribution of periods affects the result for the binary fraction. In order to do that, we use our VLT/FLAMES dataset for Sculptor. We leave the mean and sigma of the period distribution as free parameters, together with the binary fraction. In particular, we explored 11 values of f between 0 and 1, 11 values of $\overline{logP[\mathrm{d}]}$ between 3.5 and 5.8 and 11 values of σ_{log P} between 1.8 and 3, all of them equally spaced. Regarding the value of γ used for computing the PDF, we just used γ = 6 since it worked so well in the previous analysis (see confidence intervals for Sculptor in Table 5). We show the posterior distribution function obtained in Fig. 11. We can see that a degeneracy between these parameters exists. In general, we can reproduce the same difference between velocities using lower binary fractions and lower mean periods or larger binary fractions and larger mean periods. This relation is also translated in the dispersion of the period distribution, since a larger dispersion increments the possibilities of producing short period binaries, and those systems matter a lot in terms of contributing to large βs, so this implies lower binary fractions. This is well in lines with the results obtained by Minor (2013). However, there are three important points that deserve remarks. First, it is not likely to have binary fractions below 0.4 for any period distribution. Second, the different values of σ_{log P} produce almost the same shape in the probability distribution. This is because we have constrained the limits of the period distribution properly and then the width of the distribution does not affect the results that much. And third, even though we also have degeneracy in $\overline{logP[\mathrm{d}]}$ for periods larger than $\overline{logP[\mathrm{d}]}=4.8$, for the lower values we can see that binary fractions between 0.4 and 0.6 are favored for all possible σ_{log P}. If we look for the most probable case, we obtained a binary fraction of f = 0.50, $\overline{logP[\mathrm{d}]}=4.2$ and ${\sigma }_{logP}^2=1.7$. Even though with the data that we have for Sculptor, this particular case is still the one that better fits the observational CDF of β, as we checked, we can see that there is still a lot of degeneracy. In particular, for f = 0.50 almost all the possible combinations of periods provided similar results. Therefore, we can not make any strong conclusion about the period distribution if we want to constrain the binary fraction simultaneously.

Fig. 11. PDFs obtained for the FLAMES Sculptor dataset when simultaneously obtaining two parameters between the binary fraction, $\overline{logP}$ and σlog P, and keeping the third one fixed to the values in the legend. From left to right, the parameter being kept fixed is the width of the period distribution, the mean of the period distribution and the binary fraction. The probabilities have been normalized so that if we sum a whole column, we obtain 1.

We note that in the work by Moe et al. (2019), which analyzed the fraction of close binaries in Solar-type stars, found it to have a dependence on metallicity. Rescaling a log-normal distribution to a given observed value of the fraction of stars with companions can account for these results into a period distribution that is not constant with metallicity, with the period distribution shifting towards shorter periods at low metallicity. It could be argued that a distribution with $\overline{logP[\mathrm{d}]}=4.0$ and ${\sigma }_{logP}^2=2.3$ is better aligned with the period distribution for metal-poor stars in Moe et al. (2019). In this case, we would find that the binary fraction of Sculptor is $f=0.{34}_{-0.19}^{+017}$, but if we look at the CDF of β, this leads to a worse fit to the data.

We have also summed the values of each probability map, to compute the PDF for a certain f, $\overline{logP}$, or σ_{log P} for the combination of the all the other parameters simulated. The results can be seen in Fig. 12. Each bin is the sum of one probability map of the ones shown in Fig. 11, and in the x-axis we have the columns shown in this same Fig. 11. We can see that binary fractions around 0.55 are favored when we combine all possibilities of $\overline{logP}$ and σ_{log P}. Concerning the period distribution, even though there is noise, we can see that the values of $\overline{logP[\mathrm{days}]}$ around 5.25 are slightly favored, this could be a hint indicating that we are going in the correct way to constrain it if we achieve better data. About σ_{log P} we cannot say anything, it is almost a flat distribution. However, this can indicate that the limits of the period distribution are well selected. With this low-precision result, we cannot properly constrain all the characteristics of the binary star population of Sculptor at once. Nevertheless, as a tentative result, it seems that a similar distribution of parameters of the binary fraction than the one we find in the solar neighborhood for G-dwarfs could work for red giants in DGs, if limited properly.

Fig. 12. PDFs obtained by summing the colormaps of Fig. 11. Here, each bin correspond to the sum of a whole colormap. In the upper panel, we show the PDF for the binary fraction summing all the distribution of periods simulated. In the lower left panel, the one for $\overline{logP}[\mathrm{days}]$ summing all possible binary fractions and σlog P. In the lower right panel, the one for σlog P summing all possible binary fractions and $\overline{logP}[\mathrm{days}]$.

4.2. Application to other MW dSphs

Analyzing data for Sculptor, Draco, Ursa Minor, Leo II, Carina, Fornax, and Sextans, S18 found that the binary star population of MW classical dSph differ in their binary fraction, period distributions or both; for example, if the period distribution does not vary across systems, the binary fractions are spread over a range of values with a width of at least 0.3−0.4. Alternatively, if these systems share the same binary fraction, then they have different period distributions.

Given the modifications we made to the methodologies used in the literature, we reassess this issue here. Figure 13 and Table 6 display the results of our analysis on literature spectroscopic datasets for Draco, Ursa Minor, Leo II, Carina, Fornax, and Sextans. We can see that the median does not reproduce well the region of maximum probability for Carina, Fornax, and Ursa Minor. That is why we also show the mode of the distribution and the highest density intervals (HDI). The HDI is, for unimodal distributions as our PDFs, the smallest region around the mode that ensure the correct value is included within a certain confidence level. In our case, the regions are defined as 1σ_hdi, 2σ_hdi and 3σ_hdi for confidence levels of 68%, 95%, and 99.7%. In general, the distributions are Gaussian-like, which is good since we expect this behavior for a PDF, with low noise. Even though the results agree with those in the literature (Minor 2013; Spencer et al. 2017, 2018), within 1σ, the uncertainties are quite large compared to the results from S17 and S18. However, as it happens in the previous case, now the uncertainties are well derived thanks to the Monte Carlo correction (see the right part of Table 5). We can also see that the new VLT/FLAMES dataset we present for Sculptor is the one that provides the smallest uncertainties.

Fig. 13. PDFs on f for the MW dSphs indicated in the legend from our re-analysis of literature datasets, using the new method proposed in this work. The median and 1σ and 2σ confidence levels are indicated with horizontal lines with ticks.

With these new results, there is no evidence for differences in the binary fraction across MW classical dSphs (see also Appendix A). This gives some hope for being able to use dSphs as “templates” for the binary population of UFDs, since these are more challenging systems for which to assemble multi-epoch data for large samples of member stars, on which to compute binary fractions. In order to test possible variations of these parameters in the population of MW DGs, we will have to wait for future datasets (e.g., Sect. 5).

5. Outlook

In the previous section, we have seen that currently available spectroscopic datasets for MW dSphs still yield large uncertainties for the binary fraction. The situation worsens even further when allowing for variations in the properties of the binary stellar population; for example, it is well known that a strong degeneracy exists between the binary fraction, as well as the mean and the width of the period distribution (Minor 2013).

The next years will see the new generation of wide-field multi-object spectrographs, such as WHT/WEAVE, VISTA/4MOST, and DESI produce exquisite datasets for stars in the Milky Way and its neighbouring galaxies. Here, we explore the improvements that specifically designed datasets could bring to this topic, taking the characteristics and expected performance of the WEAVE spectrograph as template (Jin et al. 2023) and Draco as a possible target galaxy.

We simulated the LoS velocities for seven observations spread over two years, three of them in the first three months and the rest equally spaced in time, in which the same 750 stars are targetted¹¹. To each of these stars we assign an uncertainty on the LoS velocity that depends on the star G-magnitude, as estimated from the analysis of simulated WEAVE spectra in the low resolution mode during part of one of the operational rehearsals (Jin et al. 2023). This is, errors between 0.4 km s⁻¹ for a star with G = 16.5 mag and 2.2 km s⁻¹ for one with G = 20 mag. As targets, we choose the 750 brightest stars in Draco, selected from the list of members in Battaglia et al. (2022). All of them have a probability membership above 0.9.

We then simulate 100 mock datasets with a hypothetical binary fraction, f = 0.6, and then repeat the same analysis outlined in Sect. 3. The resulting PDFs are shown in Fig. 14. We obtain that, with a γ = 4.5, 69% of the models predict the correct binary fraction within 1σ, 90% within 2σ, and 100% within 3σ. For an example PDF, we obtain a binary fraction f = 0.59 and 1σ, 2σ, and 3σ confidence levels of [0.55, 0.64], [0.47, 0.71], and [0.43, 0.77]. With seven exposures over only two years, the 1σ uncertainties are now less than half what the current data return. Should the binary fraction of MW classical dSphs cover a range of values with width 0.3−0.4, such differences will be detectable.

Fig. 14. PDFs for simulated WEAVE datasets of Draco. Left panel: PDFs for the 100 simulated WEAVE datasets of Draco (assumed binary fraction f = 0.6) obtained using the new methodology implemented in this work. In orange, we show the cases where the correct binary fraction is obtained within 1σ (69%); in brown within 2σ (21%, up to 90%) and in black those that predict it within 3σ (10%, up to a 100%). The vertical green line indicates the median of the maxima distribution. Right panel: one single PDF of the ones shown in the left panel, this is what we expect to obtain when analyzing the future observational data. Vertical lines indicate the median and intervals of confidence of 1σ, 2σ, and 3σ, as indicated in the labels.

These expected WEAVE datasets will produce much more precise β distributions. This could help us to constrain the period distribution of binaries in dwarf galaxies, using the same comparison we did in Sect. 4.1 (see Fig. 8). This would be usefull on the degeneracy issue between the binary fraction and the period distribution. It is important to notice that we are far more sensitive to the lower part of the period distribution, so this is the regime that right now one can aim to constrain.

6. Conclusions

Recently, it has been shown that the binary fraction, f, of MW classical dSphs differ from each other and in some cases deviates from that in the solar neighborhood. Besides the intrinsic interest in understanding the properties of the populations of binary stars in low density and low metallicity environments of MW dwarf galaxies, tackling their binary fraction is also important for determinations of the dark matter content (and, potentially, distribution) of systems at the low-mass end of the galaxy mass function. In fact, LG dwarf galaxies display the largest dynamical mass-to-light ratios and these inferences rely heavily on the σ_los of the stellar component, which is (more often than not) obtained from single epoch measurements.

Whilst it is unlikely that the observed values of σ_los are significantly inflated by unresolved binary stars for dwarf galaxies, with intrinsic σ_los > 4 km s⁻¹, this might be a serious concern for systems of likely much lower intrinsic, σ_los, as UFDs, with the level of inflation being strongly dependent on f.

We found that some of the uncertainties on f in the most recent works focusing on such determinations for MW dSphs were significantly underestimated. Therefore, in this article, we propose some modifications to such methodology to amend these issues, and test the performance on mock datasets.

We applied the method to a novel VLT/FLAMES spectroscopic sample for the Sculptor dSphs, that includes repeated observations of the same 96 red giant branch stars with a time baseline of 12 yr, with the same instrument, grating (and fiber set-up for the majority of the repeated observations). We also re-analyzed the multi-epoch datasets used before by Minor (2013), Spencer et al. (2018). The results are the following:

Among the most widely assumed hypotheses for the minimum orbital distance between the components of the binary systems, name equal to the typical radius of a M = 0.8 M_⊙ red giant star (a_min = 0.21 AU) or the larger cut-off as determined by the Roche lobe overflow, the VLT/FLAMES data for Sculptor are best reproduced by the case of the RLOF approximation.

The analysis of the VLT/FLAMES dataset for the Sculptor dSph returns a binary fraction of $0.{55}_{-0.19}^{+0.17}$ (1σ). A result that is in agreement with those from the literature (Minor 2013; Spencer et al. 2018). Moreover, we have tested our methodology using mock data reproducing the VLT/FLAMES dataset. For this binary fraction the code recovers the correct result 66%, 92% and 99% of the times within 1σ, 2σ, and 3σ, respectively. Therefore the uncertainties we provide are well-estimated.

If we reproduce the methodology by S18 precisely, and apply it to the VLT/FLAMES dataset, we obtain a binary fraction for Sculptor $f=0.{62}_{-0.03}^{+0.04}$ (1σ). This result is also in very good agreement with those in the literature (Minor 2013; Spencer et al. 2018), but with our formal uncertainties improved by a factor of 4−5. This is most likely due to the homogeneity of our sample, the larger number of repeated observations and the longer time baseline. However, in this case it only recover the correct binary fraction 20%, 28%, and 42% of the times within 1σ, 2σ, and 3σ, respectively. Therefore, we concluded that these uncertainties are underestimated. Even though the VLT/FLAMES dataset for Sculptor is the best one available in the literature for the purpose of studying the binary stellar population of a dwarf galaxy, it is still not sufficient for breaking the degeneracy between the binary fraction and the mean and width of the period distribution. Our re-analysis of literature spectroscopic datasets for LeoII, Draco, Ursa Minor, Carina, Sextans, and Fornax imply that there is no evidence at present for varying properties of the binary stellar population in these galaxies or from deviations from that of the solar neighborhood.

From the point of view of the methodology, our analysis appears to return robust uncertainties. However, there are some improvements that could still be worth implementing. For example, in the current implementation, one loses all the information about which velocity variations correspond to which star, while it would be interesting to retain this information, because if a star shows large β for all the combinations of observations, they are potentially more informative of velocity variation due to binaries than an isolated one that could be caused by other effects. Moreover, we are studying only absolute velocity variations, however, we could also use the information about the sign of the difference between velocities; this could be especially interesting for obtaining periods of individual stars with future datasets.

Even though current datasets are still rather unconstraining, it seems likely that it will be possible to shed light on the universality of the binary fraction of MW dSphs in the near future. We forecast that seven repeated measurements of the same 750 brightest stars obtained over only two years in a classical dSph like Draco and an instrument like WEAVE in low resolution MOS mode should deliver a much improved precision on f, with 1-σ uncertainties about half of those from current spectroscopic data and an uncertainty of ±0.2 for a significance of 3σ. Beside WEAVE and 4MOST, as well as other future instruments (HRMOS, MSE, and WST...) will be able to provide better multi-epoch datasets with even lower velocity uncertainties. However, we are close to a level of precision in the velocities where the jitter should be taken into account too. In order to make good use of more precise measurements, the jitter should be well characterized, so it may be included in the models.

Acknowledgments

J. M. Arroyo acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia en Innovación (AEI-MICIN) and the European Social Fund (ESF+) under grant PRE2021-100638. J. M. Arroyo, G. Battaglia, G. Thomas acknowledge support from the Agencia Estatal de Investigación del Ministerio de Ciencia en Innovación (AEI-MICIN) and the European Regional Development Fund (ERDF) under grant number AYA2017-89076-P, the AEI under grant number CEX2019-000920-S and the AEI-MICIN under grant number PID2020-118778GB-I00/10.13039/501100011033. The authors acknowledge the referee for the constructive report which enhanced the quality of the manuscript.

References

Appendix A: Common binary fraction across MW dSphs

We will check whether we can reject the hypothesis that MW dSphs share the same value of the binary fraction. We define as our null hypothesis that all the dwarfs have a binary fraction of 0.5. Under this assumption, all the PDFs that we see in Sect. 4 would be the result of noise in the measurements. We decide to use 0.5 because all the PDFs have this value in the 1σ interval, and 4 of them in the 1σ_hdi. This is more or less what we expect to have for our sample of seven dwarf galaxies, given that this range indicates a 68.2% confidence level.

In order to test our null hipothesis, we will follow the next steps. (1) First, we randomly extracted one value of the binary fraction for each of the dSphs according to their PDFs. Then we subtracted the value that we want to test (f = 0.5) to each of the extracted values. (2) Next, we carried out a bootstrapping for the seven values generated, exchanging positions 10000 times. Finally, we computed the mean of the value for each of the 10000 cases. We iterated steps one and two 1000 times to sample the whole PDFs. The distribution of values generated after step three is a Gaussian distribution centered in 0.02 (which indicates that 0.52 is the exact mean of the distributions). The range with a 95% confidence level around the mean is [-0.23, 0.27]. This means that we can reject our null hypothesis for the values f < 0.27 and f > 0.77 with a 95% confidence, but not for any value inside that range. However, it is important to point out that this test is biased since our PDFs only allow binary fractions between 0 and 1. Therefore even if we only had flat PDFs the value of 0.5 would still be slightly favored.

This is why we also carried out a more qualitative analysis. We define the metric:

$$\begin{array}{c}\hfill D_f={\displaystyle \sum _{i=1}^7}(f-f_{\mathit{dwarf}}^i),\end{array}$$(A.1)

where f is a certain binary fraction that we want to test and $f_{\mathit{dwarf}}^i$ is the mode of the PDF for the i galaxy. We used this mode because this is the worst case scenario, as the value differs more than the median from 0.5. Therefore, D_f indicates how close our results are to a single value of f. We will test the binary fraction f = 0.52. The value of D_0.52 in our observational case is 1.53. Then, we simulated measurements of a constant binary fraction (f = 0.5) 7 times with the precision given by the 1σ_hdi uncertainties. So, we can compare if the result we obtain is what we expect if all the dwarfs have a common binary fraction. If we simulate 10000 cases we obtain a Gaussian centered in 1.24 and with a width of 0.35 for the expected values of D_0.5. This means, that in the case that all the dwarfs have a binary fraction of 0.5, we obtain the expected results given our precision in the measurements, even within 1σ. Therefore, we cannot neglect the null hypothesis with this method either.

All Tables

All Figures

	Fig. 1. Distribution of Gaia eDR3 sources (in gray) within a 50 arcmin radius around the central coordinates of the Sculptor dSph, with overlaid the 96 targets of the FLAMES-GIRAFFE multi-epoch observations (large dots). Left panel: distribution on the sky. Right panel: distribution on the color-magnitude diagram. The color indicates the number of repeated observations (see legend).
In the text

	Fig. 2. Distribution of time baseline and number of repeated observations for the 96 stars that are part of our Sculptor dataset.
In the text

	Fig. 3. Spectra of an example Sculptor member star, observed during three different nights, as indicated in the legend. The top panel shows the whole spectra. The bottom panels show a zoom-in view of the CaT lines, with the x-axis indicating the velocity shift with respect to the rest wavelength (indicated with dashed lines) of a given CaT line.
In the text

Fig. 4. Distributions of velocity differences and α. Upper panels: histograms of the difference between velocities for the different nights with respect to the reference night, as indicated in the legend, before the zero-point offset correction (left) and after the correction (right). Lower panels: same as the upper panels but for α instead of the velocity difference. Only three nights have been shown for ease of visualization, but all nights have been corrected.

In the text

Fig. 5. Probability distribution functions for the parameters defining the properties of the binary star populations. The black lines show the theoretical distributions while the colored histograms show the results of our sampling. Upper panels: period distribution for the two limits in mass-ratios, as indicated by the legend, assuming a constant minimum orbital distance amin = 0.21 AU (left) or the taking into account the RLOF when computing amin (right). The limits for each mass-ratio are plotted with vertical lines. Middle left panel: distributions adopted for the mass ratio (green: GS distribution from Duquennoy & Mayor 1991; salmon: uniform distribution). Middle right panel: distribution of inclinations. Lower left panel: eccentricity distribution; the differences between the uniform theoretical distribution and the one sampled are due to the dependence of e on the semi-minor axis and therefore, the orbital period, as explained in the main text. Lower right panel: distribution of the argument of the periastron.

In the text

Fig. 6. Cumulative distribution function of β for Sculptor. The black lines show the CDF obtained for Sculptor from the FLAMES dataset; the CDFs from 750 examples of mock datasets for a binary fraction of 0, 0.55 and 1 are shown in green, blue, and salmon respectively. As there is some overlap between the cases of f = 0.55 and f = 1, the median cumulative curve is highlighted in each case with brighter colors and a thicker line-width. Inset: Example for just one simulation, where in gray we highlight the area between the CDF of the simulated case and that for the Sculptor observations, that is what we use when computing the likelihood function.

In the text

	Fig. 7. PDFs obtained for Sculptor, using the VLT/FLAMES dataset. In salmon, we show the result assuming amin = 0.21 AU and in blue the one taking into account the RLOF, when computing amin. The median and confidence levels of 1σ and 2σ are indicated with horizontal lines.
In the text

	Fig. 8. Cumulative distribution function of β for Sculptor. In black, we show the observational distribution, in red the one produced by averaging the result of the 10 000 simulations for the most likely model using amin = 0.21 AU and in blue for the RLOF when computing amin. The panels show different ranges for the axes in order to better visualize the results.
In the text

	Fig. 9. Posterior distribution function for f obtained for Sculptor, Draco and Ursa Minor generated using the methodology by S18.
In the text

Fig. 10. PDFs for 100 mock galaxies simulated reproducing the characteristics of the FLAMES dataset of Sculptor, in terms of number of stars, observations of each star, and velocity uncertainties of each observation, for a binary fraction f = 0.6. Lower panel: PDFs of the 100 mock galaxies overlaid. Upper panel: distribution of the maxima of each PDF with a Gaussian fit overlaid. In both panels, the green and orange vertical lines indicate the median of the distribution of maxima and the range that encloses the 68% of the simulations (1σ confidence level).

In the text

Fig. 11. PDFs obtained for the FLAMES Sculptor dataset when simultaneously obtaining two parameters between the binary fraction, $\overline{logP}$ and σlog P, and keeping the third one fixed to the values in the legend. From left to right, the parameter being kept fixed is the width of the period distribution, the mean of the period distribution and the binary fraction. The probabilities have been normalized so that if we sum a whole column, we obtain 1.

In the text

Fig. 12. PDFs obtained by summing the colormaps of Fig. 11. Here, each bin correspond to the sum of a whole colormap. In the upper panel, we show the PDF for the binary fraction summing all the distribution of periods simulated. In the lower left panel, the one for $\overline{logP}[\mathrm{days}]$ summing all possible binary fractions and σlog P. In the lower right panel, the one for σlog P summing all possible binary fractions and $\overline{logP}[\mathrm{days}]$.

In the text

	Fig. 13. PDFs on f for the MW dSphs indicated in the legend from our re-analysis of literature datasets, using the new method proposed in this work. The median and 1σ and 2σ confidence levels are indicated with horizontal lines with ticks.
In the text

Fig. 14. PDFs for simulated WEAVE datasets of Draco. Left panel: PDFs for the 100 simulated WEAVE datasets of Draco (assumed binary fraction f = 0.6) obtained using the new methodology implemented in this work. In orange, we show the cases where the correct binary fraction is obtained within 1σ (69%); in brown within 2σ (21%, up to 90%) and in black those that predict it within 3σ (10%, up to a 100%). The vertical green line indicates the median of the maxima distribution. Right panel: one single PDF of the ones shown in the left panel, this is what we expect to obtain when analyzing the future observational data. Vertical lines indicate the median and intervals of confidence of 1σ, 2σ, and 3σ, as indicated in the labels.

In the text