© The Authors 2023

1. Introduction

Apart from the initial mass, the initial rotation rate is considered to be one of the key properties with a major influence on the evolution and end points of (massive) stars (e.g., Maeder & Meynet 2000; Brott et al. 2011; Ekström et al. 2012; Woosley & Heger 2006). Rotation can drive internal mixing processes that can have several effects on the star: new fuel is mixed into the core, the main-sequence (MS) lifetime is increased, the luminosity is significantly altered, and nuclear processed material might be brought to the surface (see e.g., Maeder & Meynet 2000; Howarth & Smith 2001). Several mixing mechanisms are thought to be at play in massive stars, including semi-convection, convection, thermohaline mixing, and mixing induced by gravity waves (e.g., Rogers et al. 2013) or rotation (for example by shear or meridional circulations, see e.g., Maeder & Meynet 2000). One-dimensional evolutionary models show that the more rapidly a star is rotating at birth, the stronger the impact of rotation-driven internal mixing (Brott et al. 2011). The effectiveness and importance of rotational mixing have, however, been questioned and alternative mechanisms to change the surface spin and surface abundances of stars have been proposed (e.g., Hunter et al. 2007; Almeida et al. 2015; Aerts et al. 2019).

High surface rotation can result both from single- and binary-star evolution: through angular momentum redistribution from the core to the envelope during the MS (see e.g., Ekström et al. 2008; Hastings et al. 2020, if the initial rotation is sufficiently high), or through accretion of angular momentum during mass transfer (e.g., de Mink et al. 2013) or by tides (de Mink et al. 2009). Recent binary population synthesis predictions by Wang et al. (2020) show that, combined with the position in the Hertzsprung-Russell diagram (HRD), the rotation rate of stars is one of the main tracers of a possible binary history (see also Wang et al. 2022a). This is most obvious for mass gainers through case-B mass transfer (Kippenhahn & Weigert 1967), which are expected to be spun up to high rotation rates. Those might appear as classical Be stars observationally (e.g., Pols et al. 1991; Hastings et al. 2021), which are rapidly rotating B-type stars surrounded by a circumstellar decretion disk that leads to strong emission lines in their spectra (e.g., Rivinius et al. 2013). In contrast, merger products are thought to be slowly rotating as they are expected to spin down quickly after an initial puffed-up stage (Schneider et al. 2020, and references therein).

Constraining the rotation rates of large, homogeneous samples of massive stars as a function of their metallicity, evolutionary stage, and physical parameters greatly helps to identify the mechanism and evolutionary scenarios that are responsible for observed distributions of surface rotation velocities. Milky Way studies indicate that cluster stars rotate, on average, faster than field stars (Slettebak 1949; Abt et al. 2002; Strom et al. 2005; Huang & Gies 2006; Daflon et al. 2007; Garmany et al. 2015). These differences are mainly explained in two different ways: on the one hand, they are attributed to the star formation process, where stars in denser environments are thought to be formed with higher initial rotation rates than stars in less dense environments (e.g., Strom et al. 2005; Wolff et al. 2007). On the other hand, they are explained as an evolutionary effect where field stars are on average more evolved, and hence rotate slower as their radius has increased toward the end of the MS (e.g., Huang & Gies 2006, 2008). Focusing on the fastest rotators, Huang et al. (2010) found that they are likely either formed by evolutionary spin-up of the most massive stars, or by binary mass transfer. Bragança et al. (2012) reported that the rotational velocities of runaway B-type stars in the field resemble those of cluster B-type stars, suggesting that they were born in clusters and later ejected from these denser environments.

Comparing rotational velocity distributions of massive OB stars in the Milky Way to those of the Large Magellanic Cloud (LMC) association N 206 (Ramachandran et al. 2018) and the Small Magellanic Cloud (SMC) wing, Ramachandran et al. (2019) show that stars in the SMC have, on average, higher rotational velocities. This is attributed to two effects: firstly, lower-metallicity stars are more compact and, at a given mass, thus generally have higher rotational velocities (e.g., Ekström et al. 2008). Secondly, lower-metallicity stars have weaker winds (if they exhibit winds in the first place, Mokiem et al. 2007) such that they lose less angular momentum.

In the LMC, the VLT FLAMES Tarantula Survey (VFTS; Evans et al. 2011) homogeneously measured the rotational velocity distributions of OB stars in the young 30 Dor region, a massive star-forming complex. Ramírez-Agudelo et al. (2013) studied the rotation rates of about 200 apparently single O-type stars in 30 Dor. They found two components in the projected rotational velocity (v sin i) distribution. A majority of stars rotate at around 80 km s⁻¹, while there is a high-velocity tail of fast rotators. This tail is qualitatively similar to expectations from binary interaction products (de Mink et al. 2013). Focusing on the primaries of O-type stars in binary systems in the same region, Ramírez-Agudelo et al. (2015) found that primaries in short-period systems are faster rotators than single O-type stars likely as a result of tidal synchronization (see also Mahy et al. 2020). Furthermore, they found a clear lack of very rapidly rotating binary O-type stars, in line with expectations that a prior interaction is needed to produce such rapid rotation.

Dufton et al. (2013) show that the 30 Dor B-star v sin i distribution is bi-modal, with one peak below 100 km s⁻¹ and a high-velocity component above 250 km s⁻¹. Given that this study was performed in a cluster environment and no spatial or radial differences were found in the rotation rates, it seems unlikely that the bi-modality is due to observations of cluster and field stars, or stars formed in different star formation episodes. Hunter et al. (2008) found the rotation rates for LMC B-type stars to be similar to the ones reported for Galactic B-type stars.

Focusing on the ∼2-Myr old cluster NGC 346 in the SMC, Dufton et al. (2019) investigated the binary properties and rotational velocities of almost 250 O- and early B-type stars. While the authors found indications that the stars in NGC 346 are on average more rapidly rotating than Galactic early-type stars, the differences in the rotational velocity distributions are not statistically significant. The possible higher than average rotational velocity would be in line with the larger number of rapidly rotating Oe and Be stars detected in SMC clusters (Martayan et al. 2007a,b).

Photometric studies of slightly older clusters, that is with ages between ∼30 Myr and 1 Gyr, reported the presence of a split MS in the cluster color-magnitude diagrams (CMDs). One of the proposed explanations attributed this to different rotational velocities of stars in the two sequences (e.g., Bastian & de Mink 2009), which was confirmed by spectroscopic observations (see e.g., Marino et al. 2018; Kamann et al. 2021, 2023, and Sect. 4.2). The origin of the different rotational velocities of the populations was, among other theories, explained as an imprint of previous binary interactions in those clusters (see e.g., Wang et al. 2020, 2022a,b; Sun et al. 2021).

In this context, the intermediate-age, massive SMC cluster NGC 330 offers the possibility to obtain the rotational velocity distribution of stars in a low metallicity environment. Previously, NGC 330 was estimated to be 20 Myr old (Keller et al. 1999), but more recent works found ages of around 35−40 Myr (Bodensteiner et al. 2020; Eldridge et al. 2020; Patrick et al. 2020). This puts it in an age range where recent theoretical computations predict a large number of binary interaction products to exist.

This is the third paper in a series. In Bodensteiner et al. (2020, hereafter Paper I), we describe our multi-epoch observations on NGC 330 obtained with the Multi Unit Spectroscopic Explorer (MUSE) and characterize the massive star content. In Bodensteiner et al. (2021, hereafter Paper II), we measure radial velocities (RVs) and investigated the current binary properties of the massive star population. For this work, we used MUSE spectra and HST photometry to quantify their physical properties. In Sect. 2 we introduce a newly implemented grid-based method to derive stellar parameters. We describe our resulting stellar parameters, in particular the rotational velocities, in Sect. 3. In Sect. 4 we compare our findings to previous works of clusters in both the LMC and SMC, and give a summary in Sect. 5.

2. Determination of atmospheric and physical parameters

Here, we describe the estimation of stellar parameters, namely the effective temperature T_eff, surface gravity log g, projected rotational velocity v sin i (where v is the stellar rotational velocity and i is the inclination angle), radius R, and bolometric luminosity L. As a majority of the stars in the sample are hotter than 15 000 K, non-LTE effects may be important in their atmospheres (Kudritzki 1979; Kilian 1994; Nieva & Przybilla 2007). We develop a grid-based algorithm based on the non-LTE TLUSTY BSTAR2006 grid (Hubeny & Lanz 1995; Lanz & Hubeny 2007). The grid provides synthetic spectra for T_eff between 15 000 and 30 000 K with steps of 1000 K, and log g from 3.0 to 4.75 dex with steps of 0.25 dex, where g is expressed in cgs units. We use the grid computed for SMC metallicity (Z = 0.2 Z_⊙), assuming a fixed microturbulent velocity of 2 km s⁻¹. The synthetic spectra cover the UV to optical wavelength range (900–10 000 Å).

Our newly developed fitting routine¹ provides two modules to derive stellar parameters: (1) the optical co-added MUSE spectra are compared with rotationally broadened TLUSTY model spectra degraded to the appropriate spectral resolution and wavelength binning of MUSE (see Sect. 2.2). (2) the available high-quality HST photometry from Milone et al. (2018) is compared to TLUSTY model SEDs to provide additional constraints on the radius of the star (see Sect. 2.3), which allows an estimate of the bolometric luminosity.

