The young massive SMC cluster NGC 330 seen by MUSE III

Context. The origin of the initial rotation rates of stars, and how a star’s surface rotational velocity changes during the evolution, either by internal angular momentum transport or due to interactions with a binary companion, remain open questions in stellar astrophysics. Aims. Here, we aim to derive the physical parameters and study the distribution of (projected) rotational velocities of B-type stars in the ∼ 35Myr-old, massive cluster NGC330 in the Small Magellanic Cloud. NGC330 is in an age range where the number of post-interaction binaries is predicted to be high near the cluster turno ﬀ (TO). Methods. We developed a simultaneous photometric and spectroscopic grid-ﬁtting method adjusting atmosphere models on multi-band Hubble Space Telescope (HST) photometry and Multi Unit Spectroscopic Explorer (MUSE) spectroscopy. This allowed us to homogeneously constrain the physical parameters of over 250 B and Be stars (i.e., B-type stars with emission lines), brighter than m F 814 W = 18 . 8mag. Results. The rotational velocities of Be stars in NGC330 are signiﬁcantly higher than the ones of B-type stars. The rotational velocities vary as a function of the star’s position in the color-magnitude diagram, qualitatively following predictions of binary population synthesis. A comparison to younger clusters shows that stars in NGC330 rotate more rapidly on average. Conclusions


Introduction
Apart from the initial mass, the initial rotation rate is considered to be one of the key properties with 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 ⋆ Full Table A.1 is only available in electronic form at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/cgi-bin/qcat?J/A+A/.⋆⋆ Based on observations collected at the ESO Paranal observatory under ESO program 60.A-9183(A) and 0102.D-0559(A).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 binarystar 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) showed 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.THhose might observationally appear as classical Be stars (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 towards 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) showed that stars in the SMC have, on average, higher rotational velocities.This was attributed to two effects: firstly, lower-metallicity stars are more compact and, at a given mass, thus have generally 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 Doradus 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 ( sin i) distribution.A majority of stars rotate at around 80 km s −1 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 shortperiod systems are faster rotators than single O-type stars likely as a result of tidal synchronisation (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 prior interaction is needed to produce such rapid rotation.Dufton et al. (2013) showed that the 30 Doradus B-star sin i-distribution is bi-modal, with one peak below 100 km s −1 and a high-velocity component above 250 km s −1 .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 on 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, andSect. 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 cluster (see e.g., Wang et al. 2020Wang et al. , 2022a,b;,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 described our multi-epoch observations on NGC 330 obtained with the Multi Unit Spectroscopic Explorer (MUSE) and characterized the massive star content.In Bodensteiner et al. (2021, hereafter Paper II) we measured radial velocities (RVs) and investigated the current binary properties of the massive star population.Here, we use 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.
Article number, page 2 of 18 Bodensteiner et al.: The young massive SMC cluster NGC 330 seen by MUSE.III.

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 sin i (where 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).
The grid provides synthetic spectra for T eff between 15000 to 30000 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 2km s −1 .The synthetic spectra cover the UV to optical wavelength range (900-10000 Å).
Our newly developed fitting routine 1 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), and 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.

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 (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 analysed 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 ratio 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/N ratios between approximately 125 and 290 at around 6100 Å) and the combined spectrum with a S/N of over 450 is shown in Fig. 1.

Spectroscopic fit
The observed MUSE spectra allow us to constrain T eff , log g and 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 1 available on github: https://github.com/jbodenst/here consider rotational velocities ranging from 0 to 500 km s −1 with steps of 20 km s −1 , which is chosen to be lower than the typical accuracy of the sin i-estimates as discussed below.We fit T eff , log g, and 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 χ 2 spec for each model as 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, and S/N is the measured S/N ratio.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 pointspread function (PSF), especially for fainter stars.
We define six diagnostic regions in the MUSE spectral range (also shown in Fig. 2): 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), the region around the He ii line at 5411 Å, the region around Hα, between 6500 and 6640 Å, the He i line at 6678.15 Å, the He i line at 7065.19 Å, and the Paschen region between 8500 and 8930 Å.
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 χ 2 spec -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 χ 2 spec 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 −1 , respectively.Previous studies based on MUSE data have, however, shown that a determination of sin i with a 1σ accuracy of 30 km s −1 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 micro = 2 km s −1 .As typical values of micro are of the order of a few km s −1 for B-type MS stars, this is much lower than the typical sin i and thus a valid assumption.Furthermore we assume that the macroturbulent velocity macro is negligible.The macroturbulent velocity for B-type stars is below 100 km s −1 in most cases (see e.g., Simón-Díaz et al. 2017).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.

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 broadband filters covering the SED from the near-UV to the red edge of the optical range, namely F225W, F336W, andF814W (centered at 2250 Å, 3360 Å, and8140 Å, respectively, Milone et al. 2018).The catalogue 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 bestfitting 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 filter2 .
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 smallscale 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 χ 2 phot for each model in the 3D-grid by comparing the three observed fluxes with the computed model fluxes according to: 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).

