© ESO 2021

1. Introduction

Observations have shown that a majority of OB stars, which are stars with initial masses higher than 8 M_⊙, are members of binary or higher-order multiple systems (see, e.g., Abt et al. 1990; Shatsky & Tokovinin 2002; Sana et al. 2012, 2014; Rizzuto et al. 2013; Kobulnicky et al. 2012; Dunstall et al. 2015; Moe & Di Stefano 2017). A high fraction of them will interact during the course of their life (Sana et al. 2012, 2013; Dunstall et al. 2015), which will drastically alter the subsequent evolution of both binary components (Paczyński 1967; Podsiadlowski et al. 1992).

In contrast to the evolutionary pathway of single stars, interactions in binary systems can lead to a multitude of post-interaction pathways. This depends most strongly on the initial mass of the primary and on the initial orbital properties of the system, especially on the initial mass ratio and orbital period (Podsiadlowski et al. 1992; Langer 2012). Very tight binaries are thought to undergo contact evolution and eventually merge (see, e.g., Pols 1994; Wellstein et al. 2001; de Mink et al. 2009). In the case of stable mass transfer via Roche-lobe overflow (RLOF), the primary is stripped of its envelope (Kippenhahn & Weigert 1967). The companion accretes material and angular momentum and is spun up to close to critical spin rates, depending on the amount of material transferred, and the ability of tides to spin it down (e.g., Pols et al. 1991; de Mink et al. 2013).

Binary interactions thus lead to a variety of binary interaction products (BiPs) that can have very different observational characteristics. If originally part of a coeval star formation event, for example, in a single-starburst stellar cluster, mass gainers and merger products can appear as rejuvenated, single stars, some of which are situated above the cluster turnoff (TO). Observationally, stars above the TO are called blue stragglers (e.g., Ferraro et al. 1997) or red stragglers if the BiP evolves to the red supergiant phase (Britavskiy et al. 2019; Beasor et al. 2019). In the case of RLOF, the secondary is spun up and can be observed as a rapidly rotating star (Pols et al. 1991), whereas the stripped core of the original primary is hard to detect (Wellstein et al. 2001; Götberg et al. 2017). If the remaining mass of the stripped star is ≲1 M_⊙, it will evolve into a subdwarf O or B star (sdO/sdB; see, e.g., Gies et al. 1998; Heber 2009; Peters et al. 2013; Wang et al. 2017, 2018, 2021; Chojnowski et al. 2018; Götberg et al. 2018). If the core of the stripped star is massive enough, it will end in a core-collapse supernova, leaving behind a neutron star or black hole. Most systems will unbind as result of the natal kick of the compact object, creating a walkaway or runaway star (e.g., Eldridge et al. 2011; Renzo et al. 2019). Systems that remain bound are possible progenitors for X-ray binaries (Reig et al. 1997; Haberl & Sturm 2016), depending on their configuration and evolutionary stage.

Because a high fraction of massive stars interact with a companion during their lifetime, the number of BiPs in a given population of massive stars is expected to become significant (de Mink et al. 2014; Britavskiy et al. 2019; Wang et al. 2020). Theoretical computations have shown that the fraction of BiPs as well as the currently observed binary fraction varies as a function of the age of the population (Schneider et al. 2015).

In order to gain insight into initial multiplicity properties, the binary properties of massive stars in young clusters both in our Galaxy and in the Large Magellanic Cloud (LMC) have been well studied through spectroscopy (e.g., Kobulnicky et al. 2012; Sana et al. 2012, 2013). These spectroscopic studies are typically sensitive to systems with periods up to about 10 years, with a decreasing sensitivity above some years. In this paper, we adopt the following nomenclature: Systems detected through spectroscopy are called spectroscopic binaries , while the corresponding intrinsic binary population after bias correction is called close binaries f_cl. Unless specified otherwise, we adopt a maximum period of 10^3.5 days to allow a direct comparison with previous studies.

Investigating O stars in young open galactic clusters, Sana et al. (2012) found a bias-corrected close binary fraction of f_cl = 69 ± 9%. This is similar to the bias-corrected close binary fraction of OB stars in Cygnus OB2 of ∼55% (Kobulnicky et al. 2014). The bias-corrected close binary fraction of young O-type stars in the 30 Dor region in the LMC was found to be 51 ± 4% (Sana et al. 2013; Almeida et al. 2017). Focusing on young B-type stars, Raboud (1996) estimated a lower limit of 52% on the observed close binary fraction in the young Galactic cluster NGC 6231. In a recent study, Banyard et al. (2021) revisit the B-star population of NGC 6231 and find a bias-corrected close binary fraction of f_cl = 44 ± 6%. In the LMC, Dunstall et al. (2015) investigated B stars in the young cluster 30 Dor and found f_cl = 58 ± 11%.

Homogeneous studies of the multiplicity of massive stars in intermediate-age clusters are still lacking, however. In clusters of ∼20−40 Myr, O stars already have evolved off the main sequence (MS) and B-type stars now dominate the cluster TO (Ekström et al. 2008; Brott et al. 2011). Furthermore, a significant fraction of the initial binary systems will have had time to interact so that the number of BiPs in these clusters is expected to be significant (Schneider et al. 2015).

Investigating the current multiplicity properties of intermediate-age clusters as well as the frequency and characteristics of BiPs thus provides an ideal laboratory for testing the current theories of binary evolution, in particular the physics controlling the interactions. One such intermediate-age cluster is NGC 330 in the Small Magellanic Cloud (SMC), which is reported to be between 25 and 45 Myr old (e.g., Sirianni et al. 2002; Milone et al. 2018; Patrick et al. 2020). While massive stars in the SMC field are reported to have a metallicity of [Fe/H] = −0.7 (e.g., Luck et al. 1998; Korn et al. 2000; Keller & Wood 2006), NGC 330 was suggested to have an even lower metallicity of [Fe/H] ≲ −1.0 (Grebel et al. 1992; Gonzalez & Wallerstein 1999; Piatti et al. 2019). The high stellar density in NGC 330 prevented a spectroscopic study of the cluster core with previous instrumentation. The advent of adaptive optics (AO) at the Multi-Unit spectroscopic Explorer (MUSE) in 2017 changed the situation (Bacon et al. 2010).

Bodensteiner et al. (2020a, hereafter Paper I) focused on a single epoch of AO-supported MUSE observations of the SMC cluster NGC 330 (epoch 4, see Table 1). Paper I described the observations, data reduction and sky subtraction. Furthermore, we gave a detailed explanation of how spectra for individual stars were extracted from the MUSE data cube using point spread function (PSF) fitting. Using an automated spectral classification based on equivalent widths of He lines, we assigned spectral types to more than 200 stars. We derived a total cluster mass of and an age of the cluster core of 35−40 Myr.

NGC 330 was previously reported to contain a large number of Be stars (see, e.g., Feast 1972; Grebel et al. 1996; Keller et al. 1999; Martayan et al. 2007). Based on Hubble Space Telescope (HST) photometry, Milone et al. (2018) reported that the color-magnitude diagram (CMD) of NGC 330 shows a split MS and a large population of Hα emitting stars, which they interpreted as Be stars. The colors of these stars indicate that they are cooler than the bulk of MS stars, therefore these stars predominantly populate a region offset from the MS region (see their Fig. 10). Using the MUSE observations, we confirmed that stars in this cooler region in the CMD are indeed Be stars. We furthermore found that the overall observed Be star fraction in the cluster core is 32 ± 3%. Around and above the cluster turnoff, it rises to ∼50 ± 10%, as was indicated before by Milone et al. (2018).

Be stars are observationally identified as B-type stars with strong Balmer-line emission (Jaschek et al. 1981). As this observational definition includes several objects of different physical nature (i.e., Herbig Ae/Be stars, B-type supergiants, and interacting binaries), the term “classical Be stars” is used to describe rapidly rotating non-supergiant B-type stars with a circumstellar decretion disk that gives rise to the characteristic line emission (Rivinius et al. 2013). The origin of the rapid rotation of classical Be stars, which is thought to be close to the critical breakup velocity, is still debated (McSwain & Gies 2005). According to the proposed single-star channels, classical Be stars could be born as rapid rotators (Bodenheimer 1995) or spin-up toward the end of their MS lifetime (see, e.g., Ekström et al. 2008; Hastings et al. 2020). Another proposed formation channel interprets classical Be stars as mass gainers in previous binary interactions (Pols et al. 1991; de Mink et al. 2013; Shao & Li 2014; Bodensteiner et al. 2020c). A way to test the relative contributions of these channels is to investigate the multiplicity properties of classical Be stars: If the single-star channels dominate their formation, Be stars are expected to have binary properties similar to those of normal B-type stars. If the binary channel is the dominating one, however, Be stars are not expected to be in close binary systems with MS companions, but to rather have stripped or evolved companions (which are hard to detect observationally) or to be runaways (Bodensteiner et al. 2020c).