2.1. Combination of MUSE spectra

The MUSE observations, the reduction of the data including sky subtraction, the extraction of spectra and their normalization is explained in detail in Paper I. We extracted four to six spectra (covering the wavelength range between 4600 and 9300 Å at a resolving power between 1700 and 3700) for 324 B and Be stars, six blue supergiants (BSGs) and ten red supergiants (RSGs) and measured RVs at each observed epoch for 282 of the B and Be stars and all of the BSGs (Paper II).

Knowledge about the RV variability of the stars allows us to shift the spectra to rest wavelength and stack them, weighted by the signal-to-noise ratio (S/N) at 6100 Å. In the following, we focus on the 191 B and 83 Be stars. We omit eight stars with composite, variable spectra indicative of double-lined spectroscopic binary systems (SB2, found in Paper II). Among the 274 analyzed stars, 26 are previously detected single-lined spectroscopic binaries (SB1s), where only one binary component can be identified in the spectra. Combining the spectra significantly increases their S/N with typical values between 200 for the fainter and over 400 for the brighter B and Be stars. An example of the individual spectra of star #654 (with individual S/Ns between approximately 125 and 290 at around 6100 Å) and the combined spectrum with a S/N of over 450 is shown in Fig. 1.

Fig. 1. Example of the co-addition of spectra for the apparently single star #654. Individual epochs are shown in color (as indicated in the legend) while the combined spectrum is shown in black. The legend also indicates the measured S/N at 6100 Å. Main spectral lines are marked.

2.2. Spectroscopic fit

The observed MUSE spectra allow us to constrain T_eff, log g and v sin i of a star, which we do based on a grid search. For each model in the two-dimensional TLUSTY grid (T_eff − log g), we apply the instrumental broadening corresponding to the MUSE (wavelength-dependent) resolving power, rebin the models to the MUSE wavelength step of 1.25 Å, and add rotational broadening using the python package PYASTRONOMY (Czesla et al. 2019). We here consider rotational velocities ranging from 0 to 500 km s⁻¹ with steps of 20 km s⁻¹, which is chosen to be lower than the typical accuracy of the v sin i-estimates as discussed below. We fit T_eff, log g, and v sin i simultaneously as the former two cannot be fitted without taking rotational broadening into account, and vice versa. Additionally, the uncertainties can be correlated.

We then compare the observed spectrum to every model in the grid by computing the spectroscopic ${\chi }_{\text{spec}}^2$ for each model as

$$\begin{array}{c}\hfill {\chi }_{\mathrm{spec}}^2=\frac{1}{S/N}·{\displaystyle \sum _{i=1}^n{({\mathrm{obs}}_i-{\mathrm{model}}_i)}^2,}\end{array}$$(1)

where obs_i and model_i are the observed spectrum and the model at a given wavelength with index i, respectively, n is the total number of spectral pixels considered. We use the measured S/N rather than individual errors on each pixel of the spectrum as these errors are typically overestimated by the spectral extraction via the fitting of the point-spread function (PSF), especially for fainter stars.

We define six diagnostic regions in the MUSE spectral range (also shown in Fig. 2): (1) the blue part of the spectrum between 4600 and 5060 Å, which contains several diagnostic lines such as Hβ, four He I lines at λλ4713.15, 4921.93, 5015.68, 5047.74 Å, and some metal lines of Si II, N II, C II and O II (which are in most cases too narrow and weak to be detected in MUSE), (2) the region around the He II line at 5411 Å, (3) the region around Hα, between 6500 and 6640 Å, (4) the He I line at 6678.15 Å, (5) the He I line at 7065.19 Å, and (6) the Paschen region between 8500 and 8930 Å.

Fig. 2. Spectroscopic fit for star #654. Each panel, one for each of the six diagnostic spectral regions, shows the combined spectrum (black) and the best-fitting TLUSTY model (red). Grayed-out regions are only shown for comparison but were not included in the fit.

Often-used diagnostic lines for (early) B-type stars, such as the ratios between He I and Mg II, or Si II and Si III, cannot be used as the lines are either too weak or outside the wavelength region covered by MUSE. The He I line at λ5875.62 Å is also not available due to a Na I-filter required when using MUSE with the laser-supported adaptive optics system (see Paper I).

The regions containing the Hα and Paschen lines were particularly affected by nebular emission, telluric absorption, and sky subtractions. Additionally the Hα line in the observed spectra turns out to be systematically stronger than in the models, which might be caused by an over-subtraction of the sky. The Paschen region further suffers from identifying the true continuum. We therefore refrain from including these regions in the ${\chi }_{\text{spec}}^2$-computation but use them as a consistency check when visually inspecting the fits. After visually inspecting every fit, we further exclude other regions for some stars when they appear noisy, contaminated, or badly normalized. The ${\chi }_{\text{spec}}^2$ was computed in all remaining regions.

The estimation of rotational velocities based on the MUSE observations is mainly hampered by the relatively low resolving power, which varies between R = 1700 − 3700 from the blue to the red part of the spectrum. This corresponds to a formal velocity resolution of about 150 to 80 km s⁻¹, respectively. Previous studies based on MUSE data have, however, shown that a determination of v sin i with a 1σ accuracy of 30 km s⁻¹ is possible for cooler stars that provide a multitude of spectral lines (e.g., Kamann et al. 2018). Here we show that it is also possible for hotter stars when using several spectra lines for our MUSE targets that have significant rotational broadening.

Given the availability of the TLUSTY models we have to adopt a microturbulent velocity v_micro = 2 km s⁻¹. As typical values of v_micro are on the order of a few km s⁻¹ for B-type MS stars, this is much lower than the typical v sin i and thus a valid assumption. Furthermore we assume that the macroturbulent velocity v_macro is negligible. The macroturbulent velocity for B-type stars is below 100 km s⁻¹ in most cases (see e.g., Simón-Díaz & Herrero 2014). While this is significant, it is lower than the typical rotational velocities and will be below or similar to the resolution capabilities of MUSE for most stars. The derived rotational velocities are thus upper limits and might be, on average, slightly overestimated. Our fitting approach also does not take into account the deformation of very rapidly rotating stars, which can have a strong impact on the derived stellar parameters (Abdul-Masih 2023). The assumptions are further discussed in Sect. 2.5.

For the Be stars, we more frequently exclude certain regions from the fit as their spectra often show strong emission, mainly in the Balmer lines, but also sometimes in He I and other spectral lines, which cannot be used to estimate the photospheric parameters of the stars. For stars that show only moderate contamination by emission lines in some of the diagnostic lines, we follow the same procedure as for B-type stars. We note that their derived parameters in general have larger uncertainties. In particular, when excluding the often contaminated Hβ line, the surface gravity is not as well constrained (see Sect. 2.4). For thirteen Be stars with strong emission lines, it is not possible to derive parameters based on the spectra and we estimate only rotational velocities as described in Sect. 2.6.

2.3. Photometric fit

In addition to the constraints from the MUSE spectra, we make use of the high-quality HST magnitudes available for almost all stars (a handful of stars fall in the CCD gap or are too close to a saturated star). HST observations were obtained in three broad-band filters covering the SED from the near-UV to the red edge of the optical range, namely F225W, F336W, and F814W (centered at 2250 Å, 3360 Å, and 8140 Å, respectively, Milone et al. 2018). The catalog also includes errors for F225W and F336W magnitudes, computed as random mean scatter (RMS) of the several independent photometric observations that are available. No errors are given for the F814W filter as there was only one epoch of observations. After conversion to fluxes, and correcting the errors for small sample statistics, we estimate the F814W error assuming the relative flux error is the same as the relative error on the F225W flux measurement.

Similarly to the spectroscopic fit, we again find the best-fitting model based on a grid search. For each model in the TLUSTY grid of T_eff and log g, which provides the Eddington flux at the stellar surface, we compute the observed flux by scaling it with the distance of the SMC of 60 kpc (e.g., Deb & Singh 2010) and apply interstellar extinction using the python package DUST_EXTINCTION assuming E(B  −  V) = 0.08 (Graczyk et al. 2014), a constant R_V = 3.1, and the reddening law of Fitzpatrick et al. (2019). The stellar radius is one of the fitting parameters and we further scale the model flux by varying it between 2 and 20 R_⊙ in steps of 0.1 R_⊙. For each model in this now three-dimensional grid, we then compute the flux of the scaled and reddened model by convolving it with the corresponding HST filter².

We refrain from letting the extinction and distance vary as this would lead to a higher degeneracy in the fit. Given the generally low extinction, the choice of E(B  −  V) has a small impact on the overall results, especially when combined with the spectroscopic fit. The distance to the SMC is well-known, and small-scale variations in the distance (for example the position of a star within the cluster), do not play a role given the overall large distance to the SMC. We do, however, note that only three photometric datapoints are available and fixing the distance and extinction might lead to systematic errors.