Combination of the two fits
We combine the two separate fits by adding up the computed χ 2 for each model with the same T eff and log g, but different sin i and radius: 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, sin i and stellar radius is found by locating the minimum in the obtained χ 2 global distribution for each parameter.We rescale the obtained χ 2 distribution such that the model with the minimum χ 2 has a χ 2 global 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 χ 2 .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 give 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 sin i are of the order of 1000 K, 0.25 dex, and 70 km s −1 , respectively.While the sin i uncertainty is in principle higher for higher rotational velocities, the dominating factor is the S/N ratio in the spectra.The errors derived for Be stars are, on average, higher than for Btype 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 =10000 K (L.Mahy, priv. comm., covering the same parameter range in log g and 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 ef-fective 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, 3, and 4) are available at the CDS.

A test case: star #654
We provide a proof of concept for one star, #654, and use it to test the assumptions on 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 χ 2 distributions derived of the combined fit can be seen in Fig. 4.
We find that # 654 has an effective temperature of 21000 +300 −2300 K a surface gravity of 3.8 +0.1 −0.3 dex, and a projected rotational velocity below 40 km s −1 .Given the estimated radius of 6.8 +1.3 −0.6 R ⊙ , the luminosity of the star is log(L/L ⊙ ) = 3.9 +0.2 −0.3 .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 = 18450 K from the spectral type, and based on this inferred log g = 3.45, and sin i = 20 km s −1 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.
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), andHunter et al. (2008) report T eff = 21200 K (again inferred from the spectral type), log g = 3.25 dex, sin i = 39 km s −1 , and log(L/L ⊙ ) = 4.33.In our approach, we find comparable parameters, namely T eff = 24000 +700 −200 K, log g = 3.5 +0.1 −0.1 dex, sin i < 50 km s −1 , 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 fitted the published FLAMES spectra of #654 and #57 with a Article number, page 5 of 18 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 Sections 2.2 and 2.3, we assume a constant E(B-V) = 0.08, A V = 3.1 and neglect 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 = 28000 K).Such a high temperature is excluded by the spectroscopic fit and the combined fit, which for this star is weighted towards the spectroscopic fit, does not change significantly.
Finally, we use iacob-broad (Simón-Díaz & Herrero 2014), which determines the rotational velocity of a star based on a combined Fourier transform and goodness-of-fit method taking into account a non-zero macro .We find sin i= 36 +46 −11 km s −1 and macro = 81 +23 −81 km s −1 .This again shows that sin i can be constrained better than its formal accuracy.Further, the derived sin i agrees with the one estimated in our approach, and the derived macro is consistent with our assumption of macro =0 km s −1 .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.

Estimation of 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 ve-Article number, page 6 of 18 locity.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 fitted 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 sin i measurements from the fitting described in Sect.2.2, we derive a FWHM-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 −1 .
The measured FWHM [Å] and sin i [km s −1 ], shown in Fig. 5, follow a roughly linear relation.Fitting it with a firstorder polynomial, weighted by the measured errors in FWHM and sin i , we find: with an RMS residual of 43km s −1 .As a consistency check, we repeat the fit for stars with fractional sin i errors from the tlusty fit <0.

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 Doradus (Dufton et al. 2013(Dufton et al. , 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 sin i-estimates for both the single and SB1 subsamples were combined using the following sources: Dufton et al. (2013Dufton et al. ( , 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 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 data sets, the implied stochastic uncertainties are significant and could lead to spurious small-scale structure.Nevertheless, we compare the intrinsic rotational velocity distribution to those derived for B and Be stars in 30 Doradus and NGC 346 (Dufton et al. 2013(Dufton et al. , 2019(Dufton et al. , 2022) ) in Sect.4.2.

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-sin icorrelation described above.In total, this provides us with stellar parameters of 251 stars and sin i-measurements of additional thirteen stars (264 stars in total), including presumably single stars and SB1s.

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.
Article number, page 7 of 18 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).

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 sin i of most B-type stars is around 100-250 km s −1 with a tail to high velocities up to 500 km s −1 .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 −1 .
The median projected rotational velocity of the B-type stars is around 170 km s −1 while it is at ∼250 km s −1 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 −3significance level.The observed projected rotational velocity distributions of B and Be stars are therefore statistically differ-ent 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 sin i >350 km s −1 ), 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 −1 .This is similar to the findings of Ramírez-Agudelo et al. (2015), who studied the primaries of Otype binaries in 30 Doradus and found a lack of stars rotating with sin i > 300 km s −1 , 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 ( sin i > 300 km s −1 ) 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 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.
The variation of the 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 will discuss the projected rotational velocities of stars in different regions further in Sect.3.4.

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 Myrs 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.
Given the large uncertainties on the individual parameters, in particular log g, we refrain from computing spectroscopic mass (M spec ∝ g/R 2 , 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 towards 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.
Article number, page 10 of 18 The HRD indicates that a majority of the rapidly rotating Be stars are towards 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.

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 sin idistributions 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 bimodality in the 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 velocties 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 −1 , which is lower than the stars in other regions.The stars in the Be star region show higher average velocities around sin i∼ 350 km s −1 .
To verify if the observed differences are statistically significant, we perform a Kuiper test between each of the groups (see Tab. 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 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 sin i above 100 km s −1 (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).
Table 1.Kuiper-probability that observed sin i distributions of different populations are drawn from the same parent distribution.
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 −1 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 sin i∼ 250 km s −1 .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.

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 Myrs, 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 iline at λ6678 to synthetic spectra (Marino et al. 2018) 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 were found to 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 was 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   ison of the rotational velocities in the different MS components will, however, be part of a subsequent paper in this series.

Rotational velocities in 30 Doradus 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 Doradus region in the LMC in Fig. 11, excluding stars classified as O9.5 (see their figure 3 for more details).We adopt a bin size of 40 km s −1 , 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.Dufton et al. (2013) reported a bi-modal sin i distribution of B-type stars in 30 Doradus, with ∼ 25% of the sample rotating slower than ∼ 80 km s −1 and a second sin i peak at around 200 km s −1 .Focusing on the O stars in 30 Doradus, Ramírez-Agudelo et al. (2013) found a similar peak at ∼ 80 km s −1 , and a tail of high velocities up to 600 km s −1 .Our results indicate that there is no bi-modal sin i distribution in NGC 330: the distribution is rather flat up to ∼300 km s −1 after which it tails of.The apparent double-peaked feature around 200 km s −1 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 Doradus.We test the statistical significance of this finding by performing a Pearson χ 2 test, which indicates that this result is statistically significant (the likelihood that both sin i-distributions are drawn from the same parent distribution is < 10 −3 ).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 Doradus 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(Dufton et al. , 2019(Dufton et al. , 2022)).While NGC 346 is at a similar metallicity than NGC 330, the stellar populations in 30 Doradus and NGC 346 are signif-   (Dufton et al. 2013(Dufton et al. , 2019(Dufton et al. , 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.
icantly 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.
The intrinsic rotational velocity distribution of the entire NGC 330 sample shows no evidence for a significant bimodality, but rather a peak at 300 km s −1 (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 −1 , with only very few stars rotating slower than 150 km s −1 .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 −1 , and more rapidly ro-tating 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 −1 ) than the one in 30 Doradus and NGC 330 (∼230 km s −1 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 −1 ) 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 −1 ).

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. 2013Dufton et al. , 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 figure 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 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 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 closeenough 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 previ-Article number, page 13 of 18 ously 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 to further disentangle the different aspects.

Summary and conclusions
In 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 by a two-component grid-search approach using the BSTAR2006 tlusty grid of model atmospheres.
The high quality of the HST observations allows 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 −1 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 sin i≈ 100−250 km s −1 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 −1 .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 to 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 Doradus region in the LMC (Ramírez-Agudelo et al. 2013;Dufton et al. 2013Dufton et al. , 2019)), and more rapidly than Btype stars towards 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, whilst 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 towards 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 provided an extensive and homogeneous data set 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. 2020Wang et al. , 2022b;;Klencki et al. 2022).Such a comparison will help to constrain currently implemented single-and binary-star evolution physics and allow 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 BSS, 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 for the binary interaction history of these stars, and provide important insights for binary merger models.Future work also includes 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 can also be provided by the investigation of additional young and intermediate-age clusters in the SMC, LMC and the Milky Way.

Fig. 2 .
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.

Fig. 3 .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.

Fig. 5 .
Fig. 5. Relation between the FWHM of the He i line at 4922Å with the sin i estimated from the tlusty fit (purple line, see Eq. 4).The purple points indicate the measured values for B-type stars .
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 non-linear relation deviates from our linear relation for lower rotational velocities (below ∼150 km s −1 ).Otherwise, the obtained rotational velocities using the relation by Ramírez-Agudelo et al. (2015) agree with ours within the RMS.

Fig. 6 .Fig. 7 .
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 non-rotating Padova isochrones of 35 (dashed line) and 40 Myr (dotted line).Article number, page 8 of 18

Fig. 8 .
Fig. 8. CMD of NGC 330 constructed from HST photometry (Milone et al. 2018), color-coded by sin i.For a reference, we provide two Padova isochrones at 35 (dashed line) and 40 Myr (dotted line).Btype stars are indicated by circles while Be stars are indicated by diamonds.The color bar only includes 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.

Fig. 9 .
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.

Fig. 10 .
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 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.

Fig. 11 .
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 Doradus (gray) (Dufton et al. 2013), adopting the same bin size as in their fig.3.

Fig. 12 .
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 Doradus and NGC 346, respectively(Dufton et al. 2013(Dufton et al. , 2019(Dufton et al. , 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.
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 ratio at 6100 Å. Main spectral lines are marked.
appear to vary as a function of color: redder stars rotate, on average, faster.A more detailed compar- Article number, page 12 of 18Bodensteiner et al.:The young massive SMC cluster NGC 330 seen by MUSE.III.