We extend our analysis here to study the multiplicity properties of NGC 330 using six epochs of MUSE observations. This paper is organized as follows. Section 2 describes the multi-epoch observations, their reduction and the extraction of individual spectra. Section 3 focuses on the measurement of radial velocities (RVs), and we discuss the correction for observational biases in Sect. 4. We derive the spectroscopic binary fraction of the overall massive-star population of NGC 330 in Sect. 5 and discuss individual groups of stars of interest. Finally, in Sect. 6 we discuss our results, and the conclusions and an outlook to future work are given in Sect. 7.

2. Observations, sample, and extraction of spectra

2.1. Observations and data reduction

Six epochs of observations of NGC 330 were obtained with MUSE between August 2017 and December 2018. The six epochs were split into four dither positions of 540 s each and were spread out over approximately 1.5 years. A journal of the observations with seeing conditions is given in Table 1. Additionally, the average full width at half maximum (FWHM) of the point spread function (PSF) as measured by fitting 2D Gaussian profiles on several bright and isolated targets in the reduced and collapsed white-light image is provided for each epoch. This gives an indication of the final image quality after data reduction, taking the effect of the AO into account.

MUSE is made up of 24 individual spectrographs that in wide-field mode (WFM) cover a field of view (FoV) of 1′ × 1′ with a spatial pixel sampling of 0.2″ in both spatial directions. At the estimated distance of 60 ± 1 kpc (see, e.g., Harries et al. 2003), 1′ corresponds to about 20 pc and 0.2″ to 0.07 pc or 1.4 × 10⁴ AU. Observations in the extended wavelength mode cover the optical spectrum between 4600 Å and 9300 Å. When the AO mode is used, the region between 5780 Å and 5990 Å is blocked to avoid stray light from the laser guide stars. The spectral resolving power λ/Δλ varies from 1700 at 4600 Å to 3700 at 9300 Å.

All epochs were reduced as described in Paper I. In short, we used the standard ESO MUSE pipeline v2.6¹ to apply bias and dark subtraction, flat fielding, wavelength and illumination correction, telluric correction and flux calibration. A sky subtraction was performed by taking the darkest pixels in each observation to create an average sky spectrum that was then subtracted from each pixel in each data cube.

2.2. Extraction of spectra

The extraction of spectra from the fully reduced MUSE data cubes was made through via PSF-fitting. In order to take the crowding in the field into account, we fit the PSF of multiple, close-by sources simultaneously to limit contamination. For this purpose, we need a catalog of accurate stellar positions and brightness, preferably from an instrument with higher spatial resolution, which is then used as input for the PSF fitting.

In contrast to the procedure adopted in Paper I, we used the photometric catalog provided by the GALFOR project², described in Milone et al. (2018), as the target input catalog. The catalog, which is based on HST wide-field camera 3 (WFC3) observations taken in 2015 and 2017, contains coordinates as well as magnitudes in different filters (i.e., F225W, F336W, F656N, and F814W). No magnitudes are available for three of the red supergiants (RSG) and for one blue supergiant (BSG). Additionally, some bright sources fall in the gap between the individual WFC3 chips. We estimated their position and magnitude by comparison to similarly bright stars and manually added them to the input catalog.

After aligning the HST catalog with the MUSE coordinate system, stars in the MUSE FoV with F336W magnitudes brighter than 17.5 (Vega magnitudes, corresponding to stars with masses ≳5.5 M_⊙ on the MS) were selected for spectral extraction. For each of these 400 stars, the spectrum was extracted using a PYTHON routine based on PSF fitting (see Paper I for more details). In short, the PSF of each star was fit at each wavelength slice (i.e., 3701 times); all stars closer than 12 pixel and brighter than 18.5 in F336W were fit simultaneously. This ensures that the high level of crowding in the dense cluster core is taken into account as close-by sources are effectively deblended. The individual flux measurements per wavelength bin were then concatenated to a spectrum for each star.

When stars are too close to the edge of the FoV, the PSF-fitting routine fails. Spectra for these stars were extracted manually using QfitsView³. The extraction in QfitsView is similar to aperture photometry, that is, the flux of a star in each wavelength slice is determined by summing the flux within a certain aperture (which was fixed to 3 × 3 pixels here). Because the stellar density toward the edge of the MUSE FoV is lower, the sources here are more isolated, and contamination by close-by sources is less likely.

Because of a small spatial offset between the epochs, the different observations do not cover the exact same region in the sky. As a result, some stars close to the edge of the MUSE FoV are not observed in every epoch and therefore have fewer than six epochs available. Furthermore, the seeing conditions varied (see Table 1), therefore the quality of the individual spectra varies from epoch to epoch.

In total, we extracted at least four spectra for 400 stars brighter than m_F336W = 17.5 mag. For 55 stars, mainly faint stars in crowded regions of the FoV, the quality of the extracted spectra, that is, the continuum signal-to-noise ratio (S/N), was ≲80 and thus was deemed insufficient for a subsequent analysis. Furthermore, there are five foreground stars in the FoV that are classified as such because while their brightness is similar to that of the B-type stars, they show the spectrum of K-M type star. The spectra of the remaining 340 stars were automatically normalized (see Paper I). In total, this provides us with normalized spectra for 324 B-type stars, ten RSGs and six BSGs.

2.3. Spectral types

Investigating the Balmer lines of the non-supergiant B-type stars for the presence of emission lines in at least one epoch, we identified 115 Be stars (i.e., 35%). The majority of them are probably classical Be stars, that is rapidly rotating B-type stars with Balmer-line emission arising in a circumstellar decretion disk (Rivinius et al. 2013). As described in Sect. 1, a similar spectral signature is observed for young stellar objects, however, that is, Herbig Ae/Be stars (Waters & Waelkens 1998), as well as for magnetic stars (Petit et al. 2013) or interacting binaries (Horne & Marsh 1986). Based on our current data, it is not possible to distinguish between these natures. Classical Be stars, however, are expected to populate a distinct region in the CMD (see, e.g., Milone et al. 2018). Additionally, given the crowding of MUSE FoV, and despite the advanced spectral extraction technique used here, blending with another source and thus contamination in the Hα line cannot be fully excluded. In order to assess this, we counted the number of faint sources (i.e., with F814W magnitudes > 18) with bright emission-line stars close by, and find that there are only ten such objects.

For all other B-type stars, spectral types were assigned based on an equivalent width (EW) calibration of diagnostic spectral lines of standard stars, which was developed and described in Paper I. As in Paper I, Be stars as well as BSGs were classified compared to standard star spectra by visual inspection. To do this, we used the spectra extracted from epoch 4 because this epoch has the best image quality and therefore the lowest contamination.

Our sample of 340 remaining stars includes 208 B stars, 116 Be stars, six BSGs, and ten RSGs. The S/N per epoch varies between ∼80 for the fainter and ∼200 for the brighter MS stars, and the BSGs have a typical S/N of about 300 at 6000 Å.

3. Radial velocities and the multiplicity criteria

This section focuses on the multiplicity properties of the B and Be star populations and of the BSGs. The multiplicity properties of the RSG population were studied extensively in Patrick et al. (2020) based on a more extended observational data set that also includes RSGs outside the MUSE FoV.

3.1. RV measurements

In order to assess the multiplicity status of each star, we measured the Doppler shift of a set of spectral lines in all available epochs. RVs are measured following the procedure described in Sana et al. (2013), where Gaussian profiles are simultaneously fit to a set of spectral lines for all epochs. The shapes of the individual line profiles are kept constant across the epochs, and the different lines fit for each epoch are required to be reproduced by the same RV. The method is thus more robust than line-by-line fitting, especially for epochs with poorer S/N than the average.

For the B and Be stars, we used a set of six He I lines as standard line list: He Iλλ 4713.15, 4921.93, 5015.68, 5047.74, 6678.15, and 7065.19 Å⁴. We first visually inspected all spectra of each star and, if necessary, reduced the number of lines used to a subset of the standard line list. Lines were excluded depending on the S/N of the spectra, which particularly affects the weak He I lines at 5015 and 5047 Å. In addition, Be stars show line infilling by emission which can significantly reduce the set of lines that can be used for RV measurements. On average, four He I lines were used, but only one line is available for a handful of stars.

While spectra were extracted for 324 B and Be stars, RV measurements were only possible for 282 of them. This was mainly due to the low S/N in some of the spectra, either because the stars in question are too faint, or because they are close to significantly brighter sources, which adds to the photon noise of the target itself. For the 282 B and Be stars, we have six epochs for 273 stars, five epochs for eight stars, and four epochs for the remaining star.

For the six BSGs in the sample we chose a different set of lines to measure RVs. Because they are all early to late A-stars, their spectra do not show strong He I lines. We therefore used a set of Fe II and Si II lines (i.e., Fe IIλλ 4923.93 and 5018.44 Å, and Si IIλλ 5055.98, 6347.10, and 6371.36 Å). For stars 327 and 586, we excluded the Fe II at 4923 Å because it is blended with a weak He I line. For the late-A supergiant 503, we only used the Si II lines around 6350 Å because of the large number of contaminating spectral lines in the spectrum.

Our final catalog, including coordinates, magnitudes, number of epochs, spectral types, RV measurements, the set of lines used, and the binary flag we assigned to them (see Sects. 3.2 and 3.3) and individual notes is given in Table C.1. This thus provides a large-scale homogeneous and self-consistent catalog of RVs for almost 300 stars, the full version of which will be made available electronically at the CDS.

3.2. Multiplicity criteria

In order to assess whether the measured RV variability is statistically significant, which allows us to classify stars as binary candidates, we adopted two multiplicity criteria (see, e.g., Sana et al. 2013; Dunstall et al. 2015; Patrick et al. 2019; Banyard et al. 2021; Mahy et al. 2021). A star is classified as a binary if at least one pair of RVs measured in different epochs satisfies the following two criteria simultaneously. The first criterion is

(1)

where v_i and v_j are the individual RV measurements and σ_i and σ_j the respective RV errors at epoch i and j. The confidence threshold of 4.0 is chosen so that the number of false positives is smaller than one given the sample size and number of epochs. As the RV errors are typically smaller for brighter stars, these are more likely to pass this criterion than fainter stars with larger RV uncertainties. This effect, however, is accounted for in the bias correction (see Sect. 4 and Table B.1).

The second criterion takes into account that there are other processes that can lead to RV variability in massive stars, most importantly pulsations (Aerts et al. 2009; Simón-Díaz et al. 2017) and wind variability (see, e.g., Fullerton et al. 1996). We therefore imposed that two individual RV measurements have to deviate from one another by more than a minimum RV threshold:

(2)

We investigate the effect of the choice of this threshold ΔRV_min in Fig. 1. The overall observed close binary fraction for MS stars is low, and small changes in the threshold value will not change the results significantly. We adopted ΔRV_min = 20 km s⁻¹, which is in line with similar studies (see, e.g., Sana et al. 2013; Dunstall et al. 2015; Banyard et al. 2021) and above the typical RV variations of B-type stars caused by pulsations (e.g., Aerts et al. 2009).

Fig. 1. Observed close binary fraction for MS stars as a function of the adopted RV-variability threshold ΔRVmin. For each ΔRVmin we compute the observed close binary fraction of the entire sample as well as the corresponding binomial error. The vertical dotted line gives the adopted threshold value ΔRVmin = 20 km s−1 indicated by the red line.

According to our two criteria, 24 out of the 189 B stars for which we were able to measure RVs show significant (Eq. (1)) and large (Eq. (2)) RV variability. Of the 93 Be stars with measured RVs, five pass these RV criteria, as well as one of the six BSGs. These stars are considered spectroscopic binaries in the following, even though the orbital period could not be constrained.

3.3. Double-lined spectroscopic binaries

We identified additional binary candidates by visually inspecting the spectra of all the B and Be stars as well as BSGs, specifically looking for composite spectral signatures indicative of double-lined spectroscopic binaries (SB2s). These are not always flagged as binary by the applied RV criteria because their spectral lines are typically not well represented by a single-Gaussian fit resulting in inaccurate RV measurements with large error bars.

We focus on the He I lines described above that we used for the RV measurements, and on the Balmer lines. Because the diagnostic spectral lines of Be stars are often filled-in by emission, the visual classification of SB2s is particularly difficult. We find eight systems (54, 55, 119, 156, 203, 212, 252, and 637) that show a composite spectrum that is indicative of two stars and are therefore classified as candidate SB2s. Two of these are Be stars while the remaining six are B-type stars. While they might be genuine SB2s, they might also be the line-of-sight superposition of two sources because the high stellar density in the cluster core is high. The stars show line profile variations, which would not be expected for two single stars that are aligned by chance, therefore this seems unlikely. Our data quality does not allow us to derive accurate spectral types for the secondary components. We detect them based on He I lines, and they contribute significantly to the flux in the optical, therefore they are most likely B-type stars.

4. Correction for observational biases

In order to assess the intrinsic close binary fraction of the massive-star population in NGC 330, the observed spectroscopic binary fraction needs to be corrected for the observational biases. We specifically considered the following two detection biases. The first bias is the possibility of detecting significant RV shifts above our detection threshold, depending on the orbital properties, temporal sampling, and RV measurement accuracy. The second bias affects SB2 binaries, or more precisely, unidentified SB2 binaries. In these systems, the presence of the companion reduces the apparent Doppler shift of the primary star, which in turn reduces the possibility of detecting a Doppler shift that is strong enough to meet our detection criteria.

4.1. Single-line detection bias

Following the approach described in Sana et al. (2013), we estimated the binary detection sensitivity of our observing campaign by simulating a population of binary systems. Corresponding to the number of stars in our sample we measured RVs, we simulated 280 binary systems taking into account our observational setup (i.e., the time coverage and number of observational epochs). For each system, we randomly selected the initial mass of the primary M₁ from a Salpeter initial mass function (Salpeter 1955) between 4.5 and 8 M_⊙. The choice of the range of primary masses was guided by the mass distribution of the MS stars in our sample, estimated in Paper I. We then paired them with a companion by randomly drawing the mass ratio q, the orbital period P, and the eccentricity e from parent distributions of orbital parameters. We further adopted a random orientation of the orbital plane in 3D space and a random reference time in the ephemeris. We first assumed a flat period distribution in log P, ranging from 0.15 to 3.5 (i.e., P from 1.4 to 3160 days), a flat mass-ratio distribution between q = 0.1 and 1, and an eccentricity distribution between e = 0 and 0.9 that is proportional to , with a simple circularization correction for short periods. These parameter ranges were adopted to allow for a direct comparison with other works. The choice of the underlying parent distributions of orbital parameters is discussed in Sect. 4.3.

After assigning RV errors, taken from the observations, we applied the two multiplicity criteria described in Eqs. (1) and (2). By repeating the simulation of 280 binary systems 10 000 times, we computed the probability of detecting these systems as binaries with our observing campaign. The number of simulations was chosen so that the associated statistical uncertainty due to the number of simulations is approximately five times smaller than the one due to the sample size. This provides an estimate of the sensitivity of our observation, and provides us with a correction factor that is multiplied by the observed binary fraction in order to estimate the intrinsic close binary fraction. A second run was then performed for which we adopted the intrinsic binary fraction as input to allow for a proper estimate of the uncertainties due to the binary population sample size. We find that the overall detection probability over the period range of 1.4−3160 days is ∼45%. As shown in Fig. 2, the detection probability strongly depends on the period and mass ratio: The detected probability is almost 90% for the shortest period systems, but it drops below 50% for periods longer than ∼100 days. Two-dimensional bias correction plots are shown in Appendix B.

Fig. 2. Binary detection probability simulations following our observational setup. Top panel: detection probability as a function of the orbital period P (computed assuming a flat mass-ratio distribution), and bottom panel: this as a function of mass ratio q (computed assuming a flat logP distribution). In both panels, the black line shows our first estimate of the binary detection probability taking the SB1 bias into account, and the orange line shows the detection probability, which also includes the reduced sensitivity for detecting near equal-mass binaries. The colored arrows in the bottom panel correspond to the mass ratio for which we compute the SB2 bias (see Fig. 3).

4.2. Double-line detection bias

SB2 systems for which both components show similar line strengths show an additional detection bias: The apparent Doppler shift of blended spectral lines can be reduced or even masked entirely. The reason is that the superposition of two similar line profiles is hard to distinguish from a single-line profile if the RV separation is not large enough and the variation of the (photo-)center of this presumably single line is reduced compared to the true orbital Doppler shift (e.g., Sana et al. 2011). This SB2 bias is lifted as soon as the SB2 nature of the system can be identified and a double-line profile fit can be used to measure the RVs of both components simultaneously. This bias is therefore stronger for broader lines (either due to rotation or spectral resolution) and lower S/N.

Because of the low resolving power of the MUSE spectra (i.e., between 1700 and 3700), we performed more detailed simulations in order to investigate the minimum RV separation at which their spectral lines are deblended and the two stars would be classified as SB2s. This allowed us to quantify the bias that affects the RVs measurements of unidentified SB2 systems. We used synthetic spectra from the TLUSTY B-star grid computed for SMC metallicity (Lanz & Hubeny 2007) and selected atmospheric models corresponding to spectral types B2 to B5, following the relation of stellar parameters and spectral types from Silaj et al. (2014). We performed this simulation for a B2 V single star, and a set of binary systems with a B2 V primary and a mass ratio q ranging from 0.5 to 1.0. We furthermore broadened the spectra for two rotational velocities v sini of 100 and 200 km s⁻¹, which is typical for stars in the considered mass range. We shifted the stellar spectra in steps of 10 km s⁻¹ with respect to one another, coadded them accounting for their light ratio, degraded the simulated spectra to the MUSE resolution and a binning of 1.25 Å (corresponding to the binning of our reduced spectra), and adopted the MUSE wavelength coverage and a typical S/N corresponding to the observations.

We then measured RVs with the same method as for our sample stars. Figure 3 shows that we can accurately measure RVs for a single star and that the RV bias remains small for binaries with a low mass ratio, that is, q <  0.6, for which the composite spectra are dominated by the lines of the primary. For mass ratios closer to unity, that is when the two stars have nearly equal masses, both stars contribute significantly to the composite spectrum (the one with a blueshifted profile; the other with a redshifted profile). We also inspected the simulated spectra visually to investigate whether they would be classified as an SB2. Because of the resolution and S/N of our MUSE data, the SB2 lines are only sufficiently deblended for RV separations above 200 to 300 km s⁻¹ (depending on the mass and luminosity ratio). At lower RV separations, the measured RVs are strongly biased, resulting in smaller measured RV variations than in reality.

Fig. 3. Measured primary RVs for simulated spectra of one single B2 V star, and SB2 binary systems with a B2 V primary and six different mass ratios as indicated in the legend. A rotational velocity of v sini = 200 km s−1 is assumed for all stars. The horizontal axis gives the assumed RV difference between the two components. The dashed line indicates the true velocity of the primary star.

Taking these simulations and the difficulties of detecting near equal-mass binaries into account, we repeated the bias correction computations. To do this, we adopted the measured RVs from the simulated spectra experiment described above in the v sini = 200 km s⁻¹ case as this is representative of most of our stars (see a subsequent paper in this series). We smoothed and interpolated it on a fine grid ranging from 0 to 600 km s⁻¹ in steps of 1 km s⁻¹ in RV separation and from 0.5 to 1.0 in step of 0.01 in mass ratio. This was then used as a correction to the simulations described in Sect. 4.1, for which the computed primary RV shift in the simulated time series is reduced as a function of the separation between the two components and their mass ratio. Systems with q <  0.5 are assumed to be unaffected by this bias. Figure 2 displays our result and shows that the overall detection probability is reduced by about 10 to 20% compared to the SB1 case. This reduction remains small up to a mass-ratio of 0.8, but dramatically affects systems with a mass ratio in the range 0.8−1.0, where we can only detect 10 to 20% of the near-equal mass binary population (compared to a detection probability of about 60% without the SB2 bias).

The SB2 bias is particularly dramatic for our present MUSE data. This results from a combination of factors, including the modest spectral resolution, wavelength binning size, and S/N of the spectra, the typical (projected) rotation rates, the exponent of the mass-luminosity relation in the B-star regime, and the involved masses and corresponding orbital velocities. We recall that the orbital velocity of a 7 M_⊙ + 7 M_⊙ binary with a 2-day period is about 200 km s⁻¹. Temporal sampling, inclination, and longer periods will inevitably push similar systems in the range of velocity separations that is strongly affected by the SB2 bias. Even though it would be interesting to investigate the SB2 bias for previously performed O- and B-type binary surveys, we remark that the situation is more favorable for most of the surveys discussed in Sects. 1 and 6, either because of the higher resolution, S/N, or mass range considered, so that there is no reason to believe that these results are as dramatically affected by the SB2 bias as in the present case.

This is also reflected by the number of detected SB2 systems in this work compared to similar spectroscopic studies of O-type stars. In this study we find only 8 out of 282 systems to be SB2s. This is in stark contrast to almost 45% of SB2s found in a sample of galactic O star (Sana et al. 2012). The number of SB2s detected in Dunstall et al. (2015) and Banyard et al. (2021), who focused on B-type stars, is similar to what we find here, implying that the mass-luminosity relation in the B-star regime plays an important role in detecting SB2s.

4.3. Parent orbital parameter distributions

Integrating the detection probability for flat period and mass-ratio distributions over the range of parameters shown in Fig. 2 yields an overall detection probability of about 35% when both the SB1 and SB2 bias discussed in the previous sections is taken into account. Adopting different parameter distributions would result in different values, however. To assess the effect of these choices, we recomputed the overall detection probability for different input distributions. In this exercise, we adopted power-law representations of the orbital parameter distributions,

(3)

(4)

(5)

and we varied their indexes π, κ and η.

The shape of the period distribution is by far the dominant factor that affects the resulting detection probabilities. Fixing all other assumptions and changing π from 0.0 to −0.5 modifies the overall detection probability from 0.35 to 0.45 because the period distribution now favors shorter-period systems that are easier to detect. Varying κ and η from −0.5 to +0.5, respectively, has an overall effect of a few 0.001 and is thus negligible in this context.

To adopt meaningful parameter distributions, we considered measurements obtained using the same parameter ranges as those we adopt here. We refer to Almeida et al. (2017) for a discussion of the effect of the boundaries of the orbital parameter intervals on the measured power-law indexes of the orbital parameter distributions. O-star studies in the Galaxy and LMC yielded π = −0.55 ± 0.22 and π = −0.45 ± 0.39, respectively (Sana et al. 2012, 2013). The LMC value was revised to π ≈ −0.2 by Almeida et al. (2017) using a more comprehensive data set. Finally, Dunstall et al. (2015) obtained π = 0.0 ± 0.5 for early B-type stars in the LMC. We therefore adopted a value of π = −0.25 ± 0.25, which encompasses all measured values in its ±1σ interval, as a plausible representative value for the period distribution. Neither Sana et al. (2013) nor Dunstall et al. (2015) were able to reliably constrain κ and η, so that we adopted the values from the Milky Way O-star sample, that is κ = −0.2 ± 0.6 and η = −0.4 ± 0.2 (Sana et al. 2012), but we note that the latter choices have no effect on our results.

Because the period index π was proposed to be mass dependent (i.e., skewed toward shorter periods for O stars, and to slightly longer periods for late-B stars; Rizzuto et al. 2013; Moe & Di Stefano 2017; Tokovinin & Moe 2020), we tested how this mass-dependent period index affects our estimate of the detection probability. We therefore assumed a variable period index that changed linearly from π = −0.25 ± 0.25 at the upper mass end considered here to π = +0.25 ± 0.25 at the lower mass end. We find that the obtained detection probability is lower by 0.06 which remains within the error bars of a mass-independent detection probability computed using π = −0.25 ± 0.25 for the entire considered mass range (see Table 2). Therefore we adopted the latter for the following computations.

In the final set of bias correction computations, we thus drew 10 000 artificial populations of 280 binaries (corresponding to our sample size). For each population, we randomly assigned values of π, κ, and η from the parent orbital distributions following normal distributions with central values and 1σ dispersion taken as described above. The detection probability and, after a second iteration, its error bars, were then obtained as the average and root-mean-square of the fraction of detected binaries across the 10 000 simulated populations. In this sense, the error bars encompass the uncertainties that are due to the choice of orbital distributions and those due to the sample sizes. By running simulations with much larger population sizes, we estimate that the adopted spread in orbital parameter distributions amounts to about 6% in the final uncertainty.

Finally, we also computed a magnitude-dependent bias correction for which we picked primary masses corresponding to the different magnitude bins and typical RV errors measured for stars in that magnitude bin (see Appendix B), as well as in different regions of the CMD (see below). Our bias correction did not take the evolutionary history of binary systems and the nature of possible companions into account, but was purely based on the assumed distributions of orbital parameters. The adopted distributions are typical of young massive-star populations, therefore they are more likely to represent the initial conditions than more evolved populations that have been affected by binary interactions.

5. Multiplicity properties

The overall observed spectroscopic binary fraction, that is, of our entire sample of 282 stars with measured RVs, is . Applying the bias correction described in Sect. 4 leads to an intrinsic, bias-corrected close binary fraction of % for the massive-star population of NGC 330.

Combing the spectroscopic information from MUSE with the photometric catalog of Milone et al. (2018), Fig. 4 displays a CMD with binarity information. It shows that most of the spectroscopic binaries are situated on the MS. However, the Be stars include a handful of binaries, and several binaries are located above the MS turnoff, as indicated by the single-star nonrotating Padova isochrones (Bressan et al. 2012; Chen et al. 2014, 2015; Tang et al. 2014; Marigo et al. 2017; Pastorelli et al. 2019). There are no magnitudes in the relevant HST filters for the binary BSG, or for the four additional RSGs, so they are not included in the CMD.

Fig. 4. Hubble Space Telescope CMD (Milone et al. 2018) overlaid with binarity information. B-type stars are marked in yellow, Be stars in green, BSGs in blue, and RSGs in red. Open small circles indicate that no RVs could be measured, filled circles are stars with no significant RV variability, diamonds indicate stars classified as binary systems based on their RVs, and triangles indicate SB2 candidates.

Section 5.1 describes our results for different spectral (sub)types. In Sect. 5.2 we use the position of stars in the CMD to discuss the multiplicity of different evolutionary stages. A more detailed description of several interesting binary systems is given in Sect. 5.3.

5.1. Binary fraction of the massive-star population

Below, we briefly discuss the binary fraction of the B-type stars, Be stars, and BSGs in our sample. In this context, we refer to Be stars as all stars that show Balmer-line emission in their spectra (see also Sect. 5.2).

B stars. Of the 209 B stars in our sample, 165 stars do not show significant RV variations, while 24 stars are classified as binaries based on their RV variation. Of the 165 non-RV variable stars, 6 show composite line profiles indicative of SB2 systems. RV measurements were not possible for 20 stars. Combining the RV variables with the 6 detected B-type SB2s leads to an observed spectroscopic binary fraction of %.

Be stars. Ninety-three of the 115 Be stars in our sample were investigated with RV measurements. Only 5 Be stars pass our RV-variability criteria and are considered as binaries. Additionally, we classify 2 Be stars as candidate SB2s. Altogether, we find an overall observed spectroscopic binary fraction of % for the Be stars in the sample.

Blue supergiants. One of the six BSGs in the sample, namely 729, shows significant RV variations and is classified as binary following our multiplicity criteria. This yields an observed spectroscopic binary fraction of BSGs in the core of NGC 330 of . The large errors result from the low-number statistics.

Figure 5 shows the distribution of spectral types and the observed close binary fraction for B and Be stars. Spectral subtype bins with fewer than ten objects provide few statistical constraints (i.e., spectral types earlier than B2 and later than B8) and were not considered further. The remaining spectral sequence, from B2 to B7, shows a clear decrease in the observed binary fraction. This trend becomes clearer when the SB2s are included as well, but is also visible from the RV measurements alone. When the SB2s are included, the observed binary fraction decreases from 30% at B2 to less than 3% at B7. The limitations of our data make it hard to decide whether the increase in the binary fraction toward the latest spectral types, that is, B8 and B9, is real or due to the low-number statistics.

Fig. 5. Distribution of spectral types and observed close binary fraction for B and Be stars. Top panel: distribution of spectral types for all MS stars (dim gray) and stars for which RVs were measured (dark gray). Bottom panel: observed spectroscopic binary fraction as a function of spectral type. The red curves indicate the binary fraction inferred from the RV measurements, and the orange curve combines the RV measurements with the SB2s found by visual inspection.

As discussed in Paper I, the spectral classification of the Be stars in the sample is difficult. On the one hand, the He I lines used for spectral classification can be affected by emission infilling from the Be star disk. On the other hand, Be stars are rapidly rotating stars, which, depending on the inclination angle, may appear cooler due to their rotation. This might lead to a classification into later spectral types, which affects the spectral type distribution. We therefore also considered the observed spectroscopic binary fraction as a function of F814W magnitude, which is less sensitive to these effects.

Figure 6 shows the distribution of F814W magnitudes of the B and Be stars as well as their binary fraction. The distribution of F814W magnitudes indicates that our sample is not complete at F814W >  18.5 mag. As shown in the bottom panel, the decreasing trend in the observed binary fraction that was visible in the spectral type distribution is also visible here. Neglecting stars brighter than 15th magnitude, where the number statistics are low, as well as stars fainter than F814W = 18.5 mag, where the sample is incomplete, we find that the observed spectroscopic binary fraction decreases from 33% to 7%.

Fig. 6. Distribution of F814W magnitudes and observed close binary fraction for B and Be stars. Top panel: distribution of F814W magnitudes of all stars for which spectra were extracted (dim gray), and stars for which RVs were measured (dark gray). Bottom panel: observed spectroscopic binary fraction as a function of F814W magnitude. The orange curve is based on RV measurements, and the red curve is based on RV measurements as well as SB2s. The bias-corrected close binary fraction is overplotted in dark red.

We applied a magnitude-dependent bias correction to the observed binary fraction, as described in Sect. 4.3. We find that the bias-corrected binary fraction of the brightest stars with F814W = 15 mag, which corresponds to about 8.5 M_⊙, is almost 60%. It drops to about 25% for the faintest stars considered here, which roughly corresponds to 5 M_⊙. While uncertainties in each magnitude bin are large due to smaller sample sizes, the downward trend remains. This means that it cannot be explained by the observational biases alone. At least part of it is likely due to genuine differences in the multiplicity properties. This decrease in the close binary fraction can either be due to a genuine lack of binaries at later spectral type or to a shift in the period distribution toward longer orbital periods.

In order to further investigate whether this decreasing trend in the bias-corrected close binary fraction is statistically significant, we performed a linear regression between F814W = 15 mag and F814W = 18.5 mag. Based on this, we can reject the null-hypothesis that there is no trend in the data at the 5% level. Our data are insufficient to decide whether the decrease in the close binary fraction with mass is linear or has a different functional form. A visual inspection of Fig. 6 may indeed suggest a step at F814W ≃ 17.5, corresponding to ≃6.5 M_⊙ on the MS.

5.2. Binarity in different regions of the CMD

In order to compare the observed close binary fraction in different evolutionary stages, we investigated the binary properties of the massive-star population from a different angle, that is, by defining different regions in the CMD based on expectations from theoretical considerations. In addition to the BSGs and RSGs, we used the Padova isochrone at 40 Myr to define four additional regions in the CMD, as illustrated in Fig. 7:

• 1.

  the MS region, which is an extended band with a width of about 0.4 mag around the 40 Myr isochrone;

• 2.

  the Be star region, which appears as a second group, redward of the MS, with a color offset of about 0.3 mag in Fig. 7;

• 3.

  the turnoff (TO) region, which is located around the cluster turnoff of the 40 Myr Padova isochrone;

• 4.

  the blue straggler star (BSS) region, which is above and to the left of the cluster turnoff.

While the BSS and the Be star region are clearly identifiable from the CMD, the definition of the TO region is more subjective. Furthermore, the assignment of stars to different groups is purely based on the position in the CMD. This means, for example, that a star that is observationally classified as a Be star, situated in the blue straggler region, is assigned to the blue straggler region. Similarly, we classified all stars as Be stars that are in the Be star region, even when they did not show emission in our spectra (this is the case for only five targets at the faint end of the Be star region). One star, which is classified as a binary based on RV variability, lies to the left of the MS. This source is discussed in Sect. 5.3.

Fig. 7. Binarity of different stellar groups, based on their position in the CMD. MS stars are marked in orange, stars in the Be-star region in dark green, stars around the cluster TO in red, and blue stragglers in blue. Additional HST sources for which no spectra were extracted are marked with gray dots, and stars without RV measurements are indicated with open circles. Colored filled circles indicate that the star is single, colored diamonds indicate that the star is classified as a binary based on RV variability, and colored triangles indicate SB2s. The individual systems of interest discussed in Sect. 5.3 are labeled.

Table 2 presents the spectroscopic and bias-corrected close binary fraction for each of these groups separately, including binaries detected through RV variations or through their SB2 line profiles. We also list the binary properties of the BSGs (see Sect. 5.1). The RSGs were studied by Patrick et al. (2020), who probed a different, complementary parameter space to ours. In addition to improving on the time coverage, they further used multi-epoch ESO-HARPS spectra (with a resolving power of R ∼ 80 000). Additionally, their study included five RSGs that are outside of the MUSE FoV.

For each of the regions in the CMD, we computed the binary detection probability, taking the approximate masses of the stars into account that populate different regions as well as typical RV errors as a function of F814W magnitude and color. The difficulty of measuring RVs for Be stars, whose spectra are dominated by emission, becomes apparent here: while they populate a similar brightness range as TO and MS stars, the detection probability is significantly lower for them, especially compared to the brighter TO stars.

In this exercise, no attempt was made to account for the effect of stellar evolution on the distribution of orbital parameters, that is, all simulations were performed adopting the orbital parameter distributions derived from young populations (see Sect. 4.3). Populations that are most affected by binary interaction, that is, Be stars and BSSs, are likely to deviate from these assumptions (e.g., Langer et al. 2020; Wang et al., in prep.).

Table 2 provides the bias-corrected close binary fraction for the entire sample (which includes the caveat described above) as well as for stars in the MS and TO region. For the latter, we provide the detection probability computed with our standard computation assuming a constant index π, as well as assuming a mass-dependent π as discussed in Sect. 4.3. For the Be stars and BSSs, we only provide the detection probability and not the bias-corrected binary fraction, so that the results presented in Table 2 still provide a valuable relative comparison of the effect of variable RV measurement accuracy for different categories of stars given their respective observational specificity discussed in Sects. 2 and 3. The different spectroscopic and intrinsic close binary fractions of the different groups are discussed in Sect. 6.

5.3. Individual systems of interest

The binary BSG. Star 729, or Cl* NGC 330 ARP 35, was first studied by Arp (1959) and classified as B9I by Feast (1964). In Paper I, we classified it as early A Ia-Ib. Feast & Black (1980) measured RVs and found values between 134 and 154 km s⁻¹, which are on a similar order of magnitude than the RVs we measure (i.e., a mean value of 152.6 km s⁻¹). Based on their data, they did not classify it as (candidate) binary. Using the MUSE observations, we find that star 729 is clearly RV variable and according to our two criteria, we classify it as a binary. Because we only have six epochs, it is not possible to derive orbital parameters such as the period of the system. We see no contribution of the secondary in the spectrum. Because star 729 was not included in the HST input catalog (probably because of saturation), we cannot add it to the CMD in Fig. 4.

A Be X-ray binary candidate. Paper I reported on the detection of a high-mass X-ray binary detected in the MUSE FoV (Shtykovskiy & Gilfanov 2005; Sturm et al. 2013; Haberl & Sturm 2016). Using the multi-epoch data, we investigated the binary status of the two stars (i.e. stars 94 and 97, both labeled in Fig. 7) that are within the 0.7″ error on the position of the X-ray source (we note that no spectrum could be extracted for star 97 in Paper I while star 113 discussed there is farther away from the position of the reported X-ray source). Both are classified as Be stars. While star 97 shows no significant RV variations, we cannot measure RVs following the procedure described in Sect. 3 for star 94 because all available He I absorption lines used for RV measurements are filled in with emission. A preliminary RV estimate, using the Balmer emission lines in star 94 shows significant RV variations, however. Based on the six epochs of MUSE observations, it is unfortunately neither possible to derive an orbit nor to prove association of either of the two Be stars with the X-ray source. Follow-up observations will further investigate this possible Be X-ray binary.

An O+B binary. One of the few O stars in the sample, star 119, classified as such because of the presence of He II absorption, is also an SB2. The spectra of star 119 are composite and show line-profile variations indicative of a B star companion. Additionally, the Hα line shows significant infilling, while the Hβ line only shows absorption. There is a Hα-emitting source close by, therefore this could be contamination. It could also imply the presence of a disk, either around one of the two stars or around the binary system. According to the position in the CMD (see Fig. 7), the system is in the blue straggler region. This region is populated by objects thought to be rejuvenated in binary interactions, and currently interacting systems (Wang et al. 2020). Because we detect it as a close binary system, and because in the case of ongoing mass transfer, material may be lost from the system, it is conceivable that star 119 has a circumbinary disk. More observations are required to further investigate the nature of this system.

A binary system hotter than the MS. Star 604, marked with a yellow diamond in Fig. 7, stands out in the CMD. It lies to the left of the MS, which implies that it is significantly hotter than stars of similar brightness. The spectral type we determine from the MUSE data is approximately B5, implying that the star is not particularly hot. This is based on a rather low S/N (between 60 and 85 in the different epochs), however, because the star is relatively faint (F814W = 18.8 mag). Because the star is close to the edge of the MUSE FoV, spectra could only be extracted in five of the six epochs. The average RV measured for this star is 58 ± 8 km s⁻¹ and thus differs significantly from the mean RV of the SMC of 153.7 ± 1 km s⁻¹ (Patrick et al. 2020), that is, by almost 100 km s⁻¹. The position in the CMD as well as the RVs, which differ from the RV of the SMC, may imply that star 604 is a foreground star. Magnitudes like this are also indicative of a hot, envelope-stripped companion (e.g., Götberg et al. 2018), however. We classify star 604 as a binary system based on the measured RV variations. In the MUSE spectra, we see no indication of a secondary, and further observations with higher resolution and better S/N are necessary to test this hypothesis.

6. Discussion

6.1. Different binary fractions in different regions of the CMD

As summarized in Table 2, we find that the spectroscopic binary fraction varies strongly in different regions of the CMD, that is, as a function of stellar mass and evolutionary status. While we did not compute the bias-corrected binary fraction of Be stars and BSS, the bias-corrected binary fraction for a given set of parent orbital parameter distributions is about 25% on the lower MS and increases to approximately 60% around the cluster turnoff. In the more evolved evolutionary phases, the BSGs and the RSGs, the binary fraction is 28% and 30%, respectively. However, as noted above, the orbital periods for RSGs probed by Patrick et al. (2020) are quite different. These changes are statistically significant and indicate a change in the binary fraction of the population, or in its distribution of orbital periods.

While it is difficult to disentangle the effects of initial mass and evolutionary status, we compared our findings to expectations from binary evolution theory (see, e.g., Wang et al. 2020). Most binary systems below the cluster TO are expected to be pre-interaction systems. In the hot part of the Hertzsprung-Russel diagram (HRD), binary interactions preferentially occur around the TO, that is, when the stars expand more significantly toward the end of their MS lifetime. Around the TO, the number of currently interacting systems is expected to be high, while the BSS region is populated entirely by post-interaction systems. As discussed in the introduction, Be stars are often interpreted as binary interaction products (Pols et al. 1991; Shao & Li 2014; Bodensteiner et al. 2020c), implying that their companions are expected to be either envelope-stripped stars (Abdul-Masih et al. 2020; Shenar et al. 2020; Bodensteiner et al. 2020b; El-Badry & Quataert 2020, 2021) or compact objects (Langer et al. 2020; Schürmann et al., in prep.). Both types of objects are difficult to detect spectroscopically, and most BiPs would not be captured in our bias correction due to the nature of the companions and the expected long orbital periods (Sect. 4).

MS region. The vast majority of stars on the MS are thought to either be single stars or binary stars that likely have yet not interacted. Their companions are also most likely MS stars, and we find a bias-corrected close binary fraction of %. On higher resolution data, depending on their mass ratio and S/N, these MS companions would be relatively easy to detect. Moreover, we have shown in Sect. 4.2 that our MUSE data show a significant detection bias against near-equal mass SB2 systems. This may explain the lack of detected binaries close to the red edge of the extended MS (the so-called binary line, Milone et al. 2018), where binary systems with a mass ratio close to unity are expected.

TO region. Around the cluster TO, the observed binary fraction reaches 33 ± 7% and % after bias correction. The detection capability of our campaign is better for stars in the TO region than for MS stars because their S/N and masses are higher. However, this is insufficient to explain the large difference compared to MS stars, which is over a factor of two, even after bias correction (see Table 2). As described above, binary interactions preferentially occur toward the end of the MS lifetime of stars which would allow them to remain longer around the TO than if they were single stars. The region around the TO is thus probably mainly populated by two groups of binaries: pre-interaction binaries (e.g., near-equal mass binary systems that appear to be brighter than single MS stars), as well as currently interacting systems (e.g., slow case-A binaries). Additionally, there might be a difference of intrinsic multiplicity properties between the TO and MS regions owing to their different mass ranges (∼7 M_⊙ vs. ∼5 M_⊙). On the one hand, the binary fraction increases with spectral types (e.g., Moe & Di Stefano 2017). On the other hand, solar-like stars typically have longer orbital periods than massive OB stars. This means that somewhere in between the two mass regimes, a transition should occur that mostly affects the period distribution, hence the detection probability of our survey.

Be star region. The low observed spectroscopic fraction of binaries in the Be-star region, which is only 2%, could have two explanations: Either the close binary fraction of classical Be stars is intrinsically low or Be binaries have long orbital periods, for which our detection probability drops (see Fig. 2). These findings are compatible with qualitative expectation of the binary scenario. Our RV study cannot distinguish post-interaction Be products from genuine single Be stars. Our simulations show, however, that if Be stars have an MS companion with a mass and period distribution similar to the MS B-type star population, we would have detected a significant fraction of Be binaries. In this sense, our results are also compatible with the conclusions of Bodensteiner et al. (2020c) that there is a dearth of Be+MS binaries.

BSS region. Stars brighter and hotter than the cluster TO, that is, blue stragglers, are typically interpreted as binary interaction products. They are either merger products, stars that were rejuvenated as mass gainers, or currently interacting systems. While merger products are now single stars (disregarding the possibility of initial triple systems), currently interacting (Algol) systems have a high probability of being flagged as binaries in our study. As mentioned before, because the BSS region is populated by BiPs, we did not apply a bias correction here. We find that the observed spectroscopic binary fraction of the BSS and TO populations are similar in our study, however.

6.2. Comparison to the outskirts of NGC 330

Our MUSE observations target the inner 1′ × 1′ core of NGC 330. The OB star population in the outskirts was investigated with VLT/FLAMES in Evans et al. (2006) down to V = 17 mag, reporting an observed spectroscopic binary faction of only 4%, which is considerably smaller than our overall fraction of 13% in the cluster core. As described in Paper I, the overlap between the two samples in terms of brightness is limited. More specifically, there are only 27 stars in our MUSE sample that are brighter than V = 17 mag, indicating that our comparison may show low-number statistics.

When we restrict our sample to the 27 stars that are brighter than V = 17 mag, we find that the observed spectroscopic binary fraction in the cluster core is %. This is in strong contrast to the 4% reported by Evans et al. (2006) for the outskirts. This might be due to the relatively short time coverage of the spectra from Evans et al. (2006) of only ∼10 days, however, which means that their study is only sensitive to the shortest orbital periods.

6.3. Comparison to other clusters

Comparing the close binary fraction in NGC 330 to the fraction measured in other clusters allows us to investigate the effects of two important additional cluster properties: the cluster age and the environment metallicity. Their respective effect on the multiplicity properties remains unclear. As described in the introduction, previous multiplicity studies have mainly focused on O- and early B-type stars in the MW and LMC environments. These OB stars are only present in young clusters. As discussed in Sect. 4.2, the SB2 bias investigated here has a less significant effect on the results of these works because the resolution of the data is higher.

Kobulnicky et al. (2014) investigated the binary properties of massive stars in the 3−4 Myr old Cygnus OB2 association in the Milky Way. Focusing on a sample of around 130 O and early-B type stars (spectral type earlier than B2), they observed a close binary fraction of ∼35%. Correcting for observational biases, they reported an intrinsic binary fraction close to 55%, which is higher than what we find for NGC 330. The reason probably is that the stars included in their sample are significantly more massive, and the binary fraction increases as the mass of the primary increases (see, e.g., Moe & Di Stefano 2017, and Fig. 6).

Banyard et al. (2021) studied the B-star population of the young, galactic cluster NGC 6231, which has an estimated age between 4 and 7 Myr. They reported an observed spectroscopic binary fraction of 35 ± 4%, which, corrected for observational biases using the same method as described here leads to an intrinsic close binary fraction of 44 ± 6%. This is slightly higher than the bias-corrected binary fraction we measure in NGC 330, also when only MS+TO stars are considered. There are two possible reasons for this: First, the population of stars studied by Banyard et al. (2021) is more massive. Alternatively, it might be due to the age difference of the two clusters. While a young cluster such as NGC 6231 allows us to probe close to the initial binary properties, NGC 330 is at an age where a significant number of binary systems had time to interact. As described above, our bias correction focused on MS companions and did not correct for post-interaction systems, which are extremely difficult to detect in observations.

Dunstall et al. (2015) measured the binary fraction of B stars in the 30 Doradus region in the LMC. Focusing here on the non-supergiant population, they reported an observed spectroscopic binary fraction of 25 ± 2%. Additionally, they provided measurements for individual clusters of different ages (Evans et al. 2015): They observed spectroscopic binary fractions of 27 ± 5% and 34 ± 5% for the younger clusters NGC 2070 and NGC 2060 (age ≃ 2−8 Myr, Schneider et al. 2018), respectively. The observed spectroscopic binary fractions reported for the slightly older clusters Hodge 301 and SL 639, with estimated ages of 15 ± 5 Myr and 10−15 Myr (Evans et al. 2015), are only % and 10 ± 9%. They estimated the binary detection probability of B-type stars to be 0.4. Correcting the observed binary fraction for the detection probability leads to an intrinsic close binary fraction of 58 ± 11% for the overall B-star population and to an intrinsic close binary fraction of 20.0 ± 20% and 25.0 ± 22.5% in Hodge 301 and SL 639, respectively. Despite the low-number statistics in Hodge 301 and SL 639 (i.e., 22 stars in total), the intrinsic binary fractions are very similar to the bias-corrected binary fraction of MS stars reported here (i.e., ).

As mentioned before, the older age of NGC 330 as well as the lower metallicity of [Fe/H] ≲ −1.0 might be the reason for the lower inferred binary fraction in NGC 330. Several previous studies have found that for massive stars, the close binary fraction as well as the orbital parameter distribution seem universal over different metallicities (see, e.g., Sana et al. 2013; Moe & Di Stefano 2013; Sana 2017; Almeida et al. 2017; Banyard et al. 2021). We cannot independently measure the orbital parameter distributions of the B-type stars in NGC 330 with the data discussed in this work, therefore our bias-corrected close binary fractions are based on the assumption of those. Our current data set therefore does not allow a clear distinction between the two scenarios.

Future studies comparing the results of NGC 330 to clusters of similar age in different metallicity environments would allow us to investigate any possible metallicity effect. Furthermore, determining the orbital solutions of the detected binary systems in NGC 330 would allow us to map the orbital parameter distributions of B-type stars at SMC metallicity. These would allow us to investigate whether the lower binary fraction in NGC 330 is indeed an age effect or due to the lower metallicity. The binary fraction of MS+TO stars obtained here agree in general with the binary fractions found for Galactic and LMC samples, however.

7. Summary and future work

In this second paper of a series, we used multi-epoch MUSE observations in wide-field mode and supported by AO to investigate the multiplicity properties of the massive-star population of the ∼40 Myr old, massive, open cluster NGC 330 in the SMC. We extracted spectra for 400 stars and measured RVs of almost 300 stars using Gaussian fitting of spectral lines. Based on two RV variability criteria, we classified stars into single or binary stars. We furthermore visually inspected the spectra of all stars in order to detect composite spectra and line profile variations indicative of SB2s. With this, we provide the first large and homogeneous RV study of B-type stars in a metal-poor environment.

We find that the overall observed spectroscopic binary fraction of the massive-star population in NGC 330, including B stars, Be stars, and BSGs, is . To correct for observational biases, we simulated two effects: the overall binary detection sensitivity of our observing campaign based on the amplitude of the primary RV motion (called SB1 bias), and the difficulty of detecting binary systems with similar components (referred to as SB2 bias). Based on mainly the low-resolution of the MUSE spectra and the masses of the stars considered here, we find that the SB2 bias is sizeable in this study, and probably much more important than in previous studies that implemented the SB1 bias alone (e.g., Sana et al. 2012, 2013; Dunstall et al. 2015; Banyard et al. 2021).

Correcting for these observational biases, we find a bias-corrected close binary fraction of % for the massive-star population in NGC 330. Considering the B and Be stars, we find that the observed and bias-corrected close binary fraction is a function of F814W magnitude: The bias-corrected close binary fraction drops from ∼55% to ∼25% in the considered magnitude range. While the brightest stars are above the cluster turnoff, and thus are most strongly affected by binary interaction, this might provide indications that the binary fraction drops as a function of mass. We cannot distinguish these two effects based on the current data, however.

Investigating the binary fraction of different evolutionary phases of stars, indicated by their position in the CMD, we find that the binary fractions of the subpopulations, that is, MS stars, stars around the TO, BSS, and Be stars, are different. Because some of these populations, in particular BSSs and Be stars, might be dominated by BiPs, we did not apply our bias correction to these groups. While the observed spectroscopic binary fraction of MS stars is about 10%, the observed spectroscopic binary fraction around and above the cluster turnoff is significantly higher, that is around 30%. The observed spectroscopic binary fraction of star in the Be star region is only ∼6%. These observational results are compatible with qualitative expectations from binary evolution and set the stage for a detailed comparison with predictions from binary population synthesis.

Comparing the overall bias-corrected close binary fraction of NGC 330 with the binary fractions determined for B stars in clusters in the LMC as well as in the Milky Way (Kobulnicky et al. 2014; Dunstall et al. 2015; Banyard et al. 2021), we find that it is generally lower. This is probably a combination of two main effects that are difficult to disentangle. First, the binary fraction might depend on the metallicity of the host environment. Second, the clusters studied previously are very young (i.e., only a few million years) and are thus more representative of the initial multiplicity conditions. Additionally, while focusing on the B-star population, previous works mostly focused on early-type B stars, while most stars in this study are around spectral type B5, that is, they have lower masses.

In a subsequent paper, we will use the derived RVs to combine the six observational epochs for each star in order to boost the S/N in the observations. These will then be used to obtain stellar parameters, such as effective temperatures, surface gravities, and in particular rotational rates. Investigating these across the CMD will allow us to further investigate the physical properties of the cluster members and to identify additional possible binary interaction products.

Acknowledgments

The authors would like to thank the referee, Max Moe, for his constructive feedback on the manuscript. The authors acknowledge support from the FWO_Odysseus program under project G0F8H6N and from the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) through the PRODEX Programme. The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement numbers 772225: MULTIPLES). Parts of the analysis in this project are based on the python code PHOTUTILS.

References

Appendix A: Two-dimensional bias correction plots

Figures A.1 and A.2 show the 2D representations of the bias correction calculated for the SB1 and as well as the SB1+SB2 bias that were computed with flat period and mass-ratio distributions and a eccentricity distribution. The comparison of these plots reveals the strong effect of the SB2 bias for all binaries with q >  0.8, but for those with the shortest periods. The staircase aspect of the e − log P plots results from our simplistic pseudo-circularization prescription that aims to avoid unrealistically high eccentricity at short periods as stars in those systems would collide at periastron.

Fig. A.1. Binary detection probability as a function of period, mass of the primary M1, mass ratio q and eccentricity e, computed for the SB1 bias.

Fig. A.2. Same as Fig. A.1, but computed taking both the SB1 and the SB2 bias into account.

Appendix B: Magnitude-dependent bias correction

As explained in Sect. 4, in addition to the overall bias correction for the entire sample, we computed a magnitude-dependent bias correction. For this, we repeated our simulation including the SB1 and SB2 bias for different magnitude bins between F814W = 14.5 and 19 mag in steps of 0.5 mag while changing the primary mass as well as typical RV errors accordingly. The magnitude-dependent detection probability as a function of initial periods P as well as mass ratios q is shown in Fig. B.1. The primary mass as well as the overall detection probability for the SB1+SB2 bias is given in Table B.1 for each magnitude bin.

Fig. B.1. Binary detection probability as a function of period (top panel) and mass ratio (bottom panel). The different dotted orange curves correspond to different magnitude bins, in which the detection probability is highest for the brightest stars. The continuous orange line shows the overall detection probability for the full simulation, and the continuous black line shows the detection probability for the SB1 bias alone (see Sect. 4.1).

Appendix C: Tables with radial velocity measurements

Table C.1 gives an overview over the coordinates, F814W magnitudes, spectral types, the measured RVs with errors for all B and Be stars as well as BSGs, and possible comments.

All Tables

All Figures

	Fig. 1. Observed close binary fraction for MS stars as a function of the adopted RV-variability threshold ΔRVmin. For each ΔRVmin we compute the observed close binary fraction of the entire sample as well as the corresponding binomial error. The vertical dotted line gives the adopted threshold value ΔRVmin = 20 km s−1 indicated by the red line.
In the text

Fig. 2. Binary detection probability simulations following our observational setup. Top panel: detection probability as a function of the orbital period P (computed assuming a flat mass-ratio distribution), and bottom panel: this as a function of mass ratio q (computed assuming a flat logP distribution). In both panels, the black line shows our first estimate of the binary detection probability taking the SB1 bias into account, and the orange line shows the detection probability, which also includes the reduced sensitivity for detecting near equal-mass binaries. The colored arrows in the bottom panel correspond to the mass ratio for which we compute the SB2 bias (see Fig. 3).

In the text

	Fig. 3. Measured primary RVs for simulated spectra of one single B2 V star, and SB2 binary systems with a B2 V primary and six different mass ratios as indicated in the legend. A rotational velocity of v sini = 200 km s−1 is assumed for all stars. The horizontal axis gives the assumed RV difference between the two components. The dashed line indicates the true velocity of the primary star.
In the text

Fig. 4. Hubble Space Telescope CMD (Milone et al. 2018) overlaid with binarity information. B-type stars are marked in yellow, Be stars in green, BSGs in blue, and RSGs in red. Open small circles indicate that no RVs could be measured, filled circles are stars with no significant RV variability, diamonds indicate stars classified as binary systems based on their RVs, and triangles indicate SB2 candidates.

In the text

Fig. 5. Distribution of spectral types and observed close binary fraction for B and Be stars. Top panel: distribution of spectral types for all MS stars (dim gray) and stars for which RVs were measured (dark gray). Bottom panel: observed spectroscopic binary fraction as a function of spectral type. The red curves indicate the binary fraction inferred from the RV measurements, and the orange curve combines the RV measurements with the SB2s found by visual inspection.

In the text

Fig. 6. Distribution of F814W magnitudes and observed close binary fraction for B and Be stars. Top panel: distribution of F814W magnitudes of all stars for which spectra were extracted (dim gray), and stars for which RVs were measured (dark gray). Bottom panel: observed spectroscopic binary fraction as a function of F814W magnitude. The orange curve is based on RV measurements, and the red curve is based on RV measurements as well as SB2s. The bias-corrected close binary fraction is overplotted in dark red.

In the text

Fig. 7. Binarity of different stellar groups, based on their position in the CMD. MS stars are marked in orange, stars in the Be-star region in dark green, stars around the cluster TO in red, and blue stragglers in blue. Additional HST sources for which no spectra were extracted are marked with gray dots, and stars without RV measurements are indicated with open circles. Colored filled circles indicate that the star is single, colored diamonds indicate that the star is classified as a binary based on RV variability, and colored triangles indicate SB2s. The individual systems of interest discussed in Sect. 5.3 are labeled.

In the text

	Fig. A.1. Binary detection probability as a function of period, mass of the primary M1, mass ratio q and eccentricity e, computed for the SB1 bias.
In the text

	Fig. A.2. Same as Fig. A.1, but computed taking both the SB1 and the SB2 bias into account.
In the text

Fig. B.1. Binary detection probability as a function of period (top panel) and mass ratio (bottom panel). The different dotted orange curves correspond to different magnitude bins, in which the detection probability is highest for the brightest stars. The continuous orange line shows the overall detection probability for the full simulation, and the continuous black line shows the detection probability for the SB1 bias alone (see Sect. 4.1).

In the text