A closer look at the binary content of NGC 1850

ABSTRACT

Studies of young clusters have shown that a large fraction of O-/early B-type stars are in binary systems, where the binary fraction increases with mass. These massive stars are present in clusters of a few Myr, but gradually disappear for older clusters. The lack of detailed studies of intermediate-age clusters has meant that almost no information is available on the multiplicity properties of stars with |$M\lt 4\, {\rm M}_{\odot }$|⁠. In this study we present the first characterization of the binary content of NGC 1850, a 100 Myr-old massive star cluster in the Large Magellanic Cloud, relying on a VLT/MUSE multi-epoch spectroscopic campaign. By sampling stars down to M = 2.5 M_⊙, we derive a close binary fraction of 24 ± 5 per cent in NGC 1850, in good agreement with the multiplicity frequency predicted for stars of this mass range. We also find a trend with stellar mass (magnitude), with higher mass (brighter) stars having higher binary fractions. We modelled the radial velocity curves of individual binaries using the joker and constrained the orbital properties of 27 systems, ∼17 per cent of all binaries with reliable radial velocities in NGC 1850. This study has brought to light a number of interesting objects, such as four binaries showing mass functions f(M) > 1.25 M_⊙. One of these, star #47, has a peculiar spectrum, explainable with the presence of two discs in the system, around the visible star and the dark companion, which is a black hole candidate. These results confirm the importance and urgency of studying the binary content of clusters of any age.

1 INTRODUCTION

Decades of observations have revealed that massive star clusters (⁠|$M\gt 10^{4} \, {\rm M}_{\odot }$|⁠) of all ages are key tools for studying the formation and evolution of stellar populations. Their detailed characterization has opened a window to previously unknown phenomena: from the existence of multiple stellar populations (e.g. Bastian & Lardo 2018), whose origin still lacks a coherent explanation, to the extended/bimodal main-sequence turn-offs (e.g. D’Antona et al. 2015; Milone et al. 2015), a common feature of young clusters originating from different distributions of stellar rotational velocities (e.g. Bastian & de Mink 2009; Kamann et al. 2020, 2023).

Star clusters are very dense environments, in which interactions and collisions between stars tend to occur much more frequently than in the field. These interactions have a massive impact on the binary populations inside the clusters (Heggie 1975). While binaries with low binding energies tend to be quickly destroyed, tightly bound binaries are subject to a variety of phenomena, including fly-by or exchange interactions. This facilitates the formation of close binary systems and exotic objects, such as blue straggler stars, low-mass X-ray binaries, millisecond pulsars, or cataclysmic variables (Bailyn 1992; Paresce, de Marchi & Ferraro 1992; Ferraro et al. 1999, 2009) as well as peculiar rapidly rotating B-type stars like Be and shell stars (Pols et al. 1991; Zorec & Briot 1997; Bodensteiner, Shenar & Sana 2020a). The multitude of interactions will cause the binary populations to strongly evolve with cluster age, resulting in evolving binary fractions, orbital parameter distributions, or pairing fractions. A clear indication of this process are the low binary fractions observed in ancient globular clusters (e.g. Milone et al. 2012), relative to those found in the Galactic field for stars of the same mass (Moe & Di Stefano 2017).

With the detection of gravitational waves from merging black holes (BHs) (Abbott et al. 2016a, b), the BH populations of star clusters have become a focus for the field. Given the high interaction rates and the tendency to form massive binaries via exchange interactions, star clusters appear as ideal environments for the dynamical formation of binary BHs (e.g. Rodriguez, Chatterjee & Rasio 2016; Di Carlo et al. 2019). In particular, merger cascades (or hierarchical mergers) appear feasible in star clusters, and offer a natural explanation for events such as GW190521 (Abbott et al. 2020), involving BHs beyond the mass limits typically expected from stellar evolution. However, the merger rates inside star clusters are still uncertain (see discussion in Mapelli et al. 2021). Open questions regarding the multiplicity properties of massive stars or the retention of BHs following supernova kicks (e.g. Atri et al. 2019) all contribute to this uncertainty.

Observationally, the direct evidence for the existence of BHs in star clusters is still limited, with a handful of objects in binaries with luminous companions having been reported thus far (Maccarone et al. 2007; Strader et al. 2012; Chomiuk et al. 2013; Miller-Jones et al. 2015; Giesers et al. 2018, 2019). In large part, this situation can be attributed to the challenges of performing spectroscopic studies in the crowded cluster environments, limiting our abilities to detect BHs (or other types of remnants) in binaries with luminous companions. However, the advent of powerful integral field spectrographs (IFSs) like MUSE (Bacon et al. 2010) is changing the landscape, as evidenced by the discovery of several BHs in the Galactic globular cluster NGC 3201 by Giesers et al. (2018, 2019).

Giesers et al. (2019) also highlighted the potential of MUSE to study the binary populations of star clusters as a whole and could constrain the period and eccentricity distributions of binaries inside NGC 3201 as well as the interplay between binarity and exotic objects such as blue stragglers and sub-subgiants. This raises the question if a similar approach can also be used to study young massive clusters, which may not only harbour a larger number of BHs compared to ancient globulars, but also show a number of peculiar features that have been linked to binarity. Examples in this respect are the unknown origin of the extended/bimodal main-sequence turn-offs (D’Antona et al. 2015; Bastian et al. 2020; Wang et al. 2022) or the rich populations of Be stars (Bodensteiner et al. 2021; Kamann et al. 2023).

Many spectroscopic studies have recently focused on characterizing the binary content of OB associations and very young open clusters, both Galactic and extragalactic. These studies have mainly sampled O- and early B-type stars (massive stars, ∼10 M_⊙ or above) and have nicely shown that the binary fraction is the highest for the most massive stars, and then it rapidly decreases with decreasing mass. This is in agreement with the observational results of Moe & Di Stefano (2017) and Sana et al. (2012). For example, a bias-corrected close binary fraction of 51 ± 4 per cent and 58 ± 11 per cent, was found among young O- and B-type stars in the 30 Doradus region of the Large Magellanic Cloud (LMC, a few Myr; Sana et al. 2013; Dunstall et al. 2015), of 52 ± 8 per cent among B-type stars in the young Galactic open cluster NGC 6231 (2–7 Myr; Banyard et al. 2022), of at least 40 per cent in the OB population of Westerlund 1 (10–20 Myr; Ritchie et al. 2022), of 34 ± 8 per cent among B-type stars in the SMC cluster NGC 330 (35–40 Myr; Bodensteiner et al. 2021) etc. Over the years, further efforts have also been devoted to try and characterize the observational properties of individual binaries (in terms of period, eccentricity, mass ratio etc.). Interesting results in this direction were recently published, for example, for 30 Doradus (Almeida et al. 2017) or Westerlund 1 (Ritchie et al. 2022). However, all binaries in these clusters are made of massive stars, of 10 M_⊙ or more. As a matter of fact, little or no information is currently available for binaries with intermediate-mass (2–5 M_⊙) stellar components.

To fill this gap at intermediate masses, we initiated an observational campaign focused on analysing the binary content of NGC 1850, a 100 Myr-old massive cluster in the LMC, using multi-epoch observations acquired with VLT/MUSE in two distinct observing runs. The target is ideal for two main reasons: (1) it is one of the most massive clusters in the LMC (M = 0.52 x 10⁵ M_⊙; Song et al. 2021), ensuring good statistics; (2) it belongs to a cluster age range, and corresponding star masses, yet to be explored in terms of their binary population. Furthermore, it is young enough so that little dynamical evolution has affected their binary properties, but yet old enough so that the population of exotic objects (e.g. BHs) should be already formed and potentially present in the cluster. This data set has already allowed us to trace the rotational distribution and relative binary fraction of stars in the blue and red main-sequences (a common feature in clusters of this age; Kamann et al. 2021, 2023) as well as to detect NGC1850 BH1, a binary system proposed to host a BH candidate (Saracino et al. 2022). This idea was later questioned (El-Badry & Burdge 2022; Stevance, Parsons & Eldridge 2022) but the nature of the unseen component is not yet established (Saracino et al. 2023).

In this work we explore the binary content of NGC 1850 by sampling stars down to three magnitudes below the main-sequence turnoff of the cluster. Thanks to a statistically significant sample of binaries (of the order of hundreds), we can investigate the distribution of their main orbital properties (for example their periods). Observational studies as the one presented here are urgent and timely as they will provide important constraints for modelling massive star clusters, with the aim of obtaining more solid predictions in terms of their stellar and non-stellar population content. The urgency of such studies can be understood, for example, in the context of the controversy raised about the nature of NGC1850 BH1. In fact, Stevance et al. (2022) used the binary evolution code BPASS to demonstrate that a binary like the one suggested by Saracino et al. (2022) should not exist in nature, being unsupported by stellar evolution. However, to reach this conclusion, the authors have used binary properties that are representative of the field, instead of massive clusters such as NGC 1850. In these systems, dynamical exchanges and interactions occur at a very high rate, with a strong impact on the evolution of the entire binary population and on the properties of individual binaries. How exactly they are influenced is poorly known.

The paper is structured as follows: in Section 2 we briefly describe the observations and data reduction, we present the methods used to measure radial velocities and elaborate the approach used to identify binary systems. In Section 3 we present the sample and discuss the intrinsic binary fraction of NGC 1850, within the context of the MUSE field. Section 4 is dedicated to the analysis of the radial velocity curves of binaries. The orbital properties and distributions of constrained binaries are presented in Section 5. In Section 6 we describe the properties of individual binary systems, and elaborate on those with peculiar characteristics. After an extensive discussion of the results, we draw our conclusions in Section 7.

2 OBSERVATIONS AND DATA REDUCTION

A detailed description of the observations available for NGC 1850 can be found in previous papers using the same data set (see Kamann et al. 2021, 2023 and Saracino et al. 2022). Here we limit ourselves to provide a brief description of the data, mainly focusing on their importance for detecting radial velocity variables (i.e. binary stars), which is the main purpose of this work. The data set of NGC 1850 consists of multi-epoch observations (Programme IDs: 0102.D-0268 and 106.216T.001, PI: N. Bastian) acquired with MUSE (Bacon et al. 2010), the IFS mounted at the very large telescope (VLT), in Chile. We used the wide field mode (WFM) configuration assisted by adaptive optics, to improve the spatial resolution of the field and facilitate the extraction of single star spectra in such a high density environment. Since MUSE covers a field of view (FOV) of 1 × 1 arcmin at a spatial sampling of 0.2 arcsec in WFM, we used two pointings to sample the central regions of the cluster [the effective radius is indeed r_e = 20.5 ± 1.4 arcsec (4.97 ± 0.35 pc), from Correnti et al. 2017] and collect a statistically high number of stars. The two pointings overlap each other slightly, in order to gain deeper observations of that region (see fig. 1 in Kamann et al. 2023).