We then compute a photometric ${\chi }_{\text{phot}}^2$ for each model in the 3D-grid by comparing the three observed fluxes with the computed model fluxes according to:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{phot}}^2={\displaystyle \sum _{i=1}^n\frac{{({\mathrm{obs}}_i-{\mathrm{model}}_i)}^2}{{\sigma }_i^2},}\end{array}$$(2)

where obs_i are the measured fluxes in the before-mentioned three HST filters, model_i the filter-convolved model fluxes in the same filters, and σ_i corresponds to the individual errors of the three flux measurements (i.e., n = 3 for most stars).

For a handful of stars, no photometric fit is possible as they have no HST magnitude available. Similarly to the spectroscopic fit, the photometric fit of Be stars is complicated by the circumstellar disk, in particular by excess flux in the near-IR (see e.g., Rivinius et al. 2003). While the F225W and F336W magnitudes should be relatively clean, the F814W magnitude is potentially affected by disk emission. As photometry and spectroscopy are not taken simultaneously, it is possible that the contamination varies between the two because of the variability of the disk. Given that we combine the photometric fit with the spectroscopic one, we generally find an acceptable fit for the Be stars, but with larger uncertainties than for B-type stars. We refrain from neglecting the F814W-magnitude for Be stars as we would be left with only two observed data points (i.e., F225W and F336W).

2.4. Combination of the two fits

We combine the two separate fits by adding up the computed χ² for each model with the same T_eff and log g, but different v sin i and radius:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{global}}^2={\chi }_{\mathrm{spec}}^2+{\chi }_{\mathrm{phot}}^2.\end{array}$$(3)

By doing this without a further normalization, more weight is given to the spectroscopic fit in comparison to the photometric fit when more spectral lines (and thus more data points) are available. For spectra with a lower S/N or highly contaminated spectral where only few lines are used, the photometric fit is weighted more strongly. For most stars, both individual fits agreed within the errors, with the exception of some of the SB1 systems. A handful of additional stars were flagged based on such a discrepancy as it might be indicative of undetected binaries.

The best-fitting model in this four-dimensional grid in terms of T_eff, log g, v sin i and stellar radius is found by locating the minimum in the obtained ${\chi }_{\text{global}}^2$ distribution for each parameter. We rescale the obtained χ² distribution such that the model with the minimum χ² has a ${\chi }_{\text{global}}^2$ equal to the degrees of freedom in the computation. We furthermore give 2-σ errors based on estimating the 95%-confidence interval on the parameter range after rescaling the χ². These are formal errors and do not take into account systematics, which could arise, for example, from fixing the microturbulent velocity. The limited availability of diagnostic lines in the observed spectra and the limited available model grid can affect the final parameters even beyond the given uncertainty ranges. Although the formal fitting errors in the derived stellar properties listed are quite small, the actual errors in the physical properties may be considerably larger because of systematic effects.

Typical uncertainties on T_eff, log g, and v sin i are on the order of 1000 K, 0.25 dex, and 70 km s⁻¹, respectively. While the v sin i uncertainty is in principle higher for higher rotational velocities, the dominating factor is the S/N in the spectra. The errors derived for Be stars are, on average, higher than for B-type stars as fewer spectral lines can be used on average, and the F814W-magnitude might be impacted by the disk. In particular, this affects the determination of the surface gravity as for Be stars the Balmer lines (mainly also Hβ, which is log g-sensitive and used as diagnostic line for B-type stars) is contaminated by emission. As mentioned above, when using several spectral lines, the observed rotational velocity has a better accuracy than the formal velocity resolution of MUSE.

For ∼30 stars, the derived effective temperature reached the edge of the TLUSTY grid, when the associated errors were included. For those stars, we computed additional TLUSTY models extending to T_eff = 10 000 K (L. Mahy, priv. comm., covering the same parameter range in log g and v sin i as the nominal models) and repeated the analysis. If the error bars still hit the edge of the grid, we only provide upper or lower limits. This occurred predominantly for Be stars with strong emission, for which the stellar parameters are difficult to constrain. Due to the infilling of the Balmer lines, in particular a diagnostic for the surface gravity is lacking, which propagates to the determination of the effective temperature. Best-fit parameters with associated uncertainties are given in Table A.1. Diagnostic plots for each star (similar to the examples given in Figs. 2–4) are available on Zenodo.

2.5. A test case: Star #654

We provide a proof of concept for one star, #654, and use it to test the assumptions on v_macro, distance and extinction. #654 is a B-type star with a F814W magnitude of 16.3 mag and a spectral type of B2 (Paper I) located in the turnoff (TO) region of NGC 330. #654, also referred to as NGC 330−095, was classified as B3 III by Evans et al. (2006) and its stellar parameters were further investigated by Hunter et al. (2008). Given that we found no significant RV variability (Paper II) we presume it is a single star (which agrees with the classification in Hunter et al. 2008).

The spectroscopic fit of the observed MUSE spectrum, shown in Fig. 2, is sensitive to the effective temperature, the surface gravity and the projected rotational velocity. The photometric fit, shown in Fig. 3, constrains mainly the effective temperature and the radius. The global χ² distributions derived of the combined fit can be seen in Fig. 4.

Fig. 3. Photometric fit for star #654. The top panel shows a comparison between the observed HST fluxes (black crosses) and the flux (colored circles) computed from the best-fit model (gray line) by convolution with the HST filters (colored lines). The bottom panel gives residuals.

We find that #654 has an effective temperature of $21{000}_{-2300}^{+300}$ K a surface gravity of $3.8_{-0.3}^{+0.1}$ dex, and a projected rotational velocity below 40 km s⁻¹. Given the estimated radius of $6.8_{-0.6}^{+1.3}{R}_{\odot }$, the luminosity of the star is log(L/L_⊙) = $3.9_{-0.3}^{+0.2}$. The effective temperature is in good agreement with the one expected for a spectral type B2. The surface gravity agrees with the position of #654 at the cluster TO, implying that it is close to the terminal-age MS (TAMS), as well as with the B3 III classification of Evans et al. (2006). The obtained radius is slightly higher than the typical radius of a B2 star on the MS, but probably in agreement with the radius of such a star at the TAMS. The rotational velocity of the star can be constrained below the formal resolution limit of MUSE because of the high S/N of the spectra and the large number of spectral lines used in the fit.

Hunter et al. (2008) assumed the effective temperature T_eff = 18 450 K from the spectral type, and based on this inferred log g = 3.45, and v sin i = 20 km s⁻¹ from a comparison to TLUSTY models. Assuming E(B  −  V) = 0.06, R_V = 2.72, the reddening law of Bouchet et al. (1985), and bolometric corrections from Balona (1994), they derived log(L/L_⊙) = 3.67. While the rotational velocity is in good agreement, the known degeneracy between T_eff and log g (also illustrated in Fig. 4) explains the lower derived log g by Hunter et al. (2008, i.e., decreasing T_eff by 1000 K decreases the log g-estimate from the Balmer lines by approximately 0.1 dex, which brings the estimate in good agreement). Despite the different assumptions on reddening and the different techniques used, the luminosities are comparable.

Fig. 4. Combined fit for star #654. The diagonal gives the χ2-distribution as a function of each parameter and the panels below show 2-dimensional χ2-maps. The red line marks the 95% confidence level.

The only other star in common with the samples of Evans et al. (2006) and Hunter et al. (2008) is #57 or NGC 330−036. It was classified as B2 II star by Evans et al. (2006) and Hunter et al. (2008) report T_eff = 21 200 K (again inferred from the spectral type), log g = 3.25 dex, v sin i = 39 km s⁻¹, and log(L/L_⊙) = 4.33. In our approach, we find comparable parameters, namely T_eff = $24{000}_{-200}^{+700}$ K, log g = $3.5_{-0.1}^{+0.1}$ dex, v sin i < 50 km s⁻¹, R = $10.3_{-0.1}^{+0.1}$ R_⊙, and log(L/L_⊙) = $4.5_{-0.1}^{+0.1}$. Again, the lower assumed effective temperature in Hunter et al. (2008) will lead to the lower log g derived. Overall, we thus find a good agreement between the studies, despite the lower resolving power of the MUSE spectra. To further test our method, we fit the published FLAMES spectra of #654 and #57 with a similar approach than the MUSE spectra, retrieving stellar parameters that agree within the errors.

We further use #654 to test our fitting assumptions. As described in Sects. 2.2 and 2.3, we assume a constant E(B  −  V) = 0.08, A_V = 3.1 and neglect v_macro. When changing A_V to a lower value, for example A_V = 2.74 as reported for the SMC (Gordon et al. 2003), the T_eff and R derived from the photometric fit remain the same within the error bars. Similarly, changing E(B  −  V) by up to ±20% does not strongly impact the results. A larger change in the extinction, for example caused by a local overdensity of material around specific sources, would lead to a significantly higher T_eff (i.e., when assuming twice the extinction, E(B  −  V) = 0.16, the photometric fit implies T_eff = 28 000 K). Such a high temperature is excluded by the spectroscopic fit and the combined fit, which for this star is weighted toward the spectroscopic fit, does not change significantly.

Finally, we use IACOB-BROAD (Simón-Díaz et al. 2017), which determines the rotational velocity of a star based on a combined Fourier transform and goodness-of-fit method taking into account a nonzero v_macro. We find $vsini={36}_{-11}^{+46}$ km s⁻¹ and $v_{\mathrm{macro}}={81}_{-81}^{+23}$ km s⁻¹. This again shows that v sin i can be constrained better than its formal accuracy. Further, the derived v sin i agrees with the one estimated in our approach, and the derived v_macro is consistent with our assumption of v_macro = 0 km s⁻¹.

The test case #654 demonstrates the working principle of the fitting approach used here. As #654 is a particularly slowly rotating star, we show the fitting of a second test case, a rapidly rotating Be star, in Appendix B.

2.6. Estimation of v sin i based on FWHM measurements

As described above, several Be stars (13 out of the 83) could not be fitted with the TLUSTY approach as some of their spectral lines are dominated by emission. We thus resort to another commonly used method for estimating their projected surface rotation velocity. It is based on measuring the full-width at half maximum (FWHM) of spectral lines (see e.g., Slettebak et al. 1975; Abt et al. 2002; Strom et al. 2005; Bragança et al. 2012; Garmany et al. 2015), which to the first order is linearly related to the projected rotational velocity. Lines that are typically used for this purpose in B-type stars are lines in the blue part of the spectrum (mainly He I at λλ4388 and 4713, C II at λ4267, or O II at λ4366; see e.g., Abt et al. 2002). These He I lines were also used to estimate the rotational velocity of Be stars (Steele et al. 1999). As described in Sect. 2.2, the FWHM-method neglects the effect of additional broadening by micro- and macroturbulence.

In Paper II, we fit Gaussian profiles to spectral lines (He Iλ4922, among others) in order to measure RVs for the B and Be stars in the sample. From all stars with both FWHM estimates and v sin i measurements from the fitting described in Sect. 2.2, we derive a FWHM-v sin i relation which we then apply to the thirteen Be stars that we could not fit with the TLUSTY approach. We here use the FWHM measurements of the He I line at λ4922 Å and exclude stars from the fit whose error bars hit the boundary of our grid, that is they are consistent with zero or reach velocities higher than 500 km s⁻¹.

The measured FWHM [Å] and v sin i [km s⁻¹], shown in Fig. 5, follow a roughly linear relation. Fitting it with a first-order polynomial, weighted by the measured errors in FWHM and v sin i, we find:

$$\begin{array}{c}\hfill vsini[{\mathrm{km}\mathrm{s}}^{-1}]=-183.3[{\mathrm{km}\mathrm{s}}^{-1}]+67.1\left[\frac{{\mathrm{km}\mathrm{s}}^{-1}}{\mathrm{\AA }}\right]·\mathit{FWHM}[\mathrm{\AA }],\end{array}$$(4)

Fig. 5. Relation between the FWHM of the He I line at 4922 Å with the v sin i estimated from the TLUSTY fit (purple line, see Eq. (4)). The purple points indicate the measured values for B-type stars.

with an RMS residual of 43 km s⁻¹. As a consistency check, we repeat the fit for stars with fractional v sin i errors from the TLUSTY fit < 0.2. The difference in the determined relation is negligible (as expected given the weighting by the errors). Ramírez-Agudelo et al. (2015) derived a similar relation for O-type stars using the He I line at λ4922 Å. Their nonlinear relation deviates from our linear relation for lower rotational velocities (below ∼150 km s⁻¹). Otherwise, the obtained rotational velocities using the relation by Ramírez-Agudelo et al. (2015) agree with ours within the RMS.

2.7. Deriving the intrinsic rotational velocity distribution

The sample sizes of the B- and Be-star in NGC 330 are large enough to estimate their current rotational velocity distribution. Assuming that their rotation axes are randomly distributed, we can infer the probability density function, p(v_e), using the iterative procedure of Lucy (1974) as implemented by Dufton et al. (2013). This procedure has been used previously for B-star samples in the VFTS (Evans et al. 2011) of 30 Dor (Dufton et al. 2013, 2022), and for NGC 346 (Dufton et al. 2019, where further details can be found). Here, convolutions were undertaken separately for both the B- and Be-type samples in NGC 330, as well as for the total sample. In all cases, the stars classified as apparently single and as SB1 candidates were included. The latter may not have a random distribution of orbital (and hence by implication rotational) axes of inclination, as discussed by Ramírez-Agudelo et al. (2015). However, given the low number of SB1s, and combined with the likely presence of undetected binaries in the apparently single star sample, the overall distribution of rotational axes of the entire sample should be close to random.

In addition, we have also re-derived rotational distribution functions for equivalent samples in the VFTS and NGC 346 surveys. The v sin i-estimates for both the single and SB1 subsamples were combined using the following sources: Dufton et al. (2013, 2022) for the VFTS apparently single stars, Garland et al. (2017), Dufton et al. (2022) for the VFTS SB1 primaries, and Dufton et al. (2019) for the NGC 346 targets. No correction has been made to the v sin i-estimates for gravity darkening (see, for example, Rivinius et al. 2013; Zorec et al. 2016), in order to maintain consistency with our NGC 330 estimates.

Care should be taken when interpreting the small-scale variation in the rotational distributions presented here. Even for the largest datasets, the implied stochastic uncertainties are significant and could lead to a spurious small-scale structure. Nevertheless, we compare the intrinsic rotational velocity distribution to those derived for B and Be stars in 30 Dor and NGC 346 (Dufton et al. 2013, 2019, 2022) in Sect. 4.2.

3. Stellar parameters of the massive star population

We attempt a combined spectroscopic and photometric fit for 274 B and Be stars. However, the spectra were too noisy or contaminated, or photometric measurements were missing for 23 stars such that no stellar parameters could be determined. For a subset of the 23, namely thirteen stars, the projected rotational velocity could be derived based on the FWHM-v sin i-correlation described above. In total, this provides us with stellar parameters of 251 stars and v sin i-measurements of additional thirteen stars (264 stars in total), including presumably single stars and SB1s.

3.1. Temperatures, surface gravities, radii and luminosities

Using the HST magnitudes of Milone et al. (2018), we construct color-magnitude diagrams (CMDs) to investigate the distribution of the stellar parameters for the combined fit. Figure 6 shows that the temperature increases with the color (m_F336W − m_F814W) of the stars, as expected, but has a significant scatter. This is especially the case for the SB1s (see below) and for the stars above the cluster TO. The radius of the stars increases for the brighter stars, as well as for the more evolved stars, in line with the expectation that stars increase their radius at the end of the MS. The measured surface gravities generally also reflect this: they decrease for the more massive and more evolved stars. They are not corrected for the effect of gravity darkening. The estimated values show again a larger scatter, which is most likely due to the significant uncertainties in the log g estimates.

Fig. 6. CMDs of NGC 330 color-coded by the best-fitting physical parameters. Circles indicate presumably single stars while triangles indicate stars classified as SB1s. The isochrones are single-star nonrotating Padova isochrones of 35 (dashed line) and 40 Myr (dotted line).

In Fig. 6 we also indicate previously detected SB1s. While their derived radii and luminosities are in agreement with expectations, their temperatures and surface gravities show a larger scatter. The surface gravities of the SB1s do not show a trend in comparison to the single stars, but the effective temperatures of the SB1s seem, on average, slightly lower than the ones of single stars at the same locations in the CMD. The latter is consistent with spectral contamination from a presumably cooler unseen secondary. The greater scatter in the gravity estimates would also be consistent with such contamination by an unseen companion in, for example, the gravity-sensitive wings of Balmer lines.

One particularly hot star stands out in the CMD at m_F814W = 15.9. This star, #420, is one of the two stars showing He II absorption in its spectrum. The spectrum is furthermore dominated by emission lines, hampering the quality of the fit. As shown in Fig. 6, the log g estimated for #420 is low given its position in the CMD. This is probably connected to the fact that the estimated effective temperature of 30 000 K is at the edge of the used TLUSTY B-star grid, and the star is probably hotter (which would imply a higher log g).

3.2. Projected rotational velocities of B and Be stars

As mentioned before, information about the projected rotational velocity is available for 264 stars. Figure 7 shows that the v sin i of most B-type stars is around 100−250 km s⁻¹ with a tail to high velocities up to 500 km s⁻¹. The projected rotational velocities of Be stars show a distribution that is shifted to higher velocities: the majority of the Be stars are rotating with inferred projected velocities between 200−400 km s⁻¹.

Fig. 7. Stacked histogram of the projected rotational velocities of presumably single stars and SB1s. The top panel shows the B-type stars (single stars in yellow, SB1s in orange) and the bottom panel shows Be stars (single stars in light green, SB1s in dark green). The gray shaded region again corresponds to the formal MUSE sensitivity limit while the indicated errors are counting errors.

The median projected rotational velocity of the B-type stars is around 170 km s⁻¹ while it is at ∼250 km s⁻¹ for the Be stars. We perform a Kuiper test and reject the null-hypothesis that both are drawn from the same parent distribution at the 10⁻³-significance level. The observed projected rotational velocity distributions of B and Be stars are therefore statistically different from each other, in agreement with findings from previous works (see e.g., Hunter et al. 2008; Dufton et al. 2022).

In Fig. 7 we also compare the projected rotational velocities of B and Be stars classified as presumably single or SB1s. Given the low observed binary fraction in NGC 330, there are only a handful of SB1s in the sample. Despite the low number statistics (or because of it), we find no obvious differences in the projected rotational velocity distributions of SB1s relative to presumably single stars. There is one rapidly rotating Be star classified as SB1 (with v sin i > 350 km s⁻¹), which might be a mass gainer in previous binary interactions with an undetected companion.

We do, however, note that there are no B-type SB1s rotating at velocities > 300 km s⁻¹. This is similar to the findings of Ramírez-Agudelo et al. (2015), who studied the primaries of O-type binaries in 30 Dor and found a lack of stars rotating with v sin i > 300 km s⁻¹, in comparison to single O-type stars. They interpret this as indication that the most rapidly rotating single O stars are actually binary interaction products (de Mink et al. 2013). While their sample is dominated by SB1s, primaries and secondaries in SB2s were identified based on the luminosity ratio, which might be misleading in the case of Algol systems (see e.g., Mahy et al. 2020). Investigating Galactic O-type stars, Holgado et al. (2022) found that the projected rotational velocity distributions of single stars and SB1s are different, and that the most rapid rotators (v sin i > 300 km s⁻¹) appear as single stars.

Another possible reason for the lack of very rapidly rotating SB1s could be an observational bias, with the more rapid rotators having broader spectral lines. In turn this leads to larger uncertainties in the RV estimates with only targets with large RV variations being identified. Further investigation of this possible observational bias is needed to confidently answer whether there are indeed no rapidly rotating SB1s in NGC 330.

The CMD colored by the derived v sin i values in Fig. 8 shows a similar trend as Fig. 7: the stars populating the Be star region rotate, on average, more rapidly than stars on the MS. Furthermore, the lower MS (m_F814W > 17.8 mag) seems to show a trend between projected rotational velocity and color: while stars at the blue part of the MS rotate on average slower, the stars at the redder part of the MS rotate more rapidly. In particular, there are also some rapidly rotating B and Be stars at a color of around −1.2 mag and m_F814W ∼ 18.5 mag. Those fainter stars usually have a lower S/N in their spectra, introducing a larger error in the derived rotational velocities.

Fig. 8. CMD of NGC 330 constructed from HST photometry (Milone et al. 2018), color-coded by v sin i. For a reference, we provide two Padova isochrones at 35 (dashed line) and 40 Myr (dotted line). B-type stars are indicated by circles while Be stars are indicated by diamonds. The color bar only includes v sin i values between 100 km s−1 and 400 km s−1 such that the few stars that rotate more rapidly fall in the last bin. Detected SB1s are indicated by red frames. The example star #654, peculiar stars and stars above the TO are marked by their ID. The four regions defined in Paper II are also indicated.

The variation of the v sin i along the MS band might indicate differences in the intrinsic rotational velocities of MS stars (as predicted by theory, see e.g., Wang et al. 2022a, and Sect. 4.1). Given that stars on the lower MS are generally unevolved, most of them are not spun-up mass gainers in post-interaction systems. However, some of them could be undetected binaries in which the spectral lines of the components are not fully deblended and thus mimic a larger observed rotational velocity.

Figure 8 furthermore shows that most single B-type stars above the cluster TO are slow rotators. Additionally, there are both slowly and moderately rotating Be stars above the cluster TO, which are both classified as single and SB1s. We discuss the projected rotational velocities of stars in different regions further in Sect. 3.4.

3.3. The Hertzsprung-Russell diagram

Despite the relatively large uncertainties in the obtained parameters, we use the effective temperatures and luminosities to construct an HRD (see Fig. 9). Overplotted are single-star evolutionary tracks and isochrones of 35 and 40 Myr from Georgy et al. (2013) assuming an initial rotational velocity of 50% critical. Given that the cluster is expected to contain a significant number of binary interaction products, they are mostly to guide the eye.

Fig. 9. HRD of B and Be stars in NGC 330. Presumably single B and Be stars are marked by a circle and diamond, respectively, while detected SB1s are indicated by symbols with red frames. The projected rotational velocity of is given by the color of the symbol. To guide the eye, we overplot evolutionary tracks from Georgy et al. (2013) assuming Z = 0.002, an initial rotation rate of 50% critical, and for different initial masses as indicated in the plot. We further overplot 35 Myr (dashed) and 40 Myr (dotted) isochrones with the same assumptions. Stars above the TO in the CMD and peculiar stars are marked by their ID. The limited availability of diagnostic spectral lines, photometric data points, and suitable models can lead to actual errors larger than the formal fitting errors indicated here.

Given the large uncertainties on the individual parameters, in particular log g, we refrain from computing spectroscopic mass (M_spec ∝ g/R², so our typical log g-uncertainties of ±0.25 dex translate into a factor of 2 in mass). A comparison of the star’s position with respect to the isochrones allows to roughly estimate their evolutionary mass.

In the HRD, we indicate the IDs of stars located above the cluster TO in the CMD (Fig. 8). This demonstrates that stars situated above the TO in the CMD are also above the TO in the HRD. In particular, the slowly rotating single B-type stars above the TO (namely #522, 289, 57, 638, 670, and 228), which are above the TO in the CMD and seem to follow a younger isochrone, also trace a similar, younger isochrone in the HRD. In contrast, the particularly faint and blue star #604 does not show peculiar stellar parameters. This is in line with the possible explanation we put forward in Paper I where we proposed that the visible star might have a stripped companion, which contributes to the colors but not to the spectral lines.

The evenly spaced values of T_eff reflect the steps in the TLUSTY grid used. The apparent lack of stars toward the isochrones at lower luminosities (i.e., log L/L_⊙) < 2.8) is probably mainly due to our selection criterion in magnitude, which preferably selects cooler stars at the same luminosity. At higher luminosities, several stars are to the left of the two isochrones (although consistent with the 35 Myr-isochrone within the errors). This might imply that rejuvenation took place. It could also be due to the choice of 50%-critical initial rotation in the isochrone and might imply that this is not a valid assumption for all the stars.

Several stars, also of lower mass, populate the post-MS region where they should not be situated according to standard single-star stellar evolution. While some of them could be binary-interaction products, the large errors in derived effective temperatures (i.e., ∼1000 K) and the uncertainties in the evolutionary tracks such as the adopted rotational velocity and overshooting most likely blur the picture. Undetected binary systems might lead to a further broadening of the observed MS given that the dilution by a companion star might light to a lower derived temperature while the derived luminosity is higher than for a single star. Therefore, and given the larger error bars, we focus on the CMD in the following, which benefits from the high accuracy of HST photometry and provides a clear overview of the cluster.

The HRD indicates that a majority of the rapidly rotating Be stars are toward the cooler end of the MS band. While the errors are larger for the Be stars, this may be due to a combination of two factors: Be star with a high projected rotational velocity are likely viewed equator-on. Due to the deformation of rapidly rotating stars, they have lower temperatures around their equator than around their poles (see e.g., von Zeipel 1924). Additionally, in the cases where the photometry fit is weighed more strongly, the lower derived temperature might reflect the impact of the circumstellar decretion disk for Be stars (Townsend et al. 2004). As noted previously noted in Paper I and by Hastings et al. (2021), a large number of Be stars are unevolved and far from the cluster TO, implying they are binary interaction products.

3.4. Rotational velocities for stars in different CMD regions

In Paper II we defined different regions in the CMD to investigate possible differences in the observed and bias-corrected binary fractions (indicated in Fig. 8). We found that the observed and bias-corrected spectroscopic binary fractions differ significantly in the different regions. Here, we investigate if the v sin i-distributions also vary from region to region.

The regions, namely the main sequence (MS), turn-off (TO), blue straggler star (BSS), and the Be star region, are solely defined based on the position of stars in the CMD in comparison to the 40-Myr Padova isochrone (see Fig. 8). The association of stars to a region does not take into account their spectral appearance. In particular, the Be star region is defined as a region offset redward of the MS, not based on the presence of emission lines in the spectra. Several spectroscopically classified Be stars fall out of it, while five faint B-type stars fall in the Be star region.

We note that this current division in regions considers the MS as one region, without further subdividing it in terms of a blue and red MS (as done for example by Wang et al. 2022a; Kamann et al. 2023), despite the fact that there might be a possible bi-modality in the v sin i of the faintest stars (see Sect. 3.2). This region is particularly interesting to investigate unevolved stars and pre-interaction binaries. In a future paper of this series, we will consider additional, fainter stars and perform this detailed investigation of the rotational velocities of lower MS.

Figure 10 shows the velocity distribution of stars in the four regions. It demonstrates that the rotational velocities of MS and TO stars are rather similar. The projected rotational velocities in the BSS region are generally below 300 km s⁻¹, which is lower than the stars in other regions. The stars in the Be star region show higher average velocities around v sin i ∼ 350 km s⁻¹.

Fig. 10. Projected rotational velocities of different stellar groups, defined based on the CMD (see Fig. 8). Top four panels: from top to bottom, the panels show the v sin i distributions of star in the MS region, the TO region, the BSS region and the Be star region. In each panel, the detected SB1s are marked by a darker color, as indicated in the legend. Bottom panel: normalized cumulative distribution of the rotation rates of different groups, as indicated in the legend.

To verify if the observed differences are statistically significant, we perform a Kuiper test between each of the groups (see Table 1). The difference in the rotational velocity distribution of Be stars compared to any other group is statistically significant (that is, the probability that the distributions are drawn from the same parent distribution is below 0.03% for any of them). Despite the appearances, the v sin i-distributions of MS, TO, and BSSs are not significantly different from one another, which may in part be due to the low number of stars (i.e., only 13) in the BSS region. When considering only v sin i above 100 km s⁻¹ (which we adopt as formal MUSE sensitivity limit), we find similar results: the projected rotational velocities in the MS and TO regions are similar, while the Be star region is statistically different. The BSS region only contains six stars after this cut and can therefore not be used for such a statistical comparison.

Similar to the MS region, the TO region probably contains stars of different natures (e.g., truly single stars, pre-interaction binaries, or around the TO also currently interacting binaries). The observed rotation rates in these two regions are similar and show a large spread. There furthermore seem to be no significant differences between presumably single stars and SB1s. We do, however, note that there are a large number of SB1s with rotational velocities between 100 and 200 km s⁻¹ in the TO region. These might coincide with a group of binary systems rotating at similar velocities detected by Ramírez-Agudelo et al. (2015), who interpreted them as being spun up through tides.

The BSS region is populated by two different types of systems: On the one hand, there are apparently single B-type stars with particularly low rotational velocities. Those could be merger products, which are predicted to be slow rotators (see e.g., Schneider et al. 2020). On the other hand, there is a group of SB1 systems and emission line stars with v sin i ∼ 250 km s⁻¹. Those could be currently interacting or post-interaction binary systems. In particular, the two SB1s with Balmer emission lines could be currently interacting binaries where the emission is formed in an accretion disk (possibly Algol-like systems in which the mass ratio was reversed). Some of the more slowly rotating SB1s might be binary systems in which tides may have led to a synchronization and therefore a spin-down.

Stars in the Be star region have a significantly different rotation rate distribution: on average, they are rotating more rapidly than stars in the other regions. This is in line with theoretical expectations: Be stars could be single stars approaching their critical rotation rate, or they could have gained their rapid rotation as mass gainers in previous binary interactions. The rapid rotation (and the additional IR emission from the circumstellar disk) leads to a redder color of the stars, populating the Be star region. We further note that there are four slowly rotating stars in the Be star region, which could appear as such because of inclination effects. Given the sample size of 47 targets in the Be star region, it seems, however, unlikely that four stars (i.e., almost 10%), are viewed under low inclination angles.

4. Comparison to previous observational studies

4.1. Rotational velocities in the LMC cluster NGC 1818

The cluster NGC 1818 is located in the LMC and has a similar age as NGC 330 (∼40 Myr, Milone et al. 2018). Similarly to NGC 330, NGC 1818 was reported to have an extended TO region. Additionally, NGC 1818 has a split MS (i.e., the main MS in NGC 1818 shows two components, in addition to the Be sequence), which is less pronounced in NGC 330. Given the lower metallicity in the SMC compared to the LMC, a comparison between NGC 330 and NGC 1818 might help assessing the impact of metallicity on the observed rotational velocity distributions.

Using VLT/FLAMES in GIRAFFE mode, 44 MS stars in NGC 1818 were studied spectroscopically and their rotational velocities were derived by comparing the observed Hα- and He I-line at λ6678 to synthetic spectra (Marino et al. 2018). All targets were selected to be on the MS (both the red and the blue MS) and span a brightness range m_F814W = 15.5 − 19.0.

The derived rotational velocities of the stars in NGC 1818 vary as a function of the brightness and color of the star: stars attributed to the blue MS rotate, on average, slower than stars in the red MS. Be stars are found to be the fastest rotators, while stars around the TO are generally slow rotators (see e.g., Fig. 7 in Marino et al. 2018). This agrees with theoretical predictions that explain the split MSs observed in young star clusters with different rotational velocities (e.g., Bastian & de Mink 2009; Wang et al. 2022a). A similar bi-model velocity distribution is also observed in the MS stars of the 100 Myr-old LMC cluster NGC 1850 (Kamann et al. 2023).

While these results are based on only a low number of stars in NGC 1818 (e.g., four stars in the TO region, and four Be stars), they qualitatively agree with the results derived here for NGC 330, despite the different metallicity of the two clusters. Here, no difference was made between the red and blue MS. Distinguishing those might have an impact on the rotational velocity distribution of the lower MS. Indeed, the rotational velocities of stars in the MS box in Fig. 8 appear to vary as a function of color: redder stars rotate, on average, faster. A more detailed comparison of the rotational velocities in the different MS components will, however, be part of a subsequent paper in this series.

4.2. Rotational velocities in 30 Dor and NGC 346

We compare our observed rotational velocities to the ones measured by Dufton et al. (2013) for early B-type stars in the 30 Dor region in the LMC in Fig. 11, excluding stars classified as O9.5 (see their Fig. 3 for more details). We adopt a bin size of 40 km s⁻¹, the same as chosen by Dufton et al. (2013), and include all B and Be stars in NGC 330 (including spectral types O9 to A0, see Paper I). We note that our spectral types are on average later than the ones studied by Dufton et al. (2013), who focused on stars with spectral types B0-B3, implying that the stars studied have different masses. We refrain from restricting our sample to those early spectral types because it would reduce the sample size in NGC 330 to only ∼40 stars.

Fig. 11. Normalized histogram of the projected rotational velocities of B and Be stars in NGC 330 (turquoise) compared to early-B and early-Be stars in 30 Dor (gray; Dufton et al. 2013), adopting the same bin size as in their Fig. 3.

Dufton et al. (2013) report a bi-modal v sin i distribution of B-type stars in 30 Dor, with ∼25% of the sample rotating slower than ∼80 km s⁻¹ and a second v sin i peak at around 200 km s⁻¹. Focusing on the O stars in 30 Dor, Ramírez-Agudelo et al. (2013) find a similar peak at ∼80 km s⁻¹, and a tail of high velocities up to 600 km s⁻¹. Our results indicate that there is no bi-modal v sin i distribution in NGC 330: the distribution is rather flat up to ∼300 km s⁻¹ after which it tails of. The apparent double-peaked feature around 200 km s⁻¹ is most likely due to the adopted binning (see also Fig. 7, where such feature is visible for the Be stars). Given the velocity-resolution of MUSE, however, a bi-modality cannot be ruled out. Indeed, it might be present in the lower MS (see Sect. 3.2).

Figure 11 furthermore indicates that the stars in NGC 330 rotate, on average, more rapidly than the stars in 30 Dor. We test the statistical significance of this finding by performing a Pearson χ² test, which indicates that this result is statistically significant (the likelihood that both v sin i-distributions are drawn from the same parent distribution is < 10⁻³). When only considering B-type stars and excluding Be stars in NGC 330 (which on average show more rapid rotation rates) we still find a similar result, though with a lower confidence.

We further compare the intrinsic rotational velocities of the B and Be stars in NGC 330 (see Sect. 2.7) with stars in 30 Dor and NGC 346 (see Fig. 12). NGC 346 is a ∼2 Myr-old, star-forming H II region in the SMC (Dufton et al. 2013, 2019, 2022). While NGC 346 is at a similar metallicity than NGC 330, the stellar populations in 30 Dor and NGC 346 are significantly younger than the one in NGC 330. We include all stars classified in the three samples as single stars and SB1s and note that the different observing approaches have different binary detection probabilities. We again note that the different studies target a slightly different range of spectral types, or stellar masses: stars in NGC 330 have spectral types from B0 to B9 (see Paper I), while the B-star sample in NGC 346 focuses on stars between B0 and B3, similar to the 30-Dor sample.

Fig. 12. Projected and intrinsic rotational velocity distributions of B and Be stars in NGC 330. The normalized histogram in each panel shows the observed projected rotational velocities in NGC 330, while the solid line indicates the intrinsic distribution of rotational velocities, that is corrected for random inclination angles. The dotted and dashed lines show the intrinsic rotational velocity distributions for stars in 30 Dor and NGC 346, respectively (Dufton et al. 2013, 2019, 2022). The number of stars in each sample is indicated in brackets. The top panel shows B and Be stars combined, the middle panel focuses on B-type stars, while the bottom panel shows the Be stars only.

The intrinsic rotational velocity distribution of the entire NGC 330 sample shows no evidence for a significant bi-modality, but rather a peak at 300 km s⁻¹ (see Fig. 12). The comparison of B and Be stars demonstrates again that the Be stars rotate, on average, more rapidly than the B-type stars: their intrinsic distribution peaks around 400 km s⁻¹, with only very few stars rotating slower than 150 km s⁻¹. The local maxima in the B- and Be-star distributions are most likely due to the binning of the rotational velocities and do not represent true features.

In contrast, the overall sample of stars in 30 Dor (and to a lesser extent in NGC 346) seem to have bimodal rotational velocity distributions. In general, there are fewer stars in NGC 330 with rotational velocities below 100 km s⁻¹, and more rapidly rotating stars, in particular compared to NGC 346. The median rotational velocity of the stars in NGC 346 is significantly lower (i.e., ∼150 km s⁻¹) than the one in 30 Dor and NGC 330 (∼230 km s⁻¹ in both clusters). Figure 12 further indicates that the Be stars, which follow unimodal rotational velocitity distributions in all three populations, are more rapidly rotating in NGC 330 (their median rotational velocity is 350 km s⁻¹) than the ones in the other two populations, and in particular than the ones in NGC 346 (with a median velocity of ∼300 km s⁻¹).

4.3. Discussion

The comparison of the derived rotational velocities in the different observed samples is complicated by several aspects. Firstly, as mentioned previously, the mass ranges targeted by the different studies differ from each other. While the samples in NGC 346 and 30 Dor are comprised mainly of early-B type stars (O9.5 to B3, corresponding to approximately 16 to 8 M_⊙, Dufton et al. 2013, 2019), the stars in targeted in NGC 330 and NGC 1818 are primarily of later spectral types (O9 to B9, corresponding to 16 to 3 M_⊙, with a peak at spectral type B5, see Fig. 14 in Paper I). The mass of a star impacts the strength of the stellar wind, which can spin-down stars if the mass-loss rates are high enough. However, as B-type stars show only weak winds, this should not lead to significant differences in the v sin i-distributions.

Another caveat in the comparison is that we consider populations of stars are at different metallicities. While NGC 330 and NGC 346 are both in the SMC (with a metallicity Z ∼ 0.2 Z_⊙), 30 Dor and NGC 1818 are in the LMC (with Z ∼ 0.5 Z_⊙). As discussed in Sect. 1, stars at lower metallicity should on average have higher projected rotational velocities because they are more compact and have weaker winds. The second effect, as mentioned above, is barely relevant for B-type stars. This is supported by the fact that we find generally higher rotational velocities for stars in NGC 330 compared to 30 Dor and NGC 346, while there is a general agreement with the stars in NGC 1818.

A third difference between the studies are the environments they target. As discussed in Sect. 1, Galactic studies found that cluster stars rotate on average faster than field stars (e.g., Strom et al. 2005; Daflon et al. 2007; Garmany et al. 2015). These differences were attributed either to differences in the star formation process (cf. Wolff et al. 2007) or to evolutionary effects (Huang & Gies 2008). While the studies on NGC 330 and NGC 1818 focus on the central core of the clusters, the study of 30 Dor targets a more extended star-forming region. In contrast, stars in NGC 346 are mostly field stars (Dufton et al. 2019). The latter could be subject to sequential or continuous star formation. The lower v sin i in NGC 346 could at least in part be explained by the stars being predominantly field stars.

The last and main difference between the samples is their age: while NGC 346 and 30 Dor contain very young stars (i.e., < 5 − 15 Myr, Dufton et al. 2013), the population of stars in NGC 330 and NGC 1818 are roughly 35−40 Myr old (which is also the reason for the later spectral types in NGC 330 and NGC 1818). Most stars and binary systems in NGC 346 and 30 Dor are unevolved, and the parameters detected there represent the initial parameters. In contrast, the more massive stars in NGC 1818 and NGC 330 have evolved off the MS. If in close-enough binary systems, those systems interacted, either by mass transfer or via a stellar merger. NGC 330 and NGC 1818 is therefore thought to be populated by such post-interaction products, either the spun-up mass gainers of mass transfer, or the (potentially slowly rotating) merger products (e.g., de Mink et al. 2013; Schneider et al. 2018). In the case of mass transfer, their previously more massive companions could now be compact objects or might have exploded as supernova, and the systems are now detected as single stars. Due to rejuvenation, such systems are expected to be situated around or above the cluster TO.

These expectations are in line with our observational findings in NGC 330: the rapidly rotating, mainly single Be stars could be such spun-up companion stars. The slowly rotating, apparently single stars above the cluster TO could be merger products. Finally, the moderately rotating stars above the cluster TO, among which some show emission lines and some are detected as binaries, could be currently interacting binary systems. While our observations and the theoretical predictions of stellar spins match qualitatively, other factors such as the metallicity of the stars might play an important role. A more quantitative comparison, which we will perform in a subsequent paper of this series, will allow the different aspects to be disentangled further.

5. Summary and conclusions

For this work, we have measured the stellar parameters of the presumably single stars and SB1s in NGC 330. For this purpose we used the available HST photometry as well as spectra that were combined from the six observed MUSE epochs. We determined stellar parameters with a two-component grid-search approach using the BSTAR2006 TLUSTY grid of model atmospheres.

The high quality of the HST observations allowed us to constrain radii and luminosities. Given the availability of only three filters and the fixed extinction to a common value for all stars, systematic uncertainties may remain. The low-resolution spectra lead to a larger scatter (accompanied by larger relative uncertainties) in T_eff and log g. Rotational velocities can be constrained even better than the formal MUSE resolution limit of around 100 km s⁻¹ when including several lines in the fit. In general, the derived stellar parameters agree with the expectations from the star’s position in the CMD. Furthermore, we found no significant differences in the determination of the parameters of previously identified SB1 systems, apart from a larger scatter and larger uncertainties that could be due to a possible contamination of the observables by the unseen companion.

Our observations show that a majority of B-type stars have v sin i ≈ 100 − 250 km s⁻¹ with a tail to higher velocities. The projected rotational velocities of Be stars are on average higher and show a peak at around 200−400 km s⁻¹. A Kuiper test demonstrates that the two distributions are statistically different. We further associated different regions in the CMD to different stellar populations: stars on the lower main sequence are associated with the MS region, stars around the turnoff to the TO region, blue straggler stars above the cluster turnoff to the TO region, and stars redward of the MS to the Be region. We found that the projected rotational velocities of the stars in the Be region are significantly larger than the ones of stars on the MS, and in the TO and BSS regions. While BSSs seem to, on average, rotate more slowly than MS and TO stars, small sample sizes prevent us from confirming the significance of this result. We note that a subset of the BSS stars are apparently single, slowly rotating stars, which could be the products of stellar mergers.

Comparing the stars in NGC 330 in other stellar populations, we find an overall agreement with stars in NGC 1818, a LMC cluster of similar age. We further find that the stars in NGC 330 generally rotate more rapidly than O- and early-B-type stars in the young 30 Dor region in the LMC (Ramírez-Agudelo et al. 2013; Dufton et al. 2013, 2019), and more rapidly than B-type stars toward the young SMC cluster NGC 346. Given that the latter two studies mostly targeted B-type MS stars, which in general have weak winds, spin-down due to stellar winds should not be able to explain the observed differences.

We propose that the higher rotation rates in NGC 330 could be rooted in the age difference between the clusters. In the young clusters (30 Dor and NGC 346), the late-O and early-B stars will be relatively unevolved, while most binary systems will not have interacted. By contrast both of these phenomena may be important for the older cluster NGC 330, where more massive stars had time to evolve and interact. Mass transfer in close binary systems might lead to a higher fraction of Be stars, which could be the mass gainers in previous binary interactions. The brightest Be stars close to the TO could, however, also be single stars with sufficiently high initial rotational velocities to spin up toward the end of their MS evolution. The slowly rotating, apparently single stars above the cluster TO could be merger products. This is in qualitative agreement with binary population synthesis predictions (see e.g., Wang et al. 2020).

In this work, we have provided an extensive and homogeneous dataset which can be directly compared to modern population synthesis models including the most recent treatment of binary interaction physics (e.g., Eldridge & Stanway 2009; Wang et al. 2020, 2022b; Klencki et al. 2022). Such a comparison will help to constrain currently implemented single- and binary-star evolution physics and allow for a better understanding of massive star evolution. A possible future step to further investigate the properties of (at least the brightest) cluster members is to study their chemical abundances. In particular the BSSs, which are excellent merger candidates, might show CNO-processed material on their surface such as nitrogen enhancement. Despite the low resolving power of MUSE, corresponding lines in the optical regime can be detected if they are strong. Detecting elevated N abundances in BSS stars will give a further hint as to the binary interaction history of these stars, and provide important insights for binary merger models. Future work could also include a more detailed investigation of the rotational velocities of different MS components, mainly the red and blue MS, which are predicted to be composed of rapidly and slowly rotating stars, respectively. Additional constraints could also be provided by the investigation of additional young and intermediate-age clusters in the SMC, LMC and the Milky Way.

Acknowledgments

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 and 945806: TEL-STARS).

References

Appendix A: Stellar parameters derived in this work

Table A.1 gives an overview of the derived stellar parameters for each star with respective errors, namely the effective temperature, surface gravity, projected rotational velocity, radius and luminosity. We also repeat here the spectral types and binary classification from Paper I and Paper II.

Appendix B: A second test case: The Be star #688

Given that our test case #654 was a slowly rotating B-type star, we here show a second star, a rapidly rotating Be star, namely #688 to demonstrate how the method works for rapid rotators. Figures B.1, B.2, and B.3 are analogous to Figures 2, 3, and 4 for #654.

Fig. B.1. Spectroscopic fit for star #688, similar to Fig. 2 for #654. Each panel, one for each of the six diagnostic spectral regions, shows the combined spectrum (black) and the best-fitting TLUSTY model (red). Grayed-out regions are only shown for comparison but were not included in the fit.

Fig. B.2. Photometric fit for star #688, similar to Fig. 3 for #654. The top panel shows a comparison between the observed HST fluxes (black crosses) and the flux (colored circles) computed from the best-fit model (gray line) by convolution with the HST filters (colored lines). The bottom panel gives residuals.

Fig. B.3. Combined fit for star #688, similar to Fig. 4 for #654. The diagonal gives the χ2-distribution as a function of each parameter and the panels below show 2-dimensional χ2-maps. The red line marks the 95% confidence level.

Given the strong contamination of the spectrum by emission, the affected lines (in this case mainly Hβ, are excluded from the fit). However, the fit of the He I line at 7065.19 Å might imply that there is some infilling as well. The Hα and Paschen lines show strong contamination, but those lines are excluded from the fit for all stars (see Sect. 2.2). Whether the photometry is impacted by the Be disk is hard to assess from three photometric points only. Given that the HST observations were not taken simultaneously, we cannot extrapolate from the MUSE spectra either.

The combined fit suggests that the Be star #688 has the following parameters with uncertainties: T_eff = 18000${}_{-1500}^{+2100}$ K, log g = 3.0${}_{-\dots }^{+0.3}$ dex, v sin i = 200${}_{-60}^{+50}$ km s⁻¹, R = 7.4${}_{-1.1}^{+0.8}$ R_⊙, and log(L/L_⊙) = $3.7_{-0.3}^{+0.3}$ L_⊙. As the best-fit value for the surface gravity hits the edge of the grid, no lower error can be given. On average, the error bars are higher than the ones determined for most B-type stars. This is due to the fact that fewer spectral lines are used in the fit. The correlation plot in Fig. B.3 further demonstrates that the log g is not well constrained when no Balmer lines (i.e., in particular Hβ in our case) can be included in the fit.

Appendix C: Stellar parameters from Carini et al. (2020)

Based on one epoch of MUSE-WFM-NOAO commissioning data obtained with the nominal wavelength range in August 2014, Carini et al. (2020) determine Johnson V and Cousins I magnitudes from the MUSE data and extract spectra via PSF for 10 bright B-type stars in the core of NGC 330. Using photometric estimates as input, they compare the MUSE spectroscopy in the wavelength range from λλ 4800-5000Å with synthetic models to estimate stellar parameters and measured He abundances. In the following, we compare our derived parameters with the spectroscopic parameters provided by Carini et al. (2020), noting that their average S/N is slightly lower than the one used in this study. Furthermore, given that the observations were taken without AO support and with an average seeing of 1.3″, crowding is expected to be a more severe problem here.

Given a similar pointing, all of the ten stars are included in our study. Two of the ten stars, namely #313 (A10) and #519 (A48) were classified as SB1s in Paper II. Because we do not detect any signature of a companion in their spectrum we fit the spectrum like the one of a single star. The comparison of the derived parameters (see Tab. C.1) shows that they are in reasonable agreement.

The largest discrepancy in the estimated temperatures and surface gravities is found for the SB1 #519 which might imply that there is some contamination of the secondary that manifests itself differently in the two different fitting processes, for example due to the different wavelength region considered. Star #420 also shows a large discrepancy in the derived parameters. This might be due to the fact that #420 is a Be star with strong emission line contamination in the diagnostic spectral lines, and might imply that at least for those stars, the errors are underestimated. Given that a strong He II absorption line is present in the spectrum of #420 the values derived in this work seem more plausible. We further note that the temperature estimated for the third hottest star, #522, is also significantly lower in Carini et al. (2020). The temperature estimates of all three stars could also be impacted by the fact that no He II was included in the analysis of Carini et al. (2020).

While Carini et al. (2020) only considered rotational velocities above 100 km s⁻¹ in steps of 50 km s⁻¹, the derived rotational velocities follow a similar trend. They furthermore do not provide any uncertainties on v sin i, which hampers the comparison.

All Tables

All Figures

	Fig. 1. Example of the co-addition of spectra for the apparently single star #654. Individual epochs are shown in color (as indicated in the legend) while the combined spectrum is shown in black. The legend also indicates the measured S/N at 6100 Å. Main spectral lines are marked.
In the text

	Fig. 2. Spectroscopic fit for star #654. Each panel, one for each of the six diagnostic spectral regions, shows the combined spectrum (black) and the best-fitting TLUSTY model (red). Grayed-out regions are only shown for comparison but were not included in the fit.
In the text

	Fig. 3. Photometric fit for star #654. The top panel shows a comparison between the observed HST fluxes (black crosses) and the flux (colored circles) computed from the best-fit model (gray line) by convolution with the HST filters (colored lines). The bottom panel gives residuals.
In the text

	Fig. 4. Combined fit for star #654. The diagonal gives the χ2-distribution as a function of each parameter and the panels below show 2-dimensional χ2-maps. The red line marks the 95% confidence level.
In the text

	Fig. 5. Relation between the FWHM of the He I line at 4922 Å with the v sin i estimated from the TLUSTY fit (purple line, see Eq. (4)). The purple points indicate the measured values for B-type stars.
In the text

	Fig. 6. CMDs of NGC 330 color-coded by the best-fitting physical parameters. Circles indicate presumably single stars while triangles indicate stars classified as SB1s. The isochrones are single-star nonrotating Padova isochrones of 35 (dashed line) and 40 Myr (dotted line).
In the text

	Fig. 7. Stacked histogram of the projected rotational velocities of presumably single stars and SB1s. The top panel shows the B-type stars (single stars in yellow, SB1s in orange) and the bottom panel shows Be stars (single stars in light green, SB1s in dark green). The gray shaded region again corresponds to the formal MUSE sensitivity limit while the indicated errors are counting errors.
In the text

Fig. 8. CMD of NGC 330 constructed from HST photometry (Milone et al. 2018), color-coded by v sin i. For a reference, we provide two Padova isochrones at 35 (dashed line) and 40 Myr (dotted line). B-type stars are indicated by circles while Be stars are indicated by diamonds. The color bar only includes v sin i values between 100 km s−1 and 400 km s−1 such that the few stars that rotate more rapidly fall in the last bin. Detected SB1s are indicated by red frames. The example star #654, peculiar stars and stars above the TO are marked by their ID. The four regions defined in Paper II are also indicated.

In the text

Fig. 9. HRD of B and Be stars in NGC 330. Presumably single B and Be stars are marked by a circle and diamond, respectively, while detected SB1s are indicated by symbols with red frames. The projected rotational velocity of is given by the color of the symbol. To guide the eye, we overplot evolutionary tracks from Georgy et al. (2013) assuming Z = 0.002, an initial rotation rate of 50% critical, and for different initial masses as indicated in the plot. We further overplot 35 Myr (dashed) and 40 Myr (dotted) isochrones with the same assumptions. Stars above the TO in the CMD and peculiar stars are marked by their ID. The limited availability of diagnostic spectral lines, photometric data points, and suitable models can lead to actual errors larger than the formal fitting errors indicated here.

In the text

Fig. 10. Projected rotational velocities of different stellar groups, defined based on the CMD (see Fig. 8). Top four panels: from top to bottom, the panels show the v sin i distributions of star in the MS region, the TO region, the BSS region and the Be star region. In each panel, the detected SB1s are marked by a darker color, as indicated in the legend. Bottom panel: normalized cumulative distribution of the rotation rates of different groups, as indicated in the legend.

In the text

	Fig. 11. Normalized histogram of the projected rotational velocities of B and Be stars in NGC 330 (turquoise) compared to early-B and early-Be stars in 30 Dor (gray; Dufton et al. 2013), adopting the same bin size as in their Fig. 3.
In the text

Fig. 12. Projected and intrinsic rotational velocity distributions of B and Be stars in NGC 330. The normalized histogram in each panel shows the observed projected rotational velocities in NGC 330, while the solid line indicates the intrinsic distribution of rotational velocities, that is corrected for random inclination angles. The dotted and dashed lines show the intrinsic rotational velocity distributions for stars in 30 Dor and NGC 346, respectively (Dufton et al. 2013, 2019, 2022). The number of stars in each sample is indicated in brackets. The top panel shows B and Be stars combined, the middle panel focuses on B-type stars, while the bottom panel shows the Be stars only.

In the text

	Fig. B.1. Spectroscopic fit for star #688, similar to Fig. 2 for #654. Each panel, one for each of the six diagnostic spectral regions, shows the combined spectrum (black) and the best-fitting TLUSTY model (red). Grayed-out regions are only shown for comparison but were not included in the fit.
In the text

	Fig. B.2. Photometric fit for star #688, similar to Fig. 3 for #654. The top panel shows a comparison between the observed HST fluxes (black crosses) and the flux (colored circles) computed from the best-fit model (gray line) by convolution with the HST filters (colored lines). The bottom panel gives residuals.
In the text

	Fig. B.3. Combined fit for star #688, similar to Fig. 4 for #654. The diagonal gives the χ2-distribution as a function of each parameter and the panels below show 2-dimensional χ2-maps. The red line marks the 95% confidence level.
In the text