The observations cover a temporal baseline of more than 2 yr (specifically 754.1 d), with a time sampling between individual epochs ranging from 1 h to several months. We secured 16 observations per pointing and up to 32 independent observations for stars within the overlapping region. This strategy allows us to have an almost complete coverage of the inner regions of NGC 1850 both spatially and temporally, opening up the possibility of studying the kinematics of the stars within this massive star cluster. More quantitatively, the number of stars observed for a given number of epochs is presented in Fig. 1. The clear peak we see at 16 epochs corresponds to the default scenario where one spectrum per star/epoch is obtained. A small peak at 32 epochs is also observed and refers to the ‘best observed’ stars (i.e. stars in the overlapping region). A strategy that adopts non-uniformly sampled epochs allows us to be sensitive to different types of binaries, and to discriminate between radial velocity signals of short and long period binaries. This is very important as it helps us to avoid confusion caused by time aliasing and to be confident in determining the orbital periods of the binaries in the sample. The reduction of the MUSE datacubes has been performed using the ESO MUSE pipeline (Weilbacher et al. 2020) in the ESO Reflex framework. By default, the pipeline subtracts the sky lines and the sky continuum separately from the data. We did only perform the subtraction of the sky lines because the sky continuum in a crowded field like NGC 1850 ultimately contains a contribution of starlight that we did not want to remove. As a result, the final datacubes in our case still contain the telluric continuum emission and we model its contribution a posteriori, since it is spatially flat across the FOV and does not affect the spectra extraction.

Individual stellar spectra have been extracted from the MUSE FOV using the latest version of pampelmuse (Kamann, Wisotzki & Roth 2013), which relies on a point spread function (PSF) fitting technique. We used a high-quality photometric catalogue of NGC 1850 obtained from archival observations from the Hubble Space Telescope (HST) as a reference to determine both a MUSE PSF model and spatial coordinates of the resolved sources as a function of wavelength. This information is crucial to resolve and extract single sources even in the most crowded regions. We refer the reader to Kamann et al. (2023) for more details on this, including links to the publicly available photometric catalogue and the combined MUSE spectra.

Interestingly, NGC 1850 is one of the LMC clusters in the WAGGS (WiFeS Atlas of Galactic Globular cluster Spectra) project (Usher et al. 2017). The WAGGS project exploits the wide-field spectrograph (WiFeS) on the Australian National University (ANU) 2.3-m telescope (Dopita et al. 2007, 2010) which provides high resolution (R ≈ 6800) spatially resolved spectroscopy, with a wide wavelength range (3270–9050 Å) covered by four gratings, of which two are observed simultaneously. The main aim of this project was to deduce the properties of star clusters through integrated spectroscopic studies however we exploited the capabilities of this instrument to obtain additional observations of NGC 1850 in order to measure the rotation of stars at its main-sequence turn-off as well as to look for radial velocity variability (PI: G. Da Costa). The observations acquired in 2018 September were specifically designed to have a perfect overlap to the MUSE FOV, in order to add further epochs to the radial velocity curves of bright targets (F438W < 18), as well as to have higher resolution spectra compared to MUSE. We extracted and analysed the WiFeS data cubes of NGC 1850 with pampelmuse and the procedure described below applies consistently to both data sets to derive single-epoch radial velocities for each star in the FOV. We note here that the WiFes radial velocities will not be used to derive the binary fraction in NGC 1850. They will only be exploited as additional epochs to constrain the orbital properties of some binaries, where possible.

2.1 Radial velocities with spexxy

As described in Saracino et al. (2022) and Kamann et al. (2021, 2023), the extracted spectra were analysed with spexxy (Husser et al. 2016), a python package which performs full-spectrum fits against a library of template spectra. The main aim of this process is to determine radial velocities and stellar parameters (e.g. effective temperature, metallicity) for every single source in the field. Rotational velocity measurements are not discussed in this paper as they have been already the subject of a dedicated study (Kamann et al. 2023). Given the presence in the field of a large number of hot stars with T_eff >10 000 K, we used the synthetic templates from the ferre library (Allende Prieto et al. 2018), which includes spectra for stars with up to T_eff ∼30 000 K. The initial guesses required by spexxy, like surface gravity log g, and effective temperature T_eff were taken from the comparison of the HST photometry to rotating and non-rotating MIST stellar models (Choi et al. 2016; Gossage et al. 2019), assuming age of 100 Myr and metallicity of [Fe/H] = −0.2, appropriate for NGC 1850. To avoid spurious measurements, we masked out any spectral ranges that could potentially be contaminated by the intense nebulosity that permeates the NGC 1850 field, including the cores of the strong Balmer lines (H α and H β).

We analysed the spexxy results a posteriori and considered as reliable all those that meet the following criteria: (1) The fit was marked as successful. (2) The derived signal-to-noise (S/N) of the input spectrum was >5. (3) The velocity determined by cross-correlating the spectrum with its best fit deviated by less than 3σ from the spexxy result. (4) The Mag Accuracy parameter returned by pampelmuse, which measures the agreement between the m_F814W magnitude derived from an extracted spectrum and the one in the HST photometry, was above 0.6. Spectra that did not meet these criteria were discarded, thus excluded from the final spexxy catalogue of MUSE radial velocities, which contains 3283 stars.

2.2 Radial velocities with the cross-correlation

In analysing the radial velocities measured with spexxy, we paid particular attention to two specific classes of objects: (1) stars for which spexxy failed to provide a good full-spectrum fit, due to the presence of very peculiar spectral features, impossible to reproduce with any stellar template (see for example the case of star #47 and star #2413, discussed in Section 6). As the analysis typically failed for such stars, they are not included in the spexxy catalogue. (2) Stars that are heavily contaminated by nebular emission. For these stars, a full-spectrum fit could underestimate any radial velocity variations. A solution in such cases consists of analysing only the red part of the spectra‘, i.e. the Paschen series, where nebula emission is negligible. For a detailed discussion of this issue, we refer to the case of star #224 or NGC1850 BH1 in Saracino et al. (2023).

To avoid unreliable or biased radial velocity measurements and also to recover the radial velocities of stars not analysed by spexxy, we decided to apply an alternative method to derive the relative radial velocities of the same NGC 1850 sample, which relies on cross-correlation (CC, hereafter) of the observations with a template spectrum created from the data themselves (Zucker & Mazeh 1994; Shenar et al. 2017; Dsilva et al. 2020). We performed the CC by exploiting the wavelength range between 7800 Å and 9300 Å. For each star, the spectrum with the highest S/N was used as an initial template. Afterwards, each remaining spectrum was cross-correlated against the template. Then, a new template was created by shifting each spectrum by the inverse of its measured radial velocity and creating a S/N-weighted average of the shifted spectra. To avoid that individual low-S/N spectra impact the process, only spectra with |$\gt 40~{{\ \rm per\ cent}}$| of the S/N of the initial template were included in the average. We repeated this process until the measured radial velocities converged. We note that this process results in relative rather than absolute velocities. When necessary, we converted to absolute velocities by analysing the final combined template with spexxy. We applied the same S/N cut and the same magnitude accuracy as for spexxy, to create the final CC catalogue of MUSE radial velocities, which contains 2662 stars.

Since the two main goals of the study are: (1) to estimate the binary fraction in NGC 1850 and (2) to constrain the orbital properties of as many reliable binaries in the cluster as possible, we decided to select a ‘high-quality’ sample of stars, for which we obtained completely consistent results with spexxy and cc. All the results shown hereafter will be based on this final, clean hybrid catalogue. A very efficient way to identify ‘high-quality’ stars is to compare the binary probability p derived from the two methods, instead of using individual radial velocities for each star. The binary probability, in fact, will be immediately correlated to how similar (or different) the results of spexxy and cc are for that specific star. This selection will be discussed in the next section, where the approach to estimate the binary probability for each star in NGC 1850 is presented.

2.3 Binary probability

The main goal of this work is to investigate the presence of binary systems in the MUSE field of NGC 1850 and to characterize their properties, both in terms of overall distributions and of individual systems. For this reason, radial velocities determined from single-epoch spectra are the main information we are interested in. When looking for any signs of variability in the data, three different groups of stars can be identified: (1) stars in binary systems, where the variation in radial velocity is due to the motion of the star around a companion; (2) stars that are intrinsically variables, also called pulsators (RR Lyrae or Cepheids are common examples); (3) stars that show small variations due to the finite accuracy of our measurements but are not true variables. To perform this study we are interested in stars belonging to the first group, therefore an accurate knowledge of the radial velocity uncertainties is essential to discard all possible interlopers.

To calibrate the velocity uncertainties returned by both the spexxy code and the CC method, we applied the approach presented in Kamann et al. (2020). For each measured radial velocity in our sample we followed three steps: (1) first, we identified a subsample of 100 stars with similar T_eff, log g, and S/N values as the target star in the extraction under consideration (comparison sample). (2) Then, we measured the velocity differences of the stars in the comparison sample relative to every remaining epoch and divided the differences by the squared sums of the velocity uncertainties. Obvious binary stars were excluded from this comparison iteratively via k-sigma clipping. (3) Finally, we determined a correction factor for the formal uncertainty of the radial velocity measurement under investigation by measuring the standard deviation of the epoch-to-epoch distribution of the normalized velocity differences in the comparison sample. The fact that an independent correction factor can be obtained by comparing to each of the remaining epochs further allowed us to estimate the stability of the velocity calibration.

Velocities of evolved stars [F438W<18.0 and (F336W-F438W)>0.0] in NGC 1850 are measured to an accuracy between 1.0 and 2.0 |$\, {\rm km\, s^{-1}}$| during each epoch. On the main sequence, the typical uncertainties (per epoch) are larger, ranging from 5.0–6.0 |$\, {\rm km\, s^{-1}}$| at F438W ∼16. to 10.0 |$\, {\rm km\, s^{-1}}$| at F438W ∼18. to 25.0–30.0 |$\, {\rm km\, s^{-1}}$| for the faintest stars in our sample. Please note that the raw (i.e. before calibration) velocity uncertainties derived from the CC routine are on average larger (by a few per cent, which is intrinsic to the method) than those obtained with spexxy, so the calibrated uncertainties will also be larger in the latter case.

After carrying out an in-depth evaluation of the radial velocity uncertainties, the next step is to determine the probability that each star in the sample is a radial velocity variable. To do so we used the method developed by Giesers et al. (2019) for NGC 3201. For a given star observed for a number of epochs, we computed the reduced χ² value for the set of radial velocity measurements and corresponding uncertainties under the null hypothesis of a single star (no radial velocity variations). The degrees of freedom are derived from the number of epochs available per star. This distribution over all observed stars is compared with what would be expected under the assumption that no radial velocity variables are available in the sample. More practically, for any given reduced χ² value, this method computes the ratio between the number of stars observed above this value and the number of stars expected above this value in the absence of radial velocity variable stars. This comparison allows us to assign a probability p of being radial velocity variable to each star. High values of p correspond to stars that are most likely part of a binary system, while low values of p are generally associated to single stars. As a general assumption, we define as likely binary stars all those with p > 0.5. This criterion will also be used later in the paper for creating the star catalogue that will be analysed with the Bayesian software the joker.

We applied this method to both the spexxy and the CC results independently, in order to identify (1) which stars are considered binaries from both methods; and (2) which stars have ambiguous results, i.e. they are considered binaries by spexxy but not by CC and vice versa. We realized that a non-negligible number of stars which were considered binaries (p > 0.5) according to the spexxy results, turned out to be single stars (p << 0.5) when processing the CC results. This is mainly caused by template mismatches (e.g. Be stars with emission lines), with the net effect of artificially inflating the fraction of binary stars in the inner regions of NGC 1850. In order to get a clean sample of stars, with reliable radial velocities, we proceeded as follows:

• First, we matched the stars in the two catalogues, finding 2541 stars in common. To clean up this sample from stars with poorly measured radial velocities, we required the variability probabilities in the two approaches to be within 0.2 of each other. For example, if a star has p = 0.7 according to spexxy, the same star analysed with CC needs to have 0.5 < p < 0.9 in order to be included in the final, cleaned, catalogue. A sample of 1441 stars met the criterion, corresponding to ∼57 per cent of the initial sample.

• Secondly, we looked at stars analysed by only one of the two approaches (spexxy or CC). In particular, 742 stars from the spexxy catalogue and 121 stars from CC. By looking at the magnitude distribution of these stars, we realized that stars in spexxy-only sample are all very faint (F814W>20). These stars have low S/N spectra, which is the reason why they are not included in the CC catalogue. On the other hand, the CC-only stars cover the entire magnitude range, with an overabundance of stars at intermediate magnitudes (F814W∼18–19). The latter are stars with peculiar spectra that cannot be fitted with standard stellar templates.

• We created the final, clean catalogue, by adding to the 1441 stars with consistent results between spexxy and CC, the 121 stars from CC-only, obtaining a final sample of 1562 stars. Although this sample does exclude a fraction of stars (mostly faint), it can still be considered a good representation of stars in the MUSE field of NGC 1850, i.e. those with reliable radial velocities over different methods. We will use spexxy radial velocities for all stars measured by the two methods, because of the smaller uncertainties, while we will adopt radial velocities provided by the CC approach in all other cases. The CC technique is in fact the most reliable for stars belonging to the two classes mentioned above or to other peculiar categories. We have verified, however, that the results presented hereafter (in terms of binary fraction and orbital properties of constrained binaries) do not change if only the CC radial velocities are used.

In Fig. 2 we present the (F438W, F336W-F438W) colour–magnitude diagram (CMD) of NGC 1850, where each coloured point is a star for which reliable MUSE radial velocities are available from the hybrid, clean catalogue. Stars most likely members of NGC 1850B are not shown in the figure. These stars were all located preferentially in the extension of the NGC 1850’s main-sequence towards brighter magnitudes, overlapping the region where blue straggler stars are located. The binary probability distribution is shown as a colourbar, where the darker colours refer to single stars, while the lighter ones identify the binaries. Stars that are certainly binaries (with p > 0.95) are also highlighted in the figure, with large red open circles.

To enable the reader to more easily follow the remainder of the analysis, we provide below a brief summary of the stellar samples we will be using and the criteria used to generate them.

1. To estimate the spectroscopic binary fraction of NGC 1850 we use:

   • only cluster members (1385 stars after removing 177 possible interlopers);

   • stars with a minimum of five valid radial velocity measurements and p > 0.5 are considered likely binaries;

   • a magnitude cut (F814W<19.5) to ensure completeness of the photometric sample;

   This selection results in 109 binary candidates, out of a sample of 1033 stars remaining after the cuts.

2. To derive the properties of individual binary systems with the joker we use:

   • all stars in the MUSE FOV with reliable radial velocities (no membership cuts, 1562 stars);

   • stars with a minimum of five valid radial velocity measurements and p > 0.5 are considered likely binaries;

   • no magnitude cuts;

   This selection results in 143 binary candidates, out of a sample of 1562 stars remaining after the cuts.

3 THE BINARY FRACTION OF NGC 1850

Based on the statistical approach used to estimate the binary probability for each star in the sample, we can now compute the observed spectroscopic binary fraction for the inner regions of NGC 1850 in the MUSE FOV. The MUSE field is characterized by stars belonging to three different environments: (1) the massive, 100 Myr-old cluster NGC 1850; (2) the 5-Myr old, low mass, cluster NGC 1850B, responsible for the strong nebulosity observed in the field (see for example the right panel of fig. 1 in Kamann et al. 2023). As this work focuses on the study of the binary fraction and the binary population of the main cluster NGC 1850, it is important to remove the contribution of the interlopers, in order to have a sample of cluster members as clean as possible. NGC 1850 stars and stars belonging to the field share similar radial velocity values, thus a distinction based on radial velocities remains ambiguous. In fact, Song et al. (2021) has measured for NGC 1850 a systemic velocity of v_sys = 248.9 |$\, {\rm km\, s^{-1}}$| and a velocity dispersion of σ = 2.5 |$\, {\rm km\, s^{-1}}$|⁠, while for the LMC field a systemic velocity of v_field = 257.4 |$\, {\rm km\, s^{-1}}$| with a dispersion of v_field = 23.6 |$\, {\rm km\, s^{-1}}$|⁠. The latter is also consistent with the uncertainties provided by other literature works studying the LMC field, for example with the GAIA data (Vasiliev 2018).

Likely members of NGC 1850B and the LMC field can actually be identified on the basis of photometric and/or spatial properties. In fact, from previous works in the literature, as well as from a visual inspection of the cluster CMD in specific filters, the contamination from field stars appears not severe. In addition, LMC stars are generally much older (∼6 Gyr) than NGC 1850 stars, so the LMC main-sequence overlaps with NGC 1850 only at rather faint magnitudes, starting at about F438W∼21.5, corresponding to the faintest stars in the MUSE sample, which are not included in the final clean catalogue. To be more specific, from Fig. 2, field stars are predominantly located at F438W>18. and colours (F336W-F438W)>0.2–0.25. The photometric argument does not really work for NGC 1850B, since its main-sequence partially overlaps with that of the main cluster, NGC 1850. However, its contribution in stars can still be inferred by looking at the (RA-Dec) parameter space: these two clusters are very close but still distinguishable in the sky. All stars within 10 arcsec of a visually determined centre (α = 05^h08^m39.3^s, δ = −68^○45′45|${_{.}^{\prime\prime}}$|5) can be considered likely members of NGC 1850B.

To determine the observed spectroscopic binary fraction of the main cluster NGC 1850 we applied the photometric and spatial cuts mentioned above to remove the most probable members of NGC 1850B and the LMC field from the sample of 1562 stars with reliable radial velocities. This is safe to do, since this result will not be significantly affected by a few misclassified stars. However, we anticipate here that we will run the joker on the entire data set, as we want to be as conservative as possible when looking for the orbital solutions of each individual binary in the cluster.

The catalogue of members of NGC 1850 consists of 1385 stars with reliable radial velocity measurements. All stars are detected at least in five MUSE epochs. Fig. 3 shows the F814W magnitude distribution of all stars in this sample, as well as the corresponding binary fraction, by adopting bins of 1 mag each. The observed binary fractions are obtained from the comparison between the number of stars with p > 0.5 in each magnitude bin and the total number of stars with reliable radial velocity measurements in that bin. The statistical uncertainties of the binary fractions are calculated using the prescriptions of Giesers et al. (2019): by the quadratic propagation of the uncertainty determined by bootstrapping (random sampling with replacement) the sample and the difference of the fraction for p > 0.45 and p > 0.55 divided by 2 as a proxy for the discriminability uncertainty between binary and single stars. In the upper panel of Fig. 3 two important aspects need to be pointed out: (1) the bins corresponding to the brightest magnitudes contain a limited number of stars, leading to larger uncertainties in the binary fractions due to limited number statistics; (2) for magnitudes F814W>19.5, the number of stars per bin rapidly decreases. This is where our spectroscopic sample starts to become incomplete. The binary fraction increases for the faintest stars: binaries are brighter, so they will be overrepresented in this incomplete sample. If we focus on the binary fraction for stars at magnitudes 15.5 < F814W < 19.5, we then observe a decreasing trend, with the binary fraction varying from 30 per cent to ∼7 per cent. An alternative way to visualize this difference is by plotting the binary probability distribution for stars in two different mass regimes: stars with high masses (M > 4 M_⊙,¹ corresponding to F814W < 18.) in the left panel of Fig. 4 and stars with intermediate-masses (2.5 < M < 4 M_⊙, corresponding to 18 < F814W < 19.65 mag) in the right panel of the same figure. As already suggested by the observed trend in magnitude of Fig. 3, the binary fraction in the high-mass end is 17.7 ± 2.8 per cent, more than double the one in the intermediate-mass regime, 7.3 ± 1.8  per cent.

Finally, to determine the overall spectroscopic binary fraction of NGC 1850 in the MUSE FOV, we take into account all stars with F814W < 19.5, to assure that our spectroscopic sample is complete. We derive an observed binary fraction of 10.6 ± 1.8 per cent, meaning that in a sample of 1033 stars, we find 109 stars with p > 0.5, i.e. likely binary systems. We also investigated how the observed spectroscopic binary fraction varies as a function of the distance from the centre of NGC 1850, since our data set samples nearly three effective radii of the cluster (∼60 arcsec). Although the uncertainties are significant, we observe that the binary fraction decreases with increasing distance from the centre, approximately from 12 per cent for the innermost bin to 9 per cent for the outermost bin. This is consistent with the mass segregation phenomenon already observed in other clusters (e.g. Giesers et al. 2019). We note here that variable stars (e.g. Cepheids or RR Lyrae) are partially inflating these numbers, as there is no way to distinguish them at this level from true binaries. The sample will be cleaned a posteriori of variable stars (mainly Cepheids) thanks to cross-matching with archival catalogues (e.g. OGLE-IV). A discussion on this aspect can be found in Section 6.3, but the low number of detected variables in the MUSE FOV does not have any impact on the overall binary fraction of NGC 1850.

3.1 The binary fraction among Be stars

In very young clusters stars are used to be classified as not-evolved, i.e. still on the main-sequence, and evolved, i.e. giant stars that have already evolved off the main-sequence. In NGC 1850, the latter group is significantly underpopulated compared to the former, so providing the binary fractions for the two groups separately is not very informative. However, unlike other clusters of similar age, NGC 1850 is known to host a significant number of Be star candidates (see Correnti et al. 2017; Milone et al. 2018), i.e. rapidly rotating B-type stars with a circumstellar decretion disc that gives rise to strong Balmer-line emissions (Rivinius, Carciofi & Martayan 2013). Exploiting the same data set we are using here, Kamann et al. (2023) have studied this population of stars in detail from a spectroscopic perspective, by selecting the sample of Be stars with reliable radial velocity measurements but also discriminating between ‘classical Be stars’ and shell stars, i.e. Be stars observed through their discs (viewed at very high inclination angles), thanks to a detailed inspection of their spectra. In this work we add another piece of information to the Be star picture, by studying the observed spectroscopic binary fraction of Be star candidates in NGC 1850 and compare it with the fractions found in other clusters in the literature.

Fig. 5 shows the CMD of NGC 1850 in a specific combination of HST filters (F814W, F336W-F814W) that allow us to easily identify Be stars in the cluster and clearly separate shell stars. The spectroscopic sample of Be and shell star candidates from Kamann et al. (2023) are highlighted by large red squares (202 stars) and large green diamonds (47 stars), respectively, while the colour code for all stars is the same as in Fig. 2. The inset in Fig. 5 shows the binary probability distribution of all Be star candidates in red and of the subsample of shell stars in green. We find 20 Be stars with p > 0.5, measuring an observed spectroscopic binary fraction of 10.0 ± 2.8 per cent, while we do not find any shell star showing signs of variability (f_{obs, shell} = 0 per cent). In order to understand how these fractions compare with the fraction of binaries for other classes of objects in NGC 1850, we did not apply a completeness correction factor,² but we performed a relative comparison to stars in the main-sequence of the cluster, at the same magnitude range of Be stars. The amount of completeness is expected to be the same for stars in the same magnitude range. For this subsample of main-sequence stars we derived an observed spectroscopic binary fraction of 11.0 ± 1.9 per cent. Within the uncertainties, Be and main-sequence stars have similar binary fractions in NGC 1850. This is interesting as it slightly differs from what was found in other clusters. In NGC 330, for example, Bodensteiner et al. (2021) found that the observed binary fraction of Be stars was nearly half that of main-sequence stars (15.9 per cent versus 7.5 per cent). It is worth noticing, however, that if we only focus on obvious binary systems (p > 0.8 or higher), we do find f_{obs, MS} exceeding f_{obs, Be}. On the other end, the very low binary fraction they found among Be stars ∼0.3 mag redder than the main-sequence (2 per cent, see their fig. 7), i.e. what we call shell stars, is in excellent agreement with what we find here. It is not entirely clear yet why shell stars are mostly observed as single stars, however there might be a few effects at play. Bodensteiner et al. (2021) for example suggested two possibilities: (i) the close binary fraction of classical Be stars is intrinsically low; or (ii) Be binaries have long orbital periods, for which the detection probability drops. Only a more detailed investigation of this class of objects can help shed light on this aspect.

3.2 Correction for observational biases

In order to compare the overall binary fraction of NGC 1850 with other studies in the literature, but also to investigate the reliability of the magnitude/mass trends observed (i.e. whether they are caused by our ability of detecting binaries at different magnitude levels, or they are linked to an intrinsic difference in the multiplicity properties of NGC 1850’s stars at different magnitudes/masses), a proper understanding of the different biases affecting our measurements is needed. The main limitations of our observational campaign are the following: (1) Not all orbital periods have the same probability of being detected. In particular, due to the temporal sampling of our survey, we are more sensitive to tight binaries than to long period ones; (2) Many binaries have orbital velocities below our detection threshold, so they may have been classified as single stars when they are not. (3) The method we use is actually biased against detecting binaries of (almost) equally bright stars. In such cases, both binary components contribute (almost) equally to the observed spectrum. In the combined spectrum, the radial velocities of the two components mostly cancel out, making it extremely challenging to detect such binary systems based on low-resolution spectroscopy. In a recent work with MUSE data on the 35–40 Myr-old cluster NGC 330, Bodensteiner et al. (2021) found that their probability to detect binaries with a mass ratio q > 0.8 drops to 10 per cent or less, due to the phenomenon just described.

To test our ability in detecting binary systems given our observational setup (e.g. in terms of time coverage, number of epochs, radial velocity uncertainties), we performed two different sets of simulations. The main difference between the two is on the distribution of the parent populations adopted for three orbital parameters, such as period, mass ratio (q = M₂/M₁), and eccentricity. We adopted the same distributions as in previous works of this kind (e.g. Bodensteiner et al. 2021 and references therein) to allow a consistent comparison of the results. For the first set we assumed a log-uniform period distribution ranging from 0.15 to 3.5 (i.e. P from 1.4 to 3160 d), an eccentricity distribution between 0 and 0.9 proportional to |$\rm \sqrt{e}$|⁠, with a circularization correction for periods < 2 d (see discussion in Section 5.2 and equation 1) and a flat mass ratio distribution between 0 and 1 (see the finding in Shenar et al. 2022 for O-type stars). In the second set of simulations we instead adopted power-law distributions for all three parameters: |$\rm (logP)^\pi$| with π = −0.25 ± 0.25 for the period (Bodensteiner et al. 2021); |$\rm q^k$| with k = −0.2 ± 0.6 for the mass ratio and |$\rm e^\eta$| with η = −0.4 ± 0.2 for the eccentricity (Sana et al. 2012). The exponents were drawn from normal distributions with central values and 1σ dispersions taken as described above.

We generated 1385 binary systems (the number is chosen to be equal to the number of stars in the cleaned sample of NGC 1850), by randomly assigning orbital parameters (e.g. period, eccentricity, mass ratio etc.) to each of them from the parent populations presented above. We further adopted a random orientation of the orbital plane in 3D space and a random reference time t₀. Of all binary systems thus created, we discarded those with unphysical solutions based on two criteria: (1) binary hardness, i.e. how much the system is bound. The equation adopted in the simulations is equation (5) from Ivanova et al. (2005), which determines whether the binary system survives the cluster environment based on its binding energy. All systems with binary hardness < 1 are treated as single stars. (2) Roche Lobe overflow, i.e. when the radius of the most massive star in a binary exceeds its Roche limit. When this happens, the system is in the common envelope phase and can no longer be treated as a binary in this context. In order to take into account the effect on the radial velocities of SB2 binaries, or binaries with almost (equally) bright components, we also implemented the damping factor formula derived by Giesers et al. (2019) (see their equation 1), as they found that the theoretical expected radial velocities are linearly damped with the flux ratio of the two components. The closer the mass ratio q is to 1, the closer the combined radial velocities are to the systemic velocity, completely erasing the signal and making these types of systems very difficult to detect.

Once all radial velocities were generated for all binaries surviving the two criteria (∼ 94 per cent of the initial sample), we adopted the same approach described in Section 2.3 to estimate the binary probability of each system in the sample, then its binary fraction, comparing the systems with p > 0.5 to the total number of surviving binaries. Since these are simulations, we know exactly how many binaries and single stars we have in the sample and this allows us to estimate the detection probability, i.e. how good we are in detecting binaries given the properties of our observations. For each of the two configurations we have iterated the simulations 10 000 times, to allow a robust calculation of the detection probability, but also to have a good idea of the associated uncertainties.

We obtained a detection probability of 43.4 ± 4.7 per cent from the analysis of the first set of simulations, while |$\rm 45.1^{\rm -15.3}_{\rm +9.2}$| per cent for the second set. These values were obtained by computing the 50th percentile of the distributions generated by the 10 000 simulations. When using these detection probabilities to correct the observed spectroscopic binary fractions, we obtain a true binary fraction for stars in NGC 1850 of 24.4 ± 5.0 per cent and |$\rm 23.5^{\rm -15.4}_{\rm +9.4}$| per cent for the two sets of simulations, respectively. The uncertainties associated with the detection probability and, as a consequence, to the true binary fraction of NGC 1850 are made of two terms: (1) a systematic uncertainty due to the number of simulations performed; (2) a statistic uncertainty due to the binary population sample size. We obtained the former by calculating the 16th and the 84th percentiles (± 1σ) of the above distributions.³ To estimate the latter, we instead ran two more sets of simulations (of 10 000 realizations each), this time using the true binary fractions as input parameters to the simulations. The number of binaries in these simulations corresponds to the intrinsic number of binaries in the NGC 1850 sample, thus the width of the error distributions (always ± 1σ) reflects the impact of the sample size.

Although the two sets of simulations have very different parent population distributions, the intrinsic, bias-corrected, binary fractions we obtain for NGC 1850 are remarkably similar to each other, considering their uncertainties. The latter are significant for the simulations with power-law distributions for period, mass ratio, and eccentricity but this is not surprising given the large variety of exponents that could be randomly assigned to these parameters. Indeed, they could switch from distributions in period that favour tight binaries (π ∼ -0.5), which are easy to detect with our observational setup, to those that favour much longer orbital periods (π ∼ 0), very difficult to detect. This has the net effect of increasing the spread in the detection probabilities, therefore increasing the uncertainties. The shape of the period distribution is by far the dominant factor that affects the resulting detection probabilities, as already noted by Bodensteiner et al. (2021).

To visually inspect the results of our simulations, in Fig. 6 we show how the detection probability varies as a function of the orbital period and mass ratio of the binaries. We consider only the first set of simulations (with flat |$\rm logP$| and mass ratio distributions), and use 10 realizations. We include all binaries in the sample, without applying any cut. We note here that by increasing the number of simulations used, the results would not change, only smoother curves would be obtained. As can be seen, the detection probability decreases with increasing orbital period as our observational setup makes them more difficult to detect. It ranges from 0.9 for periods of a few days to 0.2 for periods of thousands of days. For the mass ratio, a very low detection probability (∼ 0.1) is observed for binaries with low q values, reaching a peak around 0.5 for q = 0.5–0.7, then decreasing for q > 0.8. The observed behaviour is consistent with what was found by Bodensteiner et al. (2021) in their fig. 2 (orange lines) using a similar set of simulations. In this work, however, we observe that the drop in the detection probability is significantly less drastic than in the simulations by Bodensteiner et al. (2021). It is reasonable to link this discrepancy to differences in the observational setup between the NGC 1850 data set and the one available for NGC 330.

When we do the same exercise as in Fig. 6 but we now consider two different mass regimes, i.e. high-masses: M >4 M_⊙ and intermediate-masses: 2.5 < M < 4 M_⊙, we obtain the results presented respectively as a blue and a red curve in the same figure. Overall, the detection probability is higher in the high-mass regime (blue), than in the intermediate-mass regime (red). This is expected since more massive stars are also brighter in the CMD; they have high S/N spectra, and their radial velocities are measured with lower uncertainties, allowing binaries with significantly lower amplitudes to be detected than in the other regime. To investigate whether the trends observed in Fig. 3 as a function of magnitude and in Fig. 4 as a function of mass are still present after correcting for the observational biases just mentioned, we calculated the true binary fractions to compare with the observed spectroscopic binary fractions. The bias-corrected values are shown as red diamonds in Fig. 3 and are labelled in Fig. 4. Although the differences between the different regimes have become less important, trends in magnitude/mass can still be observed, which may actually reflect an intrinsic difference in the multiplicity properties of these stars.

3.3 Comparison to other clusters

The close binary fraction in the 100 Myr-old cluster NGC 1850 can be directly compared to the fraction measured in other clusters in the literature. As mentioned in Section 1, most of the clusters analysed so far are actually very young (of the order of a few or few tens of Myr) and still contain a large fraction of O- and early B-type stars with masses up to ∼10 M_⊙ or more, see for example 30 Doradus, NGC 6231, Westerlund 1 etc. NGC 1850 is the oldest massive cluster for which a detailed spectroscopic analysis of the binary content has been performed. Due to the massive stars involved, the close binary fraction in these clusters is much higher, over 50 per cent in 30 Doradus (Sana et al. 2013; Dunstall et al. 2015) and NGC 6231 (Banyard et al. 2022) and over 40 per cent in Westerlund 1 (Ritchie et al. 2022) than what we derived in this work for NGC 1850, but if we also include NGC 1850 in the sample we clearly observe a trend with mass. A much closer comparison in terms of true binary fraction can actually be made with the 35–40 Myr-old cluster NGC 330 in the low metallicity environment of the SMC, Z∼1/5 Z_⊙ (Bodensteiner et al. 2021). This cluster has a binary fraction of 34 ± 8 per cent but the observations sample stars down to 4 M_⊙. Although the overall binary fraction in NGC 1850 is still lower than what the authors find in NGC 330, if we compare the binary fraction of NGC 1850 in the high-mass regime (left panel of Fig. 4), where only stars with masses of 4 M_⊙ and above are included (thus setting the same low-mass threshold in the two clusters), the two close binary fractions come into agreement (34 ± 8 per cent versus 29.5 ± 3.7 per cent). This result confirms that the binary fraction in a cluster is indeed mass-driven, with other cluster properties such as metallicity or age potentially influencing these values as well, but much less strongly.

4 DETERMINATION OF BINARY PROPERTIES

The statistical approach applied in the previous sections has provided us with a sample of binary candidates, i.e. all stars with binary probability p > 0.5. For each of these objects we have an observed radial velocity curve which, although only sparsely sampled, contains much information about the orbital properties of the binary itself. Several methods can be used to infer these properties, and one of the most used in recent years is the generalized Lomb–Scargle periodogram (GLS; Zechmeister & Kürster 2009). It is very effective for finding the orbital period of binary systems but usually works best when uniformly sampled data sets are available. Considering the nature of our data, this approach would only be successfully applied to a small subset of stars, particularly those with a populated radial velocity curve, where the orbital period is fairly well sampled. In all other cases, unfortunately, the GLS method would not be ideal.

For that reason, we exploit a slightly different approach here, which was only recently developed and has already proven to be very successful with similar sets of data. Called the joker, this tool is a custom Monte Carlo sampler for sparse or noisy radial velocity measurements of two-body systems, and can produce posterior samples for orbital parameters even when the likelihood function is poorly behaved (Price-Whelan et al. 2017, 2020). We applied this software to our NGC 1850 data set with the goal of constraining the properties of individual systems, whose sampling and number of epochs are good enough to unequivocally determine their orbits.

4.1 THE JOKER

Literature studies have demonstrated that at least five epochs are needed to constrain the orbit of a binary system (Price-Whelan et al. 2017). Therefore, in order to run the joker, we focused only on stars with at least five reliable MUSE spectra, corresponding to five radial velocity measurements, and, as mentioned previously, with a binary probability p > 0.5. We used the entire data set, without any a priori selection based on the probability of the stars to belong to the main cluster NGC 1850.

To make the joker work properly on our data set, we also made two main assumptions: (1) all binaries are made of two stars, so no triple or multiple systems are included in our modelling. (2) in our binary systems, one of the two stars always dominates in light with respect to the other. This means that we have modelled all binaries as being SB1, even though we know this assumption is not entirely correct. SB2 stars, i.e. stars in which the two components are equally bright and contribute almost equally to the system’s light, in fact, do exist in clusters like NGC 1850 but can be identified only after a proper inspection of the stellar spectra and in general with higher spectral resolution.

The idea behind the joker is to create a huge library of possible orbits, based on a set of input parameters, such as period, eccentricity, semi-amplitude velocity K etc. When the radial velocity curve of a star is provided, the software scans the defined parameter space to find the family of orbits that best match the observed radial velocity curve. If this procedure is successful and the orbit is well constrained by the data, only one set of possible orbits (all with very similar orbital periods) is shown, providing a unique solution. Otherwise, depending on the behaviour of the data, the joker can provide a few possible solutions or a set of solutions randomly distributed over the entire allowed period range. The latter means that the orbit of that specific binary is poorly constrained or not at all constrained, respectively.

In particular we generated 2²⁹ = 536 870 912 prior samples, uniformly distributed within a period range between 1 d and 1000 d. We did not include any orbital solution shorter than 1 d as this would have introduced an artificial bias toward unphysical, very short period binaries. We checked whether some of the stars in our sample showed a preference for periods below the lower limit we imposed but could not find any within the constrained sample (see Section 6 for details). For the upper limit we assumed a period slightly larger than our temporal baseline, in order to give enough room to the software to find the appropriate solution. Binaries with orbital periods longer than 1000 d may indeed be present in the cluster but would probably have such a low semi-amplitude velocity K that it would be impossible for us to detect them as binaries.

As for the eccentricity and the velocity semi-amplitude K, we have followed the prescriptions given in table 1 of Price-Whelan et al. (2020), i.e. a beta distribution for e and a normal distribution centred in 0 |$\, {\rm km\, s^{-1}}$| and with σ = 30 |$\, {\rm km\, s^{-1}}$| for K. However, we have carried out several tests by modifying the distribution for some of these parameters (e.g. log-normal distribution versus uniform distribution in P or standard versus custom normal distribution in K) and we get very consistent results for the constrained binaries, which gives us confidence that the results do not depend on the type of distributions adopted.

We assumed a normal distribution for the systemic velocity of the binaries, centred on |$v_0=250.0\, {\rm km\, s^{-1}}$|⁠, which is close to the systemic velocity of NGC 1850 and with |$\sigma =20.0\, {\rm km\, s^{-1}}$|⁠, which is consistent with the velocity dispersion of LMC field stars (Vasiliev 2018). The systemic velocity of the binaries in NGC 1850 is expected to follow the distribution with the dispersion of the cluster (⁠|$\sigma =5.0\, {\rm km\, s^{-1}}$|⁠), however, this choice is conservative because it ensures that we cover all cluster binaries but do not bias our prior against velocity of field binaries, thus avoiding the misinterpretation of some results. On the other hand, it gives us the possibility to identify a posteriori those binaries that are likely part of the cluster, by assigning them a membership probability through photometric and/or spatial location in the MUSE FOV.

We requested 512 posterior samples (solutions) per star and, as mentioned above, we analysed only stars with a minimum of five epochs and a binary probability p > 0.5 or higher. As soon as one star gets fewer than 256 posterior samples (half of the total), two different strategies are applied, depending on how constrained the solution is. (1) If the solution is nearly unimodal, we use a dedicated Monte Carlo Markov chain (MCMC) run as described in Price-Whelan et al. (2017, 2020) to get 512 new samples. (2) If the solution is more widespread, we perform an iterative procedure, generating new samples until reaching the minimum number of 256. The results of the entire process, both in terms of the binary population properties and of individual peculiar systems are presented in the following sections.

5 THE BINARY POPULATION OF NGC 1850

We ran the joker on a sample of 143 binary systems fulfilling the criteria described in Section 4. This analysis can result in three very different outcomes: (1) a unimodal distribution, i.e. the solution is unique and the orbital parameters are well determined; (2) a bimodal distribution, i.e. two sets of solutions are possible, thus the orbit is close to be constrained; (3) a multimodal or continuous distribution, i.e. the orbital parameters are poorly constrained or completely unconstrained. Among those, we only focused on objects with unimodal or bimodal posterior distributions.

In order to discriminate between unimodal, bimodal, or unconstrained solutions, we applied the method described in Giesers et al. (2019). For each star we computed the standard deviation of the posterior periods P returned by the joker on a logarithmic scale. Stars with σ_ln P < 0.5 are identified as unimodal and represent the most constrained systems in our sample. We find that 20 out of 143 binaries have unimodal solutions. An example of a binary with a unimodal solution is presented in Appendix  A (Fig. A1) for star #260, where the phase-folded radial velocity curve⁴ is also presented, to show how good the agreement between the data and the best-fitting model is. Bimodal posterior samples were instead determined by using a k-means clustering algorithm with k = 2 from the scikit-learn python package (Pedregosa et al. 2011) to separate two sets in the period posterior distribution, where each of the two sets has to fulfill the criterion above. The set which includes the largest number of samples is used as the preferred solution in these particular cases. To determine the best-fitting orbital period per binary, we take the 50th percentile of its distribution as the median and the 16th and 84th percentiles as 1σ respectively below and above the median. Starting from this period range (P_min – P_max), we then extract the corresponding distribution in all other parameters (e.g. eccentricity, semi-amplitude velocity, systemic velocity etc.). The best-fitting values for each of these parameters are finally determined as the 16th, 50th, and 84th percentiles of the newly created distribution. We find bimodal solutions for 1 out of 143 binaries in our sample. The binary with a bimodal orbital solution is star #303 and is presented in Fig. A2.

The results of the joker in terms of star position (RA and Dec coordinates), m_F438W magnitude and best-fitting values of the fitted parameters (period, eccentricity, semi-amplitude K, systemic velocity v_sys, etc.) are listed in Appendix  A in Table A1 for all binaries with constrained solutions. The last column helps identifying which binary belongs to which group (unimodal versus bimodal). To verify the accuracy of the constrained solutions derived by the joker, we analysed the data using two more approaches. In particular we applied the GLS code mentioned above, as well as the ultranest software (Buchner 2021) and we obtained for all stars listed in Table A1 consistent results within their uncertainties.

All the remaining 116 binaries in the MUSE sample have multimodal or continuous distributions (their properties are not constrained), hence they cannot be used to describe the properties of the binary population in NGC 1850.

5.1 Binaries constrained with MUSE+WiFeS data

In the sample of binaries analysed with the joker, 6 out of 143 deserve particular attention. These binaries are all characterized by orbital periods which are preferentially longer than 100 d. Although such a period is well within the temporal baseline of our MUSE observations, the strategy we have used to design the MUSE campaign is such that the data are much more sensitive to binaries with orbital periods of the order of a few days or tens of days with respect to hundreds. For this reason, the solutions provided by the joker for these systems are all multimodal when considering only MUSE epochs. Adding new epochs to the sample is really the only way to extend the temporal baseline of our data set and more importantly its sampling, in order to be able to constrain even the longest orbits of these systems. By exploiting the complementary WiFeS programme of NGC 1850 (PI: G. Da Costa), we were able to extract the spectra of some bright stars in the WiFeS FOV. More specifically, the brightness of our targets (F438W < 18) combined with the perfect overlap between the WiFeS data and the MUSE FOV made it possible to add further epochs to the radial velocity curves of 64 binary stars. In this way, much better constrained orbital properties were derived for each of the six binaries mentioned above. The remaining 58 binaries can be easily divided into two categories: (1) binaries whose orbital properties were already well constrained using the MUSE data; and (2) binaries that remained unconstrained even after the WiFeS data were added.

5.2 Orbital properties of constrained binaries

Almost 17 per cent of the total number of binaries detected in the MUSE FOV of NGC 1850 have well constrained orbital properties, in terms of orbital period, eccentricity, semi-amplitude velocity K etc. These properties are listed, along with their uncertainties, in Table A1 and include all 27 binaries with unimodal and bimodal solutions, both with MUSE and MUSE+WiFeS. The number of binaries with constrained orbits was actually larger, but the final sample was created by removing photometric variable stars, e.g. pulsators, which we identified by cross-matching our data set with the OGLE-IV catalogue (a discussion on this topic is presented in the next sections).

We tried to assign each of the 27 constrained binaries a probability of being a member of NGC 1850, NGC 1850B, or the LMC field, and found that 25 out of 27 are very likely members of the massive cluster NGC 1850. Indeed, based on the HST photometry, there are no stars with colour (F336W-F438W)>0.5 and magnitude m_F438W>18, so we can safely conclude that none of them are part of the LMC field. However, there are two stars (#32 and #51) that are spatially within a radius of less than 10 arcsec from the centre of NGC 1850B. This leads us to believe that these stars are probably members of NGC 1850B, although this is not certain, as many members of NGC 1850 may still occupy that region.

The eccentricity–period distribution of the 27 constrained binaries in NGC 1850 is presented in a linear-log plot in the main panel of Fig. 7. The black dots identify binaries with unimodal posterior samples, while the red dot is the binary with bimodal posterior samples, as obtained by using the MUSE epochs. The orange hexagons instead indicate binaries with unimodal solutions obtained thanks to the addition of further epochs from the WiFeS data set. Eccentricity and period distributions are also shown as grey histograms in the vertical and horizontal panels, to better visualize the results.

We find binaries in a rather large period range, from ∼1.3 to ∼500 d, but their period distribution does not uniformly cover the entire range. It shows, instead, multiple peaks, corresponding to the period intervals in which our observations, due to the temporal sampling of the MUSE campaign, were more sensitive. The peak at ∼ 400–500 d is almost entirely produced by the binaries constrained thanks to the WiFeS data. Overall, the ±1σ uncertainty in period is rather small, except for a few binaries which have significantly larger uncertainties (though we note that the x-axis is logarithmic). The latter, despite being classified as unimodal or bimodal based on the criterion presented at the beginning of the section, show posterior samples clustered in more than one group at periods very close to each other. This affects the width of the distribution, thus the uncertainties.

The eccentricity distribution of the constrained binaries, on the other hand, covers only a narrow range of values from 0 to 0.4, with a broad distribution and perhaps the hint of a peak around 0.1–0.15. This suggests that the binaries in our sample have relatively low eccentricities. However, eccentricity is the least constrained orbital property of all. Consequently, with the exception of a very few outliers, any value of e between 0 and 0.4 is allowed for most binaries. Finally, there are no highly eccentric binaries (e > 0.5) in the sample. It is not surprising to see that the eccentricity distribution is biased towards low eccentricities, because for such orbits fewer radial velocity measurements are necessary to get a unique solution.

In Fig. 7 we also show as a cyan dashed line a maximum eccentricity e_max power law derived from a Maxwellian thermal eccentricity distribution:

for a given period P (Moe & Di Stefano 2017). This represents the binary components having Roche lobe fill factors ≤70 per cent at periastron. They predict that all binaries with P < 2 d should have circular orbits due to tidal forces. It is interesting to note that, with the exception of one star (star #1828, see Table A1), all the others behave accordingly. The period distribution presented in Fig. 7 is now plotted against the peak-to-peak radial velocity distribution (ΔV_r = 2K) of the 27 binaries, in a log–log plot in Fig. 8. The colours are the same as in Fig. 7. Here we observe a large range of amplitude values among the binaries, specifically from tens of |$\, {\rm km\, s^{-1}}$| to hundreds of |$\, {\rm km\, s^{-1}}$|⁠, with one system even exceeding 300 |$\, {\rm km\, s^{-1}}$|⁠. The properties of this peculiar binary system NGC1850 BH1 have been already discussed (see El-Badry & Burdge 2022; Saracino et al. 2022; Stevance et al. 2022; Saracino et al. 2023).

According to Clavel et al. (2021), this plot can be used to identify the region where stars with massive companions are possibly located. A main-sequence turn-off star in NGC 1850 is of ∼5 M_⊙, which is the maximum mass a star in such a cluster can assume, according to its absolute age (100 Myr). We used this assumption to derive the position, in the (P–ΔV_r) parameter space, of a binary made up of two components: (1) a primary star as massive as a 5 M_⊙ star (the star we observe), and (2) a secondary component with exactly the same mass, so that the mass ratio of the system (q = M₂/M₁) is equal to 1. The locus of points with such a property is identified by the cyan dashed line in Fig. 8.

All binary systems having q = 1 but masses for the individual components lower than 5 M_⊙ define lines parallel to the left of the cyan dashed line (the lower the mass, the more they move to the left). Conversely, the cyan shaded region to the right of the cyan line is populated by binary systems in which the secondary (or invisible) component is significantly more massive than the primary (or observed) star (for M₁ ≤ 5M_⊙, which is an appropriate assumption for stars in NGC 1850). Binaries with massive invisible companions are very interesting as they are the best candidates for hosting dark compact objects, such as NSs and BHs. In our sample of 27 constrained binaries, there are four that fall in this region, with two being very close to the dashed line. As labelled in the plot, these binaries have HST ID #23, #47, #157, and #224. We label the latter as BH1 instead, to adopt the same nomenclature as the papers that have discussed its properties (Saracino et al. 2022, 2023). We will discuss the properties of star #23, star #47, and star #157 in detail in a dedicated section (6). The binary mass function, f, is an important orbital property we can derive directly once we know the period, the eccentricity, and the semi-amplitude velocity K of a binary system, thanks to the formula:

where G is Newton’s gravitational constant. Without knowing a priori the mass of the visible source, this quantity can instruct us on the properties (e.g. mass) of the source that we do not detect. In fact, using Kepler’s third law, the binary mass function f(M) can also be written in the form:

where M₁ and M₂ are the masses of the primary visible star and the secondary unseen component, respectively, and i the inclination angle of the binary with respect to the line of sight. In massive star clusters, where we know the maximum mass a star can assume, the binary mass function f(M) is a very informative parameter since, without making any assumptions on the mass of the observed star and/or the inclination of the system, it is able to predict which of the two components should be more massive. In NGC 1850, for example, where the maximum mass for a star is ∼ 5 M_⊙, if we consider the most conservative scenario, i.e. the binary system is nearly edge-on (i ≈ 90°), we obtain that the secondary (unseen) component is surely more massive than the primary (visible) component if f(M) >1.25 M_⊙. The lighter the observed star, the greater the mass difference, for a given value of f.

Fig. 9 shows the histogram of the distribution of the binary mass function for our sample of 27 constrained binaries. As can be seen, most of the stars exhibit a low binary mass function (f(M)≪1 M_⊙). This is consistent with what we expect from the visual inspection and analysis of their spectra. In fact, they all appear to be SB1 binaries, i.e. binaries in which the flux contribution of the observed star is significantly larger than that of the companion, which is not detected. In most cases, this is due to the fact that the invisible source is so faint (therefore of low mass) that its contribution is negligible. Looking at the distribution, however, we note that there are four systems which have a mass function exceeding 1.25 M_⊙, with f(M) values between 1.25 M_⊙ and ∼10 M_⊙. With significantly massive companions, these systems are the same ones we identified in Fig. 8 and are very interesting systems as they are the most likely candidates to host a compact dark object, such as a neutron star (NS) or a BH. We sound a note of caution here that there may be several reasons why the contribution of a second bright star might not be detected in stellar spectra. For example, if it rotates so fast that all the spectral lines become very broad (e.g. LB-1; Shenar et al. 2020 and HR 6819; Bodensteiner et al. 2020b) or if it has an accretion or decretion disc that partially blocks its light. Further observations are needed to characterize the nature of the unseen sources in the binaries with f(M) >1.25 M_⊙, however a compilation of all the information we have is presented later in the text. The main results of this work can be well summarized by modifying the (F336W-F438W, F438W) CMD of NGC 1850 in Fig. 2 to include all the information we obtained from the analysis performed with the joker, as well as the comparison with other photometric catalogues in the literature. Fig. 10 shows all stars observed in the MUSE FOV for multiple epochs and for which we have estimated a probability to be in binary systems. Dark colours indicate a high probability that they are in binaries, while light orange colours suggest that they are single stars. The 27 binaries with constrained orbital properties are shown as large red dots and are found mostly in the upper part of the main-sequence, being bright. The binary probability is larger than 90 per cent (p > 0.9) for all constrained binaries, except for the one with bimodal solution. Fig. 10 also highlights different types of stars/binaries to emphasize some classes of objects that we have identified in our sample. Eclipsing binaries are shown as large green squares (including star #51, that we have classified as a member of NGC 1850B) while binaries with massive companions (f(M) >1.25 M_⊙) are indicated by large cyan diamonds. The large yellow star denoting the BH candidate BH1, also belongs to the latter category (its properties are discussed in Saracino et al. 2022, 2023). The orange hexagon identifies a binary system that could harbour a stripped subdwarf O (sdO) star candidate. Finally, large pink triangles highlight the variable stars (mainly Cepheids) found in NGC 1850, all of which are located in the evolved portion of the CMD. Each of these object classes will be discussed in the next section, where more details will be provided.

6 INTERESTING OBJECTS IN NGC 1850

The detailed analysis of the binary content of NGC 1850 has shown the presence of a few very interesting binary systems, of very different natures. Here we discuss those cases separately.

6.1 Binaries with massive companions

The orbital properties derived by the joker suggest that a sample of 4 out of 27 binaries have a binary mass function f(M) > 1.25 M_⊙. They are the binaries with HST ID #23, #47, #157, and #224. The latter has been already extensively discussed in Saracino et al. (2022, 2023) and will not be mentioned here. Instead we focus on the others later in the text.

The first interesting property we note is when a system has a value of f(M) > 1.25 M_⊙. In the context of NGC 1850, this tells us that the component of the binary we do not see is more massive than the star we observe. In the first place, this finding only means that the mass ratio q of the binary is above 1, if it is defined as q = M_unseen/M_visible, but it does not provide any information about the nature of the object we do not observe.

Indeed, a source could be massive but not observed because: (1) it is rotating so fast that the blurring of the absorption lines hamper the lines from the two components to be observed; (2) there is an accretion disc which blocks a significant fraction of its light, making it difficult to detect; (3) it is a dark object, such as a NS or a BH, so it emits no observable light.

We have looked more into the binary system with the highest binary mass function, f(M)=10.53 M_⊙ or 14.46 M_⊙ and we report the details below.

6.1.1 Star #47

For star #47 we have a total of 25 radial velocity measurements, spanning a range of more than three years. In particular, nine observations come from the WiFeS instrument (WAGGS survey), while the remaining 16 epochs come from MUSE observations. This star has been observed only in one of the two MUSE pointings of NGC 1850, i.e. the central pointing, meaning that we were able to extract a high S/N spectrum for every image secured.

According to the criterion we have adopted so far in the paper, star #47 has a unimodal solution, however, when looking closer at the posterior samples, the joker finds for the system two possible orbital solutions, with periods P = 508.4 d and P = 984.3 d. These solutions, shown in Fig. 11, are very well defined when WiFeS and MUSE observations are analysed together, however they are already well recognizable when only the MUSE data are used. The phase-folded radial velocity curve of star #47 using the first and the second best-fitting solution are presented in Figs 12 and 13, respectively, where also residuals are shown. The two solutions are equally good, as demonstrated by the residuals as well as the reduced χ² value reported. It is clear that star #47 belongs to the class of long-period binaries. There is no evidence to suggest that any short-period solutions can be missed. Furthermore, MUSE observations were designed to be extremely sensitive to binaries with orbital periods of the order of days, so given the available measurements, if the system was a short period binary, we would definitely have detected it.

As presented in Table A1, both orbital solutions derived by the joker imply a binary mass function f(M) well above 1.25 M_⊙. In particular, f(M) = 10.53 M_⊙ for the P = 508.4 d solution and f(M) = 14.46 M_⊙ for the P = 984.3 d solution. This means that the unseen companion in this system is considerably more massive than the star we observe and this is reported in Fig. 14, where the two curves (black-dashed and red-dotted), computed based on the measured mass functions, show how massive the unseen component of the system would be as a function of the assumed mass of the visible component. On the x-axis, a range 0–10 M_⊙ is presented, indicating all the possible masses the observed star can assume (0 M_⊙ on one extreme, 10 M_⊙ on the other, if the star belongs to the 5 Myr-old cluster NGC 1850B). This plot seems to suggest that the system containing star #47 needs to host a massive BH candidate, in order to explain such a trend. While the available data are not sufficient to fully understand the nature of the system, we have analysed the spectra of this peculiar system, also applying the disentangling technique (see Saracino et al. 2023 and references therein for details on the methodology) and we summarize below all the information acquired:

• The spectra of star #47 show prominent emission in H α and H β, as well as a double-peaked emission in the Paschen series. The shape of the Balmer lines varies as a function of the orbital phase, while the shape of the Paschen series do not (see the insets in Fig. 15).

• The He i line at ∼4920 Å, and the Paschen lines are moving in phase, while the Balmer lines (H α and H β) result in close-to-constant radial velocities.

• The disentangling of the He i line and the Paschen lines yields no notable features for the presence of a secondary in the system, showing a close-to-flat spectrum. On the other hand, the disentangling of the Balmer lines identifies two components: a clear single-peak emission component, superimposed on a double-peaked emission component. The contribution in flux of the single-peak emission component is not negligible and we estimated it to contribute ≈30–40 per cent on the entire flux in H α and H β.

• We measured the rotational velocity of the observed star from the He i line at ∼4920 Å in the MUSE spectrum and we found a value close to 300 |$\, {\rm km\, s^{-1}}$|⁠. The observed star is a rapid rotator, and its rotational velocity is close to the break-up velocity observed for other stars in NGC 1850 (the detailed analysis is presented in Kamann et al. (2023) and their fig. 11).

By combining all the information, it seems clear that the spectrum of the visible star is rather similar to that of a Be-type star (also corroborated by the position of star #47 in the region of the CMD where blue stragglers are located, see Fig. 10), surrounded by a decretion disc. Another disc also appears to be present in the system, around the massive and unseen object. It is responsible for the single-peaked emission in the Balmer lines. The fact that this emission is observed at rest is in line with the expectations, because the centre-of-mass of a system made up of two components in which one is significantly more massive than the other coincides almost perfectly with the position of the massive source itself. Much information has been recovered from the MUSE spectra, but the nature of the unseen component is still unclear. One of the most plausible interpretations, however, could be that star #47 is part of a binary or a tertiary system. If it is a binary, a Be-type star with a decretion disc is orbiting with a long period around a massive object (e.g. a BH candidate, given its mass), which has its own accretion disc. In a tertiary configuration instead, star #47 would orbit, with a short period, around a donor (low-mass and faint) star from which it has gained most of the mass, and the whole binary orbits around a massive object, such as a BH candidate, on a long orbit. The accretion disc may have been generated by accretion of material onto the BH e.g. by strong winds coming from star #47.

6.1.2 Star #23 and star #157

Star #47 is not the only peculiar source observed in the inner parts of NGC 1850. In fact, there are two other sources, star #23 and star #157, which have similar properties to star #47. For example they have comparable luminosity in F438W (see Fig. 10), belong to the class of long-period binaries and are characterized by a binary mass function f(M) >1.25 M_⊙. Unlike star #47, however, the values we are dealing with here are much less extreme and that is why we treat them together.

Star #23 is located in the brightest part of the main-sequence of NGC 1850, and has 58 radial velocity measurements in total, the first 26 from the WAGGS spectra, while the remaining 32 from MUSE. This star falls in the overlapping region between the two MUSE pointings, so it has as many epochs as possible (see Fig. 1). Considered well constrained (i.e. it has a unimodal solution) based on the criterion in Giesers et al. (2019), this star has an orbital period P = 109.1 d according to the joker. Fig. 16 shows two main clusters of solutions for the star, but the one around P∼130 d is significantly less populated than the other. That is the reason why we do not treat the two clusters separately, but take them both into account when calculating the median and the uncertainty associated with its orbital parameters (see Table A1). In the bottom-left panel of the figure, we also show the phase-folded radial velocity curve of the binary, which fits well the behaviour of the observed velocities. Star #23 has an estimated f(M) = 1.26 M_⊙.

On the other hand, star #157 is rather off the main sequence, located at intermediate colours between the unevolved and the evolved sequences. We constructed the radial velocity curve using four observations from WAGGS and 16 observations from MUSE. The best-fitting orbital solution derived by the joker for this system has a period of 491.0 d, as can be seen in Fig. 17, and its phase-folded radial velocity curve also well reproduces the observations (see the bottom-right panel of the same figure). For star #157 we obtained a binary mass function of 1.45 M_⊙. Their spectra do not show striking features like that of star #47 (with the exception of the non-standard Paschen lines in the spectra of star #23), so we have not performed further investigations. The only information we have so far is that their invisible companions must be quite massive. However, the latter could be either dark remnants such as NSs and BHs, or stellar components which due to some properties of the system are not easy to detect.

6.2 Star #2413–O-type subdwarf candidate

Star #2413 is located in a particular region of the CMD in Fig. 10. It has a magnitude of F438W = 19.58 but it is significantly bluer than the main sequence of NGC 1850 (by ∼1 dex in F336W-F438W). Furthermore, it has more than 95 per cent probability of being a binary system. The star is relatively faint, so only five radial velocity measurements from MUSE survived the quality cuts we imposed to obtain the final catalogue of radial velocities in the cluster. For the same reason, no WAGGS spectra are available. The number of radial velocity measurements therefore is not sufficient for the joker to constrain the orbital properties of the system.

The information we have for star #2413 are limited, but its spectra are very interesting. No metal lines can be detected, but the He ii 4859 Å/H β blend, the He ii 6560 Å/H α blend, several other pure He ii lines and a weak but clearly detectable line of He i at 7065 Å. An example is shown in black in Fig. 18.

Interestingly, this spectrum appears strikingly similar to those presented in Götberg et al. (2018), which are the typical spectra of B- and O- type (sdB-O) subdwarf stars (i.e. very hot stars that have been stripped of most of their material due to mass transfer and are now substantially left with a He core). These stars go through a short-lived phase (Irrgang et al. 2020), in which they are quite bright, and end up in this magnitude range once they are relaxed and have reached the stripped phase. To investigate whether this might be a stripped sdO-type star, we adopted very hot stellar templates (Reindl et al. 2016) from which we estimated effective temperature T_eff, surface gravity log(g), and He abundance log(He/H) of star #2413 via full-spectrum fitting. The stellar parameters obtained from the comparison with the best-fitting stellar template (shown in red in Fig. 18) are T_eff = 50 653 ± 791 K, log(g) = 5.20 ± 0.04, and log(He/H) = −0.58 ± 0.04 (number fraction). When comparing these numbers with theoretical tracks of stripped hot subdwarf stars from Götberg et al. (2018), the location of our source in NGC 1850 matches the predictions quite well (shown in Fig. 19). This is a strong indication that star #2413 is indeed a stripped sdO star candidate orbiting a companion. We note that most of the known sdO stars are located at higher T_eff and surface gravity (Culpan et al. 2022). The T_eff and log(g) of star #2413 seem a bit too low for a post-extreme horizontal branch star (Dorman et al. 1993) and also a post-asymptotic giant branch nature seems unlikely since at T_eff-= 50 653 K, the log(g) would be required to be lower (Renedo et al. 2010). Interestingly, the Götberg et al. (2018) tracks actually predict that – for this mass – the companion should dominate the optical spectrum. In conclusion, more data are needed to characterize the evolutionary history of this star.

6.3 Pulsating stars

As mentioned above, the aim of this work is to detect radial velocity variations due to the motion of a star around a companion, rather than an intrinsic variability physically linked to the star itself (i.e. pulsations). These two classes of objects are in principle indistinguishable, as they manifest themselves in a similar manner, however there are ways to discriminate between them.

There are compilations of variable stars (e.g. RR Lyrae or Cepheids) detected via photometric monitoring with dedicated instruments for the Milky Way, as well as for the LMC and SMC. The optical gravitational lensing experiment (OGLE; Udalski et al. 1992), for example, is one of the surveys that has identified and characterized the largest number of variable systems in the last decades, even in star clusters.

To clean up our sample of these objects, we compared our final catalogue of binaries with the OGLE-IV release. We found a list of four sources in common between the two catalogues, including three Cepheids (star #96 or CEP-1082; star #306 or CEP-1078, and star #506 or CEP-1080) and one anomalous Cepheid (star #1916 or ACEP-024). Of those, three out of four have orbital properties constrained by the joker. They are not included in the 27 systems with constrained orbits and are not in Table A1 as they are not binary systems.

The period derived from the joker for these stars is in excellent agreement with the pulsational period derived from the OGLE-IV light curves, confirming the accuracy and reliability of our radial velocity measurements and results. The fourth pulsator in the sample (star #96) has a constrained solution according to the joker (P≈ 5.5 d), but the orbital period does not match with the one provided by OGLE-IV (P = 7.86 d). We have used both the spexxy and CC radial velocities for this star and the solutions we got were always the same, hence we are not able to explain where the discrepancy comes from.

6.4 Eclipsing binaries

By comparing OGLE-IV with our catalogue of binaries in NGC 1850, we identified five other sources, this time true radial velocity variables, which are classified as eclipsing binaries (ECL), i.e. binaries with an inclination such that as one star passes in front of the other, it creates a dip in the light curve of the system. ECL light curves are all very different from each other and their shape really depends on the configuration of the system, as well as the mass ratio of the components involved. Only 3 out of 5 have unimodal solutions from the joker (star #51 or ECL-08397; star #191 or ECL-08469 and star #224 or ECL-29851). The latter, also called NGC1850 BH1, was the focus of two dedicated works (El-Badry & Burdge 2022; Saracino et al. 2022) which demonstrated that the source is not an eclipsing but an ellipsoidal binary. In all three cases, the orbital period measured by the joker almost exactly matches that listed by OGLE-IV, providing independent confirmation that they are indeed binaries. Given its spatial position in the MUSE FOV, there is a high probability that star #51 is not member of the main cluster NGC 1850, but of the much younger cluster NGC 1850B. An example of an eclipsing binary is shown in Fig. A3 of Appendix  A. The phase-folded radial velocity curve of #191 as derived from MUSE (using P = 6.46 d) is presented in the left panel, while the phase-folded light curve by adopting the same orbital period is in the right one.

Eclipsing binaries have huge potential. Indeed, when both radial velocity curves and light curves are available for the same binary system, the degeneracy between the inclination of the system and the mass ratio of the two components can be broken and these properties can be conveniently inferred. By adopting i = 90° (given the observed eclipses), and visible star masses inferred from the comparison to best-fitting isochrones (M₁ ∼ 14.79 M_⊙ for star #51, since it belongs to NGC 1850B and M₁ ∼ 4.98 M_⊙ for star #191, member of NGC 1850), we used equation (3) to estimate a value of the mass ratio q for the two eclipsing binaries in our sample. In particular, we derived q = 0.14 for star #51 and q = 0.45 for star #191. These values are an approximation, since they are based on assumptions that may not be entirely true but are only used here to make a visual comparison to the CMD in Fig. 10. The position of both stars in the CMD seems consistent with the q values we have derived. In star #51, for example, the observed star dominates in light, hence its position in the CMD looks basically indistinguishable with that of a single star of the same mass (in agreement with isochrones). In star #191, the contribution of the unseen companion is more significant, but still places the binary on the turn-off of NGC 1850.

7 DISCUSSION AND CONCLUSIONS

In this paper we presented the first characterization of the binary content of NGC 1850, a 100 Myr-old and massive star cluster in the LMC. This work took advantage of a multi-epoch spectroscopic campaign performed with the IFS MUSE mounted on the VLT. Thanks to 16 epochs of observations spanning a temporal baseline of more than two years (754.1 d), we were able to identify a sample of 143 possible binaries in the MUSE FOV. This allowed us to derive an observed spectroscopic binary fraction of 10.6 ± 1.8 per cent for all stars with F814W<19.5 in NGC 1850. This fraction was corrected for incompleteness (reflecting the limitations of the observational setup of our survey) thanks to the comparison with two sets of simulations, which differed for the distributions adopted for period (uniform in log P versus power law), mass ratio (uniform versus power law) and eccentricity (power law but with two different exponents). We obtained a detection probability of 43.4 ± 4.7 per cent from the analysis of the first set of simulations, while |$\rm 45.1^{\rm -15.3}_{\rm +9.2}$| per cent for the second set. We estimated a true binary fraction of 24.4 ± 5 per cent and |$\rm 23.5^{\rm -15.4}_{\rm +9.4}$| per cent, respectively, which is smaller than what derived for other clusters in the literature (e.g. over 50 per cent in 30 Doradus, Sana et al. 2013; Dunstall et al. 2015 and NGC 6231, Banyard et al. 2022, over 40 per cent in Westerlund 1, Ritchie et al. 2022 and about 34 per cent in NGC 330, Bodensteiner et al. 2021) where much more massive stars were sampled. NGC 1850 nicely fits into this picture, proposed by Moe & Di Stefano (2017), that more massive stars tend to show higher binary fractions.

We also investigated whether the intrinsic binary fraction of NGC 1850 changed as a function of mass and after correcting for our efficiency in detecting binaries at different magnitudes, we still observed a trend, with a decrease in the intrinsic binary fraction from 48 per cent to 17 per cent, for magnitudes in the range 15.5 < F814W < 19.5. The same behaviour was also observed when splitting the NGC 1850 sample in two mass regimes. In fact, stars with intermediate-masses (2.5 < M < 4 M_⊙) have an intrinsic binary fraction of |$f_{\rm SB, \rm corr}$| = 18.1 ± 4.2 per cent, while stars with high-masses (>4 M_⊙) have |$f_{\rm SB, \rm corr}$| = 29.5 ± 3.7 per cent. When comparing the high-mass value with that found by Bodensteiner et al. (2021) in NGC 330 (34 ± 8 per cent), where they sample stars down to the same mass level (M = 4 M_⊙) we get similar results, confirming that the trend in mass observed in other clusters is also manifesting in NGC 1850.

By exploiting the high flexibility of the joker, a software designed to infer binary properties from sparse radial velocity curves, we were able to analyse our clean data set of binaries and to directly constrain the orbital properties (period, eccentricity, semi-amplitude velocity K etc.) for a sample of 27 binary systems in the MUSE FOV, all of which are likely cluster members, except for two which seem to belong to NGC 1850B. By tracing the distribution of these parameters we found that the sensitivity of our survey is not the same for all orbital periods but significantly depends on the sampling of our epochs, with the intervals between 10 and 30 d and 200 and 400 d the least sensitive. The peak at ∼ 400–500 d, i.e. an excess of points in this range, is almost entirely produced by the binaries constrained thanks to additional epochs from the WiFeS data. Of all the Keplerian parameters we have derived, the eccentricity is the most uncertain, although a preference for low eccentric orbits appears to be present.

Despite the limited number of epochs available, the MUSE survey has allowed the characterization (in terms of individual orbital properties) of approximately 17 per cent of all binaries included in the final clean data set of NGC 1850 stars with reliable radial velocity measurements. Among the sample of binary systems we identified several interesting objects such as two eclipsing binaries, a subdwarf O-type star candidate and four objects showing a binary mass function f(M) >1.25 M_⊙, i.e. NGC 1850 systems in which the invisible component is more massive than the observed star. One of these systems is NGC 1850 BH1, an ellipsoidal binary with a peak-to-peak radial velocity variation exceeding 300 |$\, {\rm km\, s^{-1}}$| which has been extensively discussed in Saracino et al. (2022, 2023). Of the three remaining systems, star #47 is really peculiar. It shows double-peaked emission both in the Balmer lines (H α and H β) and in the Paschen series and its radial velocity curve is consistent with a rather long orbital period (P = 508.4 d or P = 984.3 d). Both the very high binary mass function of the system (f(M) = 10.54 M_⊙ or f(M) = 14.46 M_⊙) and the results of the spectral disentangling indicate that this system hosts a BH candidate. Further data in different wavelength ranges will be crucial in constraining the nature of this and the other systems.

These findings confirm that we need to continue working in this direction for two main reasons. In fact, determining the orbital solutions of the detected binary systems in NGC 1850 and other clusters will allow us (1) to map the orbital parameter distributions of these stars so that the binary fractions we measure no longer need to rely on an adopted parent orbital parameter distribution as done in this work and in all other literature studies so far; (2) to investigate the presence of peculiar systems hosting dark objects such as NSs and BHs, and finally constrain their populations in clusters, so far totally unknown. The observational campaign we are conducting with MUSE will see an extension in the future, with new data and new clusters to be added to the list with the aim of shedding more light on these questions.

ACKNOWLEDGEMENTS

SS acknowledges funding from STFC under the grant no. R276234. SK acknowledges funding from UKRI in the form of a Future Leaders Fellowship (grant no. MR/T022868/1). MG acknowledges support from the Ministry of Science and Innovation (EUR2020-112157, PID2021-125485NB-C22, CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033) and from AGAUR (SGR-2021-01069). TS acknowledges support from the European Union’s Horizon 2020 under the Marie Skłodowska-Curie grant agreement No 101024605. VHB acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through grant RGPIN-2020-05990.

DATA AVAILABILITY

The data underlying this work are available in Table A1 of the paper. The MUSE observations analysed here are stored in the ESO archive, while the WAGGS images can be retrieved from the WiFeS website. The radial velocity measurements underlying this study will be shared upon reasonable request to the authors.

Footnotes

References

APPENDIX A: ADDITIONAL MATERIAL

This Appendix contains supplementary material (Table A1 and Figs A1, A2, and A3) which help us to better understand the content of the paper, without interrupting its flow.