The effects of stellar rotation along the main sequence of the 100 Myr old massive cluster NGC 1850

Young star clusters enable us to study the effects of stellar rotation on an ensemble of stars of the same age and across a wide range in stellar mass and are therefore ideal targets for understanding the consequences of rotation on stellar evolution. We combine MUSE spectroscopy with HST photometry to measure the projected rotational velocities (Vsini) of 2,184 stars along the split main sequence and on the main sequence turn-off (MSTO) of the 100 Myr-old massive (10^5 M_sun) star cluster NGC 1850 in the Large Magellanic Cloud. At fixed magnitude, we observe a clear correlation between Vsini and colour, in the sense that fast rotators appear redder. The average Vsini values for stars on the blue and red branches of the split main sequence are ~100 km/s and ~200 km/s, respectively. The values correspond to about 25-30% and 50-60% of the critical rotation velocity and imply that rotation rates comparable to those observed in field stars of similar masses can explain the split main sequence. Our spectroscopic sample contains a rich population of ~200 fast rotating Be stars. The presence of shell features suggests that 23% of them are observed through their decretion disks, corresponding to a disk opening angle of 15 degrees. These shell stars can significantly alter the shape of the MSTO, hence care should be taken when interpreting this photometric feature. Overall, our findings impact our understanding of the evolution of young massive clusters and provide new observational constraints for testing stellar evolutionary models.


INTRODUCTION
Contrary to the Milky Way, the Large and Small Magellanic Cloud are known to host populations of young (< 1 Gyr) and intermediateage (few Gyrs) massive stellar clusters. Due to their high masses (∼ 10 5 M ) and low extinction, along with their relative proximity allowing for high precision photometry, these clusters have opened ★ Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 0102.D-0268 and 106. 216T.
† E-mail: s.kamann@ljmu.ac.uk new windows into a plethora of unexpected phenomena. Studies of their colour magnitude diagrams (CMDs) revealed the extended main sequence turn-off (MSTO) phenomenon (e.g., Mackey & Broby Nielsen 2007), high fractions of Be stars (Feast 1972;Grebel et al. 1992;Bastian et al. 2017;Milone et al. 2018), as well as split or dual main sequences (e.g., Milone et al. 2017). The characteristics of these features strongly depend on cluster age (e.g., Niederhofer et al. 2015), and the fact that some of them, such as the extended MSTO (Cordoni et al. 2018), have also been identified in lower-mass open clusters in the Milky Way, suggests that the underlying mechanisms work irrespective of cluster masses, metallicity or birth environments. While there are similarities to the CMDs of the ancient globular clusters, which host multiple populations, the origin of these complex features appears to be different between the old and young clusters.
In the ancient clusters, the complex CMDs are caused mainly by star-to-star chemical abundance variations (see reviews by Bastian & Lardo 2018;Gratton et al. 2019). Similar variations have been observed in massive Magellanic Cloud clusters with ages of ∼ 2 Gyr or above (e.g., Saracino et al. 2020;Martocchia et al. 2020;Cadelano et al. 2022), while they appear to be absent in the stars observed in even younger clusters (e.g., Mucciarelli et al. 2014). Instead, stellar rotation has been suggested as the dominant cause of the unexpected features in the young clusters (e.g., Bastian & de Mink 2009;D'Antona et al. 2015).
To date, most of the spectroscopic observational works on young clusters have been focused on the MSTO, due to the relative brightness of the stars there (e.g., Dupree et al. 2017). Studies of the extended MSTOs in Magellanic Cloud and Galactic open clusters have been able to directly link the spread in the CMD with the rotation rates of the stars (e.g., Marino et al. 2018b;Sun et al. 2019;Kamann et al. 2020). The cause of the split main sequence is less well constrained in such clusters. Marino et al. (2018a) measured the projected rotation rates ( sin ) of 31 stars along the blue and red main sequences in the young LMC cluster NGC 1818 (∼ 40 Myr) and found that the average sin values of the two sequences were approximately 70 km s −1 and 200 km s −1 , respectively.
The physical mechanism creating the different rotation rates is still under debate. The idea that all stars are born as fast rotators before tidal torques in binaries create the population of slow rotators (e.g., D'Antona et al. 2015) appears in conflict with observations of massive clusters in the Magellanic Clouds that have found comparable binary frequencies for fast and slowly rotating stars (Kamann et al. , 2021. Alternatively, Bastian et al. (2020) proposed that bimodal rotational velocity distributions originate from differences in the lifetimes of the stars' pre-main sequence disks while very recently, Wang et al. (2022) suggested stellar mergers as the cause of the blue main sequence.
In addition to the split main sequence observed in many young massive clusters (YMCs), such clusters also contain a high fraction of classical Be stars, i.e., rapidly rotating stars with a decretion disk surrounding them (e.g. Lee et al. 1991;Rivinius et al. 2013). Such stars can be identified photometrically or spectroscopically through the bright H emission from the disk. The high Be fractions of up to ∼ 50% (Bastian et al. 2017;Milone et al. 2018;Bodensteiner et al. 2020b) within YMCs with ages up to a few 100 Myr are consistent with expectations of a large rapidly rotating population of stars within the clusters as inferred through the split main sequences. Such a large population of Be stars, all with the same age and distance, offers the opportunity to study the phenomenon in more detail, and with greater ease, than possible in the field. For example, it allows us to study a special class of Be stars, known as shell stars, which are seen nearly edge-on and can be used to estimate the opening angle of the decretion disks.
Additionally, while there is general agreement that rapid rotation is a necessary condition for the Be star phenomenon, additional processes appear to be required for a star to become a Be star (e.g., Baade & Rivinius 2020). There are thought to be two main channels to form Be stars, the single star evolutionary path (e.g., Fabregat & Torrejón 2000) and a path through mass transfer binaries (e.g., Pols et al. 1991). The mechanism underlying the single star evolutionary path is the contraction of the stellar core during the main sequence lifetime, in combination with a radially outward transport of angular momentum to the expanding stellar envelope. This causes the star to approach its critical rotation velocity as it approaches the MSTO (e.g., Ekström et al. 2008;Granada et al. 2013;Hastings et al. 2020). At some point the critical velocity approaches the rotational velocity of the star and the formation of a Keplerian decretion disk is triggered. Alternatively, during the mass-transfer phase of a binary system, angular momentum is transferred as well, which can result in a rapidly rotating star (e.g., de Mink et al. 2013;Shao & Li 2014;Klement et al. 2019).
Observational studies have shown that binaries do play a significant role in creating Be stars. While optical spectroscopy struggles to confirm the binary nature of Be stars owing to the low masses and small sizes of their companions (e.g., Bodensteiner et al. 2020b, but also see Nazé et al. 2022), other methods, such as spectral energy distribution modelling of the disks (Klement et al. 2019), far-UV spectroscopy , the search for runaway Be stars (Boubert & Evans 2018), or near-infrared interferometry (Klement et al. 2022) have shown that a significant fraction of Be stars lives in binaries. Furthermore, spectroscopic observations of Be stars did not reveal the enhanced surface nitrogen abundance predicted for stars stemming from the single star evolutionary path (e.g., Lennon et al. 2005). For the YMCs in the Magellanic Clouds, Hastings et al. (2021) estimated that by adopting rather extreme assumptions regarding the initial binary population, non-canonical initial stellar mass functions, and the efficiencies of mass transfer, binaries can account for the large fractions of Be stars observed near the turn-offs of the clusters (e.g. Bastian et al. 2017;Milone et al. 2018). For a sample of Galatic open clusters, McSwain & Gies (2005) estimate that up to 73% of the observed Be stars originated in binaries, while the rest were contributed by the single evolutionary path. Models that invoke the single or binary star path to form Be stars make explicit predictions that can be tested with targeted observations (e.g., Granada et al. 2013).
Improving our understanding of the role of rotation (and Be stars in particular) is not only of interest regarding the morphology of cluster CMDs. Decretion disks provide interesting constraints for angular momentum transport in stellar models (Rímulo et al. 2018, e.g.). Furthermore, Be stars are known to occur in Be/X-ray binaries (Reig 2011), systems where a compact object (typically a neutron star) orbits a Be star and emits X-rays each time it passes through the disk. Such systems are of renewed relevance because they represent an intermediate step in proposed pathways for the formation of gravitational wave sources.
With the aim to better understand the stellar populations within YMCs, we have initiated a multi-epoch observational campaign targeting NGC 1850, a ∼ 100 Myr, ∼ 10 5 M cluster in the LMC. For this we use the Multi-Unit Spectroscopic Explorer (MUSE Bacon et al. 2010) instrument on the VLT, which has the power to study thousands of stars, from the main sequence through the red giant branch (RGB), in a single pointing, with high enough spectral resolution to find and study binaries through their radial velocity variations as well as derive the sin values of individual stars. Parts of this large dataset have already been used to determine the binary frequency on each of the arms of the split main sequence in order to test its origin (Kamann et al. 2021) as well as to study the peculiar system NGC 1850-BH1 El-Badry & Burdge 2022). In an independent analysis of the same data, Sollima et al. (2022) determined a dynamical mass of the cluster of 10 4.84 and found a link between oxygen abundances and line widths among MSTO stars that the authors interpreted as evidence for different stellar rotation rates. Here, we present the full data set and use it to study the rotational velocity distribution of the stars as well as the Be star population within the cluster. This paper is organised as follows. We describe the MUSE obser-vations and their reduction in Sect. 2, followed by a summary of the photometric and spectroscopic data analysis in Sect. 3. In Sect. 4, the measurement of the stellar parameters is described before the sin values are discussed in Sect. 5. Sect. 6 is devoted to the Be and shell stars discovered in the data. We conclude in Sect. 7.

OBSERVATIONS AND DATA REDUCTION
NGC 1850 was observed with MUSE (Bacon et al. 2010) at the ESO Very Large Telescope as part of programs 0102.D-0268 and 106.216T (PI: Bastian). The observations were carried out in the MUSE wide-field mode (WFM), which provides low-to mediumresolution spectroscopy ( ∼ 1700−3500) across a wide wavelength range ( = 4 800 − 9 300 Å) and over a continuous field of view of 1 × 1 , with a spatial sampling of 0.2 . We made use of the adaptive optics (AO) system of MUSE in order to improve the spatial resolution of the data and facilitate the extraction of single-star spectra.
The observations were carried out during 16 visits in 13 individual nights between March 2019 and February 2021. During each visit, the same two fields were observed, a central field, situated on the cluster centre and an outer field, located about 1 to the south-east of the cluster centre (cf. Fig. 1). The exposure times per visit were 2 × 400 s for the central field and 3 × 500 s for the outer field. In between the individual exposures, derotator offsets of 90 degrees and small spatial offsets of ≤ 0. 4 were applied in order to homogenise the data quality across the field of view.
The data reduction was performed using versions 2.6 and 2.8.3 of the standard MUSE pipeline (Weilbacher et al. 2020). All basic reduction recipes were carried out as explained in Kamann et al. (2020). However, owing to the strong nebular emission across NGC 1850, the sky subtraction procedure had to be optimized compared to our previous work. Briefly, under the skymethod=model setting that we adopted, the MUSE pipeline identifies the requested fraction of faintest spaxels in the field of view and combines their spectra. The result is split up into a continuum and a telluric emission-line component. While the former is subtracted globally from the data, the latter is fitted to the data of each slice individually, in order to account for variations in the line spread function (LSF) across the field of view (see Weilbacher et al. 2020, for details). In our case, however, the continuum component always contains a contribution from the nebular emission, which would be subtracted from the data when using the standard procedure. Therefore, we manually replaced the spectrum of the continuum component with zeros. This implied that our final data cubes still contain telluric continuum emission. This, however, does not affect the extraction of the spectra described in Sect. 3 below. The reason is that the telluric component is spatially flat across the field of view, and such a component can be readily accounted for when performing spectrum extraction via point spread function (PSF) fitting.
We created individual data cubes for each pointing and each visit. For illustration purposes, we further created one data cube combining all individual exposures from program 0102.D-0268. In Fig. 1, we show two colour images created from this combined data cube, using either the SDSS broadband filters or custom narrow-band filters centred in the O , N , and H emission lines. The latter illustrate the strong extended gas emission originating from the young cluster NGC 1850B, visible in the western corner of the northern (central) pointing.

Photometry
The HST data were taken in 2015 as part of programs 14069 (PI: Bastian) and 14174 (PI: Goudfrooĳ). PSF photometry was performed on the flat-field corrected, and bias-subtracted WFC3 images in the filters F275W, F336W, F343N, F438W, F656N, and F814W. This was carried out with DOLPHOT (Dolphin 2016), a modified version of HST (Dolphin 2000) using Jay Anderson's PSF library. For more details on the photometry we refer the reader to Gossage et al. (2019); Kamann et al. (2021) and references therein. Bright stars were saturated in almost all the long exposures available, hence to recover their magnitudes we made use of the shortest exposure (7 sec) in F814W. This photometric catalog was then used as a reference for the spectra extraction (see Sec. 3.3 below). A F336W − F438W vs F438W CMD generated from the HST photometry is shown in Fig. 2. Following our earlier work described in Kamann et al. (2021), we identified two samples of red and blue main sequence stars in the magnitude interval 18.0 < F438W < 20.0, where the two sequences are most clearly separated. The reader is referred to Kamann et al. (2021) regarding the details of this process. In essence, we assumed that the blue main sequence stars constitute a constant fraction of 0.2 of all stars at a given magnitude level. The dividing line between stars we consider to lie on the red or the blue main sequence is depicted as a solid blue line in the right panel of Fig. 2. Wang et al. (2022) recently reported that the ratio of red to blue main sequence stars in a sample of similarly aged massive clusters is magnitude dependent (see also Milone et al. 2018). However, visual inspection of Fig. 2 suggests that at least for the magnitude range considered in this work, assuming a constant ratio of blue to red main sequence stars does not result in a significant number of misclassified stars.
When comparing the radial distributions of the red and blue main sequence stars, we found that the former were centrally concentrated compared to the latter. A two-sided Kolmogorov-Smirnov test of the two distributions yielded a probability of 5 × 10 −12 that the two were drawn from the same parent sample. This finding is at odds with the analysis of Correnti et al. (2017), who did not find any differences in the concentrations of the two populations.
In order to clean our sample of red main sequence stars from photometric binaries, we imposed a second selection criterion, illustrated by the dashed blue line included in the right panel of Fig. 2. This second line was constructed by shifting the dividing line between the two main sequences by Δ( F438W − F336W ) = 0.1 + 0.025 × ( F438W − 18) and follows the drop in stellar density visible to the red of the red main sequence. Stars lying redwards of this line are considered to be in binary systems with two luminous companions. We note that this is a purely photometric selection of binaries. Stars showing radial velocity variations were treated as described in Sect. 3.3 below.
We note that photometric binaries originating from the blue main sequence would overlap in CMD space with the red main sequence. Depending on the mechanism adopted in order to explain the split main sequence, the fraction of binaries among blue main sequence stars is expected to be lower, higher, or comparable to the fraction among the red main sequence stars. In Kamann et al. (2021), we detected similar fractions of binaries for both main sequences. Given the relative low number of blue main sequence stars, we expect any contamination of the red main sequence from "photometrically migrating" blue main sequence binaries to be small.

Isochrones
In order to derive stellar parameters from the HST photometry, we compared the data to isochrones from the MIST database (Dotter 2016; Gossage et al. 2019). The MIST isochrones are derived from the Modules for Experiments in Stellar Astrophysics (MESA Paxton et al. 2011). Following Yang et al. (2018), we adopted an age of 100 Myr, a metallicity of [Fe/H] = −0.24, a distance modulus of 18.45, and an extinction of = 0.301. We verified that the isochrones for this set of parameters provided a good by-eye fit to the CMD shown in Fig. 2.
Of particular relevance for the present work is the treatment of stellar rotation in the isochrone models, which is detailed in Gossage et al. (2019). Rotation is parameterized via the parameter Ω/Ω crit , specifying the fractional angular velocity of a star relative to the critical value at the zero-age main sequence (ZAMS). In the MIST isochrones, Ω crit is defined as the limit where the centrifugal force equals the gravity of the star. The isochrones are available for 10 discrete steps of Ω/Ω crit , ranging from 0 to 0.9. Gravity darkening in the models is treated following Espinosa Lara & Rieutord (2011) and its effects are accounted for via a surface-averaged modification of the luminosity and effective temperature. In reality, gravity darkening causes a viewing-angle (i.e. inclination) dependence of the observed colours of a rotating star, which is not included in the models.
In Fig. 2, we compare the HST photometry to four isochrones with different rotation parameters. We note that in the magnitude range displaying the split main sequence, i.e. 18 < f438W < 20, some features of the CMD are not accurately represented by the isochrones. For example, the blue main sequence shows bluer colours than even the isochrone without rotation (i.e. Ω/Ω crit = 0) predicts. Recently, Wang et al. (2022) argued that stellar mergers can account for this colour offset, as they lead to a rejuvenation of the merger product, resulting in bluer colours. However, it is unclear if mergers during the first few Myr following the formation of a cluster, as advocated by Wang et al. (2022), can rejuvenate stars sufficiently to account for the observed colour shift. Furthermore, at higher rotation rates, the isochrones predict slightly different slopes of the main sequence compared to what is observed. These deviations likely stem from remaining uncertainties in the modelling of massive stars and the treatment of rotation in those models. The treatment of rotation varies between stellar evolution codes, highlighting the need for more observational constraints suited to verify or reject assumptions made in the models.

Spectroscopy
We extracted individual stellar spectra from the final MUSE data cubes using P M (Kamann et al. 2013). The code works by using a reference catalogue of sources to determine a MUSE PSF model and spaxel coordinates of the resolved sources as a function of wavelength. This information is subsequently used to optimally extract the spectra of the resolved stars from the data. The reference catalogue used for NGC 1850 was the same as discussed in Sect. 3.1 and the stars included in the extraction were selected based on their F814W magnitudes. The extraction was performed on each individual data cube. As de- In both panels, only stars with MUSE spectra are shown and colour-coded according to their average spectral signal-to-noise ratio per pixel. Green lines show the predictions from MIST isochrone models for different rotation rates, relative to the critical rotation rate Ω crit . In the right panel, the blue solid line indicates the adopted division into red and blue main sequence stars, and the blue dashed line shows the adopted division between red main sequence stars and photometric binary stars.
scribed in Saracino et al. (2022), the extracted spectra were thereafter analysed with S (Husser et al. 2016), a code which determines stellar parameters via full-spectrum fitting against a library of templates. As in our previous work on NGC 1850, we used the synthetic templates from the library presented in Allende Prieto et al. (2018) when running S .
The individual radial velocities derived from the S fits were used in Kamann et al. (2021) and Saracino et al. (2022) to study the binary properties of the stars in our sample. In this work, however, our aim is to obtain a single spectrum per star in our sample at maximum signal-to-noise (S/N). For this reason, we combined the spectra obtained for the individual visits on a star-by-star basis. Before combining them, each spectrum was corrected for its radial velocity as measured by S and thereby converted to restframe. This correction ensures that the widths of the spectral lines, which are used below to infer the stellar rotation of our sample stars, are not dominated by the orbital motions in binary systems. We sound a note of caution that our binary detection via radial velocity variations is biased against systems composed of equally luminous companions (e.g. Giesers et al. 2019;Bodensteiner et al. 2020a) with blended spectral lines. Such SB2 or double-line binaries can still result in artificially high sin measurements. We discuss their impact in Sect. 5.1 below.
Furthermore, we corrected the individual spectra for atmospheric absorption prior to the combining. This is possible because S fits the telluric absorption bands simultaneously with the stellar features, using an internal library of atmospheric spectra (see Husser et al. 2016, for details).
In the combining process, we weighted the individual spectra by their signal-to-noise ratio (S/N, measured per wavelength bin and averaged across the entire wavelength range). For each star, spectra with a S/N lower by a factor < 0.5 compared to the one with the highest S/N were discarded altogether, as they tended not to improve the quality of the combined spectrum. In addition, we discarded any spectra for which the results from the S analysis were deemed unreliable for any of the following reasons. (1) The fit was marked as unsuccessful by S .
(2) The S/N of the input spectrum as estimated by S was < 5. (3) The velocity determined by crosscorrelating the spectrum with its best fitting template deviated by > 3 from the actual S result.
The gaseous emission across NGC 1850 (cf. right panel of Fig. 1) poses a major challenge for the extraction of spectra. It varies over spatial scales that are comparable to the resolution of the MUSE data ( 0.5 ), hence it cannot be easily deblended from the stellar emission. We experimented with setting up a fine background grid (using distances down to 10 spaxels between the individual grid points) and including the flux of each background component in the extraction process. However, residuals from the strong nebular emission lines were still visible in some of the extracted spectra. Therefore, we designed an approach to account for any contamination by nebular lines during the analysis of the spectra. For each extracted stellar spectrum, we created a mask as follows. We first selected a comparison sample of 100 spectra of stars with similar photometric magnitudes as the star linked to the target spectrum. Then, we calculated the median absolute deviation (MAD) of the spectral fluxes in the comparison sample and divided the result by the square root of the median spectral fluxes in the comparison sample. The latter is done in order to account for the fact that the noise in spectral absorption lines is typically lower than in the continuum. Afterwards, we removed the continuum from the result of the previous step, using a median filter of 200 wavelength bins width. On the continuum-corrected MAD spectrum, we determined the 84th percentile in a rolling window of size 200 bins and flagged all bins as contaminated where the actual values exceeded the smoothed ones by a factor of > 3. Finally, we removed isolated masked bins by processing the mask with a minimum filter with a window size of 3 bins. The masks created this way were included in the spectrum analysis of the combined spectra using S . We show the distribution of our targets in an HST ( F336W − F438W , F438W ) colour-magnitude diagram in Fig. 2, where they are colour-coded according to their S/N, determined using the method of Stoehr et al. (2008). The low S/N ratios of the evolved stars visible to the top right of the CMD are very likely a shortcoming of the method when applied to low-resolution spectra of cool stars, as the numerous molecular bands are mistaken as noise. Given that the evolved stars only play a very minor role in the present work, we did not make an effort to correct for this effect. We also note that for some of the brightest stars no photometry in said filters was available, because the stars were saturated in the HST images. In such cases, we made use of the F606W and F814W magnitudes recovered from the extracted spectra. They were compared to one of the isochrones overplotted in Fig. 2, and the missing magnitudes were copied from the nearest data point in ( F606W − F814W , F606W ) space.

Field stars
In order to clean our sample from field stars, we utilised the radial velocities derived from the S fits. For each star, we averaged the velocity measurements obtained across the individual epochs, using inverse-variance weighting. Prior to this process, the uncertainties of the single-epoch velocities were calibrated using the method presented in Sec. 3.4 of Kamann et al. (2020). Velocities which were deemed unreliable according to the criteria outlined in Sec. 3.3 above were discarded before averaging the results. Following this process, our kinematic sample consists of 4 207 stars with available radial velocity measurements.
Cluster membership probabilities were determined under the assumption that the observed sample of stars can be described by a cluster population and a field population. Regarding the cluster population, we further made the assumptions that its surface density follows a King (1962) profile with the structural parameters obtained by Correnti et al. (2017) and that its velocity dispersion profile can be modelled by a Plummer (1911) profile with central dispersion 0 and scale radius 0 . For the field population, we adopted a Gaussian velocity field with mean velocity back and velocity dispersion back . We note that because we did not try to clean our sample from Milky Way foreground stars, the fit parameters obtained for the field population could be biased towards low mean velocities and high dispersion values.
For each star, a membership prior was calculated based on its distance to the cluster centre and the surface density profile we adopted. The coordinates for the cluster centre were taken from Milone et al. (2018). Then, we used (Foreman-Mackey et al. 2013), a P implementation of the affine-invariant Markov-chain Monte Carlo (MCMC) sampler presented by Goodman & Weare (2010), to determine the model parameters ( back , back , 0 , 0 , and the systemic cluster velocity 0 ) in a maximum likelihood approach. We used 200 walkers in the process and the chains were prop-agated for 500 steps each. Discarding the first 200 steps of each chain as burn-in, we found the following set of parameters to maximise the likelihood of the model given the MUSE radial velocities, −19 arcsec, back = 252.2 ± 0.8 km s −1 , and back = 20.2 ± 0.8 km s −1 . For each parameter, we adopted the 50th percentile of the distribution returned by the chains as best-fit parameter, while the confidence intervals were obtained from the 16th and 84th percentiles of the same chains. When comparing our results to Song et al. (2021), who performed a similar analysis on NGC 1850, we find reasonable agreement in the parameters that appear in both models (see their  Table 7).
Given our set of model parameters, we are able to assign posterior membership probabilities to all stars with radial velocity measurements available, using the method described in, e.g, Watkins et al. (2013). The distribution of posterior membership probabilities is clearly bimodal, with 11% of the sample having probabilities < 0.1 and 60% having probabilities > 0.75. Because the vast majority of the stars with F336W − F438W > 0.8 and F438W < 21 (see left panel of Fig. 2), which constitute the red giant branch of the LMC field population, fall into the former group, we decided to consider the 3 737 stars with membership probabilities > 0.1 in the subsequent analyses. We note that because the mean velocity of the cluster and field populations only differ by 5 km s −1 , any velocity-based separation into field and cluster stars remains somewhat uncertain.

NGC 1850B
As mentioned above, the MUSE footprint covers the young cluster NGC 1850B, visible in the right edge of the central pointing in Fig. 1. In order to identify stars belonging to this cluster, we drew a circle of 10 arcsec radius around the visually estimated cluster centre ( = 05 h 08 m 39.3 s , = −68 • 45 45. 5) and considered all stars within this circle as members of NGC 1850B. Out of the 903 stars in the HST photometry that were identified this way, MUSE spectra are available for 140 stars. In contrast to the field stars, we decided to keep the stars associated with NGC 1850B in our sample and discuss their impact on the analysis when appropriate. Given that our selection results in a mixture of stars from NGC 1850 and NGC 1850B, we did not try to perform a dedicated analysis of the stellar content of NGC 1850B.

STELLAR PARAMETERS
During the spectral analysis of the combined spectra with S , we measured the projected rotational surface velocity sin and the effective temperature eff , and the surface gravity log of every star. The initial values for the latter two were determined from the comparison between the HST photometry and the MIST isochrones as outlined in Sect. 3.1. In addition, for stars for which isochrone comparison suggested a value eff < 8 000K, we also included the metallicity [Fe/H] as a free parameter in the analysis. As stars above this temperature threshold do not show any significant metal lines in the MUSE spectra, the metallicity was fixed to the isochrone value for such stars. Note that the inclusion of log in the fitted parameters deviates from our analysis for NGC 1846 (c.f., Kamann et al. 2020), where this parameter was fixed to the value obtained from the isochrone comparison. While the main conclusions of our work are unaffected by this choice, we found that the sin values determined with varying log agreed better with the values obtained via individual line fits (cf. Sec. 5.1).
Following the analysis of the combined spectra, we applied several quality cuts to our sample. Besides discarding results from formally unsuccessful S fits, we applied a S/N cut at 20, and also dropped results from spectra for which the recovered 814 magnitude deviated strongly from the corresponding value in the HST catalogue. Because the latter could be a sign that the spectrum is contaminated by nearby stars (as a result of PSF mismatches or inaccuracies in the underlying HST astrometry), we discarded spectra for which the Mag Accuracy parameter used by P M was < 0.5. In combination with the cleaning for field stars as described in Sect. 3.4.1, these criteria resulted in a final sample of 2 184 stars with valid results that will be discussed in the following.
For eff , we find a large range of values, as expected given the large range in spectral types covered by our observations. It is worth noting that we measure a median temperature offset of 1 003 K between the samples of red and blue main sequence stars determined as outlined in Sec. 3.1. For the same range in magnitudes (18.0 < F438W < 20.0), the MIST isochrones predict temperature changes of up to 1 500 K when increasing sin from zero to almost critical.
The median metallicity determined from the cooler stars ( eff < 8 000, see above) is −0.33, with the 16th and 84th percentiles of the distribution being located at −0.40 and −0.20. Note that the scatter in our measurements does not imply an intrinsic metallicity spread in NGC 1850, but mostly reflects our measurement uncertainties. Recently, Song et al. (2021) measured a metallicity of [Fe/H] = −0.31 using high resolution spectroscopy, in good agreement with the value derived in this work. In addition, our value also agrees with the metallicity obtained by Sollima et al. (2022) in their analysis of the MUSE data (−0.31 ± 0.01).

STELLAR ROTATION
We studied the stellar rotation of the NGC 1850 stars in two ways. First, we looked at individual stellar spectra. This analysis, which is described in Sect. 5.1, itself rests on two pillars. On the one hand, the spectral fits mentioned in Sect. 3.3 above, which provided us with a value for the (Gaussian) line broadening required to match the spectra. We converted this value into a sin measurement as described in Appendix A. On the other hand, for stars that are hot enough to show He lines in their spectra, we also obtained sin measurements by directly fitting He and He lines covered by the MUSE spectral range. We note that while there are four He lines available (at 4 922 Å, 5 016 Å, 6 678 Å, and 7 065 Å), we only used the reddest two lines, as the bluer lines are blended with Fe lines (cf. Sec. 5.2). While in Be stars of early spectral type, the He lines often show emission-line components, we do not expect such complications for the cooler stars residing in NGC 1850. The only He line available is at 5 411 Å. However, only four stars in our sample are hot enough to show He absorption, and all of them belong to NGC 1850B according to the criterion of Sec. 3.4.2.
The second way in which we investigated stellar rotation was based on stacked spectra, obtained by summing up the MUSE spectra extracted for stars that are expected to be fast or slow rotators based on their positions in the HST colour-magnitude diagram. We present this analysis in Sect. 5.2 below.

Individual stars
As our parent sample for obtaining individual sin measurements, we considered the same set of 2 184 spectra mentioned in Sect. 4 above. For each spectrum in the parent sample, the fitting of individual lines was performed as outlined below. However, the S sample was cleaned from any Be stars or likely members of NGC 1850B. Be stars, which show hydrogen line emission that are likely to impact the spectral fitting, will be discussed in Sec. 6. NGC 1850B members were discarded because our isochrone fitting was not matched to the young age of this cluster. Excluding these two types of stars left us with a sample of 1 873 stars with sin measurements based on the S fits. Note that both Be stars and NGC 1850B members were not removed from the line-fitting sample, as their He lines may still result in useful sin measurements.
When performing the line fitting, each He line was fitted individually. Prior to the fits, the continuum of each spectrum was determined via a polynomial fit (from which the spectral lines were iteratively excluded using kappa-sigma clipping) and used to normalise the spectra. The fits utilise the line spread function (LSF) of MUSE at the wavelength of each line. This is possible because the MUSE LSF as a function of wavelength is measured by the pipeline and provided in the LSF_PROFILE calibration file (see Sect. 4.10 in Weilbacher et al. 2020). We performed a non-linear least squares fit in which the initial LSF profile was broadened with a kernel accounting for stellar rotation with a given sin . To calculate the latter, we used the rotBroad function available in the P A package 1 package (Czesla et al. 2019), which determines the impact of stellar rotation using the prescriptions provided in Gray (2008). The fits were carried out using (Newville et al. 2016) with a Levenberg-Marquardt optimisation and during each fit, we varied sin , the line centre c , and the total line flux . We only considered a line as successfully fitted if c was within 0.5 Å of the tabulated value and the value of exceeded the noise level of the continuum around the line by a factor of 7. If more than one line was successfully fitted per spectrum, the results from individual line fits were averaged. Furthermore, we used such cases to calibrate our measurement uncertainties. To this aim, the difference between each pair of sin values derived from the same spectrum was obtained and normalised by the squared sum of their uncertainties. In the case of correctly calibrated uncertainties, the resulting distribution should be Gaussian with a standard deviation of unity. Otherwise, the calibration is performed via multiplication of the uncertainties with the actual standard deviation of the distribution. We found that the uncertainties returned by underestimated the true uncertainties by a factor of 1.5.
We obtained a sample of 306 stars with line-based sin values. The considerably lower number of stars compared to the S approach can be explained by the gradual disappearance of the He lines at lower effective temperatures. As a consequence, the faintest stars for which we can still analyse the He lines are at F438W ∼ 18.7, corresponding to eff ∼ 13 000 K in the S fits and the isochrones.
Where sin measurements from both methods were available for the same stars, the two measurements were averaged, resulting in a final sample of 1 963 stars with individual sin results available. The remaining 221 stars from the parent sample were either Be or NGC 1850B stars for which the sin measurement via He lines failed. Where possible, we compared the results from the two methods and found that on average, the sin values derived using S where larger than those derived from the single-line fits by 16 km s −1 . The standard deviation between the results from the two methods is 50 km s −1 , which is in agreement with the expected accuracy of our sin measurements, taking into account the spec- tral resolution and wavelength coverage of MUSE as well as the added complication of the nebulosity impacting the S fits in the Balmer lines.
In Fig. 3, we show the measured sin values for stars within NGC 1850 in a F336W − F438W vs. F438W CMD, focusing on the main sequence (MS) and main sequence turn-off (MSTO). Overall, we observe that for a given magnitude there is a correlation between colour and sin , in the sense that redder stars rotate faster. To illustrate this, we show in Fig. 4 the colour dependence of our sin measurements for different magnitude bins. For this purpose, we define a colour distance or "pseudo-colour" Δ F336W, F438W which measures the horizontal difference of a star relative to the median colour of main sequence stars at a given F438W magnitude. Besides the individual values, we also show in the left panels of Fig. 4 the median sin values in pseudo-colour bins of 0.05 mag as brown diamonds. For each magnitude bin, we finally show histograms of the sin measurements in the right panels, separately for the bluest 20% of the stars and the remaining 80%. This division is motivated by the photometric analysis, which resulted in a fraction of 20% of blue main sequence stars (cf. Sec. 3.1).
The sin distributions shown in the upper two rows of Fig. 4 look quite similar to those obtained for the MSTO of the 1.5 Gyr old cluster NGC 419 in Kamann et al. (2020). They confirm that indeed, stellar rotation plays a dominant role in shaping the MSTO of YMCs.
It is interesting to note that while for the brightest bin ( F438W < 17), the relation between sin and pseudo-colour appears continuous, a drop starts to appear for the other magnitude bins, in the sense that the bluest ∼20% of the stars have substantially lower sin values on average than the remaining stars. In particular at fainter magnitudes (18 < F438W < 20), where the two main sequences can be most easily distinguished, a clear bimodality is visible, with the first three bins (with Δ F336W, F438W < −0.1) having median sin values of 100 − 120 km s −1 , while the following bins have median sin values of 190 − 220 km s −1 . This confirms previous suggestions that the split main sequence in young (< 300 Myr) clusters is primarily due to stellar rotation (e.g., Marino et al. 2018b). Based on the lower three panels of Fig. 4, we adopt median sin values of 110 ± 20 km s −1 and 210 ± 20 km s −1 for the blue and red main sequences, respectively. The uncertainties that we assign to the median values are based on the scatter between the values of the individual pseudo-colour bins and the mean difference between the line-based and S -based sin values reported above, which we consider as representative for the strengths of the systematic errors involved in our analysis.
As will be discussed further in Sec. 5.3, where we compare our measurements to stellar evolutionary models, the critical (breakup) velocity crit for the stars included in Fig. 4 is ∼ 450 km s −1 (±50 km s −1 , depending on mass). We note that a small fraction of the measurements (23 stars) shown in Fig.4 exceed this value. While some of these outliers can be explained by the limited accuracy of our measurements for the faintest stars in the sample (see error bars included in the left panels of Fig. 4), visual inspection of the spectra of some of these sources also reveals SB2 binaries (i.e. binaries with two luminous companions that both contribute lines to the combined spectrum). These spectra have multiple component absorption lines, meaning that the integrated profiles will be broader than expected for a single star, which results in the assignment of high sin values for these sources. Given the small number of stars with sin > crit , we did not make an effort to remove them from the computation of the median values included in Fig. 4.
There is still a possibility, however, that the sin distribution shown in Fig. 4 is skewed by SB2 binaries. However, only binary stars with a mass ratio ∼ 1 and orbital periods 100 d will produce combined spectra with a period-induced line broadening that is comparable to the observed sin values. Such binaries are expected to be rare, so that we do not expect a significant impact of SB2 binaries on our measured sin distribution. We note that for stars on the binary main sequence identified in Sec. 3.1, which roughly corresponds to the range in Δ F336W, F438W 0.05 in the lower two panels of Fig. 4, we do not observe a trend towards higher sin values that would hint towards an impact of binary orbital motions on the observed line profiles. On the contrary, we observe a decline in the median sin values (albeit with substantial scatter). This may be expected if the binary main sequence is mainly populated by binaries on relatively wide orbits ( 100 d) that have been (partially) braked by tidal interactions, as tidal interactions increase with increasing companion masses (and hence towards larger Δ F336W, F438W ). Indeed, the stars in binaries with orbits 500 d are expected to be slowly rotating due to tidal interactions (Abt & Boonyarak 2004). In that case, however, the individual binary components would be blue main sequence stars. Their combination in a binary would result in a redder colour and therefore push the binary towards the red main sequence, rather than onto the binary main sequence. Hence, it is not surprising that most of the data points with sin > crit in Fig. 4 have pseudo-colours similar to the red main-sequence stars. We note that in our analysis, we treated systems on the binary main sequence in the same way as all other identified cluster members. Given the expected small impact of their orbital motions on the observed line profiles, this approach appears justified in hindsight.
In Fig. 3, it is further visible that the region of the CMD populated by blue stragglers (with F336W − F438W −1 and F438W 18.5) hosts stars with a large range in sin values. However, care must be taken as this region is also populated by main sequence stars of NGC 1850B. We highlight members of the younger cluster (according to Sect. 3.4.2) via green circles in Fig. 3. When we omit the likely NGC 1850B members, we find that a majority of the blue stragglers are relatively slow rotators (with sin 150 km s −1 ), but that some blue stragglers with sin 200 km s −1 are also observed. The latter include a spectroscopically identified Be star, visible in the right panel of Fig. 3 at F336W − F438W = −1.07 and F438W = 17.1 (as well as several Be star candidates identified photometrically, cf. Sec. 6.1). Wang et al. (2020Wang et al. ( , 2022 predict that the blue straggler region should be predominantly populated by slowrotating merger products. In this scenario, the fast rotators might be considered as products of recent mergers that did not yet have time to spin down -and their number could be used to estimate spin down times. However, besides mergers, mass transfer provides an alternative pathway to blue straggler formation, in which case a fast rotating product is expected. Along these lines, it is interesting to note that the three blue stragglers with sin > 300 km s −1 all show evidence for radial velocity variations in the MUSE data. The binary properties of our sample will be the topic of a separate publication (Saracino et al., in prep.). Given that many Be stars have been found to reside in binary systems with stripped companions, the presence of Be stars among the blue stragglers of NGC 1850 seems unsurprising.
We will further discuss the implications of the results shown in Figs 3 and 4 below in Sect. 5.3, following a summary of the results derived from the combined spectra. The rotation rates of Be and shell stars, which are highlighted in the right panel of Fig. 3, will be discussed in Sec. 6.

Analysis of combined spectra
In total, our MUSE spectroscopic sample contains 337 blue main sequence stars and 988 red main sequence stars with S/N>10. 2   made further use of the large sample sizes by creating S/N-weighted mean spectra for both populations, with the aim of detecting the same He lines we were able to identify in the individual spectra of turn-offs stars, or any metallic lines that are sufficiently narrow to enable direct measurements of sin via single-line fits. Indeed, the combined spectra showed a number of such lines, namely the He lines at 4 922 Å, 5 016 Å, 6 678 Å, and 7 065 Å, as well as Si lines at 6 347 Å, 6 371 Å, and the O line at 7 774 Å. In order to verify whether these lines are truly isolated at the spectral resolution of MUSE in the eff range under consideration, we downloaded a spectrum of the star HD 196426 from the 2nd data release of the X-Shooter spectral library (XSL, Gonneau et al. 2020). Based on the spectral analysis performed by Arentsen et al. (2019), the star has similar stellar parameters compared to the average main sequence star entering our combined spectra. We found that the two He lines at 4 922 Å, 5 016 Å are blended with Fe lines (at 4 924 Å and 5 018, Å, respectively), while the O line at 7 774 Å is a doublet. Hence we are left with four truly isolated lines, for which the line profiles for the different spatial distributions of the two populations in combination with the complicated spectroscopic selection function (as S/N depends on both brightness and location of a star). the two combined spectra are compared in Fig. 5. In all cases, it is clearly visible that the line profile in the spectrum of the red main sequence stars is broader compared to the spectrum of the blue main sequence stars, confirming the results from the individual stars of a strong difference in sin between the two sequences. We note that the difference in effective temperature between the red and blue main sequence stars (∼ 1 000 K) are very unlikely to be responsible for the observed differences, given that its impact on the lines shapes is below the resolving power of MUSE and that, as we will show below, we obtain consistent results for different elements (He and Si ). In addition, note that the red main sequence stars are cooler than the blue main sequence stars. Therefore, temperature-dependent effects would counteract the observed differences. In order to quantify the line width differences visible in Fig. 5, we carried out single-line fits in a very similar fashion to those described for the individual spectra in Sec. 5.1. Again, we operated on normalized spectra and fitted for sin by convolving a model for the MUSE LSF with a broadening kernel generated using the rotBroad function available in the P A package (Czesla et al. 2019). The results are summarized in Table 1. The uncertainties included in Table 1 were obtained by repeating the fits for 100 bootstrap realizations, in which we randomly picked red and blue main sequence stars from our sample and averaged their spectra in the same manner as before.
The weighted averages of the results listed in Table 1 are 99 ± 5 km s −1 and 188 ± 7 km s −1 , respectively, for the blue and red main sequences. The results are in very good agreement with those obtained from the analyses of the individual spectra, giving us further confidence into the validity of our approach.
We performed another check in which we took the XSL spectrum of HD 196426, for which Arentsen et al. (2019) provide a low sin of 20 km s −1 , and converted it into a mock MUSE spectrum. This was achieved by convolving the XSL spectrum with the MUSE LSF model provided by the pipeline as a function of wavelength, adding Gaussian noise such that we achieved similar S/N as in the combined MUSE spectra, and finally rebinning the spectrum to a constant sampling of 1.25 Å per pixel. For the resulting mock spectrum, we performed the same line fits as for the two combined spectra. For all four lines listed in Table 1, we found low sin values 40 km s −1 , consistent with the detection threshold we expect in our data. This test shows that stellar rotation is indeed the dominant line broadening mechanism in the high sin regime we are concerned with.

Comparison to models
In Sec. 5.1 and 5.2, we found strong evidence for very different sin values along the blue and red main sequences of NGC 1850. In order to compare our results to the latest stellar evolutionary models, we need to correct the former for the effects of inclination. In this work, we assume a an isotropic distribution of spin axes, i.e. sin = /4. While it has been proposed that in star clusters the spin axes could be aligned as a result of cluster formation (Rey-Raposo & Read 2018), observational evidence for anisotropic spin distributions in clusters is still sparse (e.g., Lim et al. 2019;Healy et al. 2021). We did not make an effort to search for a possible deviation from isotropy in NGC 1850, but consider this a worthwhile endeveour for a future publication.
Under the assumption that the spin axes are distributed isotropically, our median values derived from the combined spectra correspond to mean equatorial velocities of surf, blue = 127 ± 7 km s −1 and surf, red = 233 ± 9 km s −1 . Using the values from the individual fits instead, we find slightly higher values of surf, blue = 140 ± 26 km s −1 and surf, red = 276 ± 26 km s −1 .
For the magnitude range 18.0 < F38W < 20.0 and the adopted cluster properties of NGC 1850, the MIST isochrones predict stellar masses between 2.3 and 4.2 M . The median predicted masses for our samples of red and blue main sequence stars are 3.01 M and 3.07 M , respectively. As mentioned above, the SYCLIST models published in Georgy et al. (2013) predict a critical velocity of crit = 450 km s −1 for a 3 star at the age and metallicity of NGC 1850 (cf. Fig. 6). Adopting this value, our median velocities correspond to ranges in / crit of 25 − 30% and 50 − 60% for the blue and red main sequences, respectively. To enable a better comparison with isochrone predictions, we also express our results in terms of the critical angular velocity, Ω crit . Note that the relationship between / crit and Ω/Ω crit is non-linear, given the deformation of a star as its spin increases (e.g., Granada et al. 2013). We find ratios of 35 − 40% and 67 − 79% for the median Ω/Ω crit of the two main sequences.
Our results suggest that rotation velocities close to the critical value 3 are not required in order to explain the split main sequence. At face value, this is in agreement with the latest predictions from isochrone fitting: The MIST models discussed in Sect. 3.2 suggest that blue main sequence stars have Ω/Ω crit 0.3, whereas red main sequence stars have Ω/Ω crit 0.6. The models presented by Wang et al. (2020) are instead parametrized in / crit and Wang et al. (2022) found that values of 0.35 and 0.65 match the observed main sequences in the young massive cluster NGC 1755. However, one must be careful in that both of the aforementioned models use the critical velocity at the zero-age main sequence (ZAMS). As shown by Hastings et al. (2020), / crit (and equivalently Ω/Ω crit ) can vary substantially over the main-sequence lifetime of a star. To illustrate this, we show in Fig. 6 the time evolution of the surface velocity, its critical value, and the ratio of two as predicted by the SYCLIST models by Georgy et al. (2013) for a 3 M star in NGC 1850 (see Hastings et al. 2020, for similar plots showing the predictions for the models used by Wang et al. 2022). Even though the model was initialized with a / crit = 0.9 at the ZAMS, its value at the age of NGC 1850 (100 Myr) has decreased to ∼ 0.7, as a result of a steep drop in surface velocity in the first 10 Myr.
In Fig. 6, it can be seen that apart from the initial drop, the SYCLIST models predict the surface velocity to barely change during the main sequence evolution (see also Bastian et al. 2020). On the other hand, a substantial decrease in the critical velocity, from an initial value ∼ 500 km s −1 at the zero-age main sequence to ∼ 320 km s −1 at the turn-off age of roughly 300 Myr, is predicted. 3 Our definition of "near critical" or "close to critical" is 80% of the break-up velocity This is caused by the expansion of the stellar envelope in response to a strong chemical gradient between the convective core and the radiative envelope. Comparing the value at the MSTO to the sin distribution of red main sequence stars in Fig. 4 suggests that a substantial fraction of the stars will be close to critically rotating when reaching the end of their main sequence lifetimes. Interestingly, our data do not show an increase in / crit when approaching the MSTO. For an age of 100 Myr and = 0.006, the crit predicted by the SYCLIST models for stars of 4 M and 5 M are ∼ 450 km s −1 and ∼ 400 km s −1 , respectively. As 5 M corresponds to the turn-off mass of NGC 1850, we can compare the latter value to the sin measurements shown in the top row of Fig. 4. Barely any values exceed 250 km s −1 , so that, unless the observed distribution is significantly compressed by inclination effects, the bulk of our MSTO sample is restricted to / crit 0.6. At face value, this could be taken as evidence for a lack of near-critically rotating stars at the MSTO of NGC 1850. However, the sin measurement for (almost) critically rotating stars can be biased towards lower values because of the rotationally induced darkening of the equatorial regions of the star (Townsend et al. 2004). In addition, NGC 1850 contains a large fraction of Be stars among its MSTO population (cf. Sec. 6), which are considered to be the outcome of near-critical stellar rotation.
Comparing the upper two rows in Fig. 4, we observe a shift towards lower sin values for the brightest stars. We found a similar trend in our previous work on the 1.5 Gyr old cluster NGC 1846 . This might indicate that the stars at the tip of the main sequence are already being braked as they start evolving into red giants. Fig. 3 shows that the trend that at a given magnitude, bluer stars rotate slower than redder stars is not restricted to the magnitude range of the split main sequence, but persists all the way to the MSTO. When comparing to the isochrone tracks shown in Fig. 2, it becomes evident that this sin dependency of the colour is only predicted for the fainter stars in our sample (i.e. in the magnitude range showing a split main sequence), yet not at the MSTO, where the isochrone tracks for different rotation rates merge. Interestingly, this is not the case in the models presented by Wang et al. (2022), where a sin dependency of the observed colour is predicted even for MSTO stars. This illustrates how our observations can be used to scrutinise predictions from stellar evolutionary models.
Overall, our observations appear to be in reasonable agreement with current stellar evolutionary model predictions. In particular, the predicted difference in the rotation rates of blue and red main sequence stars is confirmed. However, an in-depth comparison is hampered by the different reference times. While observations naturally reveal the stellar rotation properties at a given cluster age, models typically refer to these properties at the zero-age main sequence. This is particularly problematic as different models predict different relations between rotation at the zero-age main sequence and at a given cluster age. For example, the drop in visible at early ages in Fig. 6 appears to be absent in the models used by Hastings et al. (2020). Model predictions of the surface velocity as a function of cluster age will offer a promising venue for future research, enabling the combined analyses of photometry and sin measurements in order to better understand star cluster populations (e.g., Lipatov et al. 2022).

Be AND SHELL STARS
As established above, NGC 1850 hosts a substantial population of stars that are rapidly rotating and previous studies have found large populations of Be stars within the cluster using HST F656N narrowband photometry (Bastian et al. 2017;Milone et al. 2018). A common method to identify Be stars photometrically is by detecting outliers in a colour calculated from F656N , which is centred on H , and a nearby broadband filter. As our photometry is lacking any -band equivalent and using a bluer (redder) broadband filter could result in cool (hot) stars being misclassified as Be, we calculated the photometric index 438, 656, 814 = F438W − 2 × F656N + F814W and plotted it as a function of F438W (cf. Fig. 7). For the magnitude range 16 < F438W < 19, we defined the ridgeline of the main sequence (after removing all stars with 438, 656, 814 > 0.5) and selected as Be stars all sources that deviated from the ridge line by more than 6× their uncertainty in 438, 656, 814 . 4 This resulted in a sample of 397 Be star candidates. MUSE spectra are available for 218 of them, with most of the remaining candidates being located outside of the observed MUSE field of view.
In order to identify Be stars in the MUSE spectra, we extended the fitting of individual lines described in Sect. 5.1 to the H line, using a double Gaussian profile in order to account for a potential emission component. Following the fits, we summed up the equivalent widths of the two Gaussian components. In Fig. 7, we show the resulting H equivalent (EW H ) width as a function of the photometric index 438,656,814 , with the lower panel showing an almost linear relation between the spectroscopic and photometric indices. Fig. 7 demonstrates that there is good agreement between the stars showing H emission in the MUSE data and those showing F656N excess in the HST data. We measured EW H < 0 in the spectra of 202 stars, out of which 185 have also been flagged as Be star candidates using the photometric approach. For another 33 photometric Be star candidates, we measured EW H > 0. As can 4 lower thresholds would result in the misidentification of normal MSTO stars as Be stars, given the spread of the MSTO in F438W − F814W colour. be verified from Fig. 7, the latter stars have low F656N excesses, indicating that their emission components are not strong enough to fill in the entire H lines. This was confirmed by visual inspection of some of the spectra.
The relation between EW H and 438,656,814 shown in the bottom panel of Fig. 7 does show significant scatter. Part of it can be explained by the time gap between the HST and MUSE observations of about 5 years, as Be stars show variability on smaller timescales (Labadie- Bartz et al. 2018). In addition, some of the H line profiles of the spectra underlying Fig. 7 show contamination from the nebulosity in the field of view of the cluster (see Fig. 1). We note that the nebulosity varies on scales similar to the resolution of MUSE (∼ 0.5") and its H emission line extends across seven wavelength bins (corresponding to 9 Å or 400 km s −1 ). Hence, even following a PSF-based extraction of the stellar spectra, a fraction of the stellar H lines remained affected by the nebulosity. Fig. 7 suggests a trend of a decreasing fraction of Be stars with increasing F438W magnitude. This trend, which has been noted previously for NGC 1850 (Bastian et al. 2017) is unlikely to be caused by incompleteness due to decreasing emission line strength, even though the latter is generally expected to decrease for later spectral types given the lower ionizing fluxes. However, in the magnitude range F438W 18.3, the few Be stars that are observed still cover a similar range in F656N excess and H as the more numerous brighter Be stars. Also, we do not expect the diffuse nebulosity to play a significant role here because otherwise, we would expect to see a significant increase in the width of the main sequence in the top panel of Fig. 7 towards fainter magnitudes.

Be star demographics
In order to have a closer look at the fraction of Be stars as a function of magnitude, we revert to the photometric selection. This is because out of the 397 Be star candidates identified based on their F656N excess, spectroscopy is only available for 218 stars. Given the good agreement between the spectroscopic and photometric selection illustrated in Fig. 7, we do not expect any significant contamination of the Be star sample when using the latter.
In Fig. 8, we show the position of the Be stars in a ( F336W − F814W , F814W ) CMD. In the magnitude range 16 < F438W < 19, we constructed six selection boxes of 0.5 mag width and determined the ratio of Be stars to all other cluster members for each box. It can be seen from Fig. 8 that the fraction increases by an order of magnitude between F438W = 19 and F438W = 17. According to the MIST isochrones, the stellar mass increases from 3.2 M to 4.8 M in this range. We note that similar behaviour has been observed in a sample of LMC/SMC stellar clusters by Milone et al. (2018), where few or no Be stars were found on the main sequence but large fractions were found near the MSTO.
Such a trend with magnitude is in agreement with the model predictions of Hastings et al. (2020), who treat the Be phenomenon as being due primarily to single star evolution. In their models, the trend is explained by the dropping of the critical velocity of a star as it approaches the MSTO (cf. Fig. 6). Indeed, based on the results we obtained for the split main sequence, it seems plausible that a large fraction of the stars will be rotating close to their break-up velocities towards the end of their main sequence lifetimes. On the other hand, as discussed in Sec. 5.3, we do not find a significant number of near critical rotators among the current MSTO stars of NGC 1850. As shown in Sec. 6.4 below, this finding also applies to the Be stars. While on average, they rotate faster than normal MSTO stars, their sin distribution does not extend to significantly higher values. In order to better understand the efficiency of forming Be stars via the single-star channel, it would be insightful to compare the sin distributions of main sequence stars with comparable masses of 3 M across clusters of different ages (up to ∼300 Myr) and look at the evolution of / crit . A trend of Be star fraction with magnitude as shown in Fig. 8 can also be reproduced by models where Be stars are formed exclusively through binary channels (Hastings et al. 2021). However, such models require rather extreme assumptions regarding the initial binary fraction (near unity), a non-canonical initial mass function and very non-conservative mass transfer.
Note that the Be star fractions shown in Fig. 8 are based on the assumption that the Be star have comparable F438W magnitudes to normal B stars of the same mass. As shown by Haubois et al. (2012), the presence of a disk can alter the broadband magnitudes of such stars, to an amount that depends on wavelength, disk viscosity, and inclination angle. Haubois et al. (2012) found differences between +0.2 and -0.35 magnitudes for the band. As the impact of the disk increases with wavelength, we expect no strong F438W magnitude differences between B and Be stars in NGC 1850.
We further looked at the spatial clustering of the Be stars and checked if they were more or less centrally concentrated than normal B stars of the same magnitude range in. When performing a twosample Kolmogorov-Smirnov test on the distributions of projected distances to the cluster centre, we obtained a probability of 17% that the two were drawn from the same parent distribution. One might argue that if a large fraction of Be stars live in binaries, they should be centrally concentrated compared to normal B stars. However, NGC 1850 is only ∼150 Myr old, so that barely any evolutionary mass segregation is expected.

Shell stars
Visual inspection of some of the Be-star spectra revealed that they could roughly be divided into two classes, as illustrated in Fig. 9. Besides the "normal" Be-stars spectra, which are characterized by H emission, broad Paschen lines, and few He lines, we identified a group of spectra that showed a pronounced double peak in the H line profile, a series of very narrow Paschen lines extending to much higher orders than for the "normal" Be-star spectra, and a number of narrow Fe and Si absorption lines. These features are characteristic for shell stars, Be stars observed (almost) equator-on such that that the star is seen through the decretion disk (e.g., Rivinius et al. 2006).
We noted that the strengths of the various shell features (doublepeaked H line, high-order Paschen lines, Fe and Si lines) varied from star to star. This gradual transition can be understood in two ways. First, for decreasing inclinations, the fraction of the photosphere that is observed through the disk decreases. Second, the amount of matter in the disks is time-dependent. In order to obtain an unbiased estimate of the number of shell stars in our sample, we performed a cross-correlation of all spectroscopic Be stars against the two prototypical spectra shown in Fig. 9. Then, we compared the strengths of the two cross correlation signals using the parameter as defined by Tonry & Davis (1979) and classified as shell stars all stars for which shell > 5 and shell / Be > 1.5. The 47/202 stars classified this way are highlighted as green diamonds in Fig. 7. It is evident that compared to normal Be stars at the same magnitude, shell stars show higher H equivalent widths. Furthermore, the shell stars appear fainter in F438W than the normal Be stars on average. These trends appear to be in agreement with the expectations from absorption in a disk.
As mentioned previously, in our sample of 202 spectroscopic Be Fe ii Figure 9. Comparison of a "standard" Be star spectrum (star number 92) and a prototypical shell star spectrum (star number 5 000 003) from our sample. Vertical dashed lines represent absorption features used to identify the shell stars.
stars, we classified 47 as shell stars, giving a shell-star fraction of 23%. This fraction can be used to estimate the opening angle of the disk, under the assumption of a random distribution of spin axes (as found by Hummel et al. 1999, for another young massive cluster, NGC 330), resulting in a half-opening angle of = 13 degrees. The same value was derived by Hanuschik (1996), who used a sample of 114 Galactic Be (including 26 Shell stars) stars to estimate the shellstar fraction and obtained 22.8%. The good agreement suggests that there is not a strong effect of environment or metallicity on the Be star disk opening angle (at least in the comparison between the Milky Way field and NGC 1850 in the LMC). We note that, as shown by Cyr et al. (2015), the inferred opening angles depend on the diagnostic applied and the disk geometry adopted. The value derived for NGC 1850 is still within the confidence intervals provided by Cyr et al. (2015).

Implications for age spreads within clusters
The location of the shell stars in a ( F336W − F438W , F438W ) CMD is shown in Fig. 10. It can be seen that the shell stars lie redder and fainter than the nominal MSTO within NGC 1850. As discussed above, this is likely due to the extinction of the host star when seen nearly equator on. Given the position of the shell stars within the CMD of NGC 1850, it appears that the effect of self-extinction of Be stars can be larger than the change in colour and magnitude due to the effects of stellar rotation on the stars themselves.
To date, stellar models that include rotation do not account for the influence of the decretion disks for rapidly rotating stars. Furthermore, disks are a transient phenomenon and build up and dissipate on timescales comparable to the stellar activity cycles. As a consequence, Be stars have no fixed positions in a CMD, but instead follow loops with outlines that are mainly determined by their inclinations (de Wit et al. 2006;Haubois et al. 2012;Rímulo et al. 2018). Hence, the current positions of Be or shell stars in the observed CMDs are not expected to be matched by models, even those including rotational effects. The expectation is therefore that the width of the MSTO will be larger than predicted by models that include a wide range of rotational velocities. Correnti et al. (2017) have studied NGC 1850, comparing the observed HST CMDs with isochrones that include stellar rotation. These authors conclude that the MSTO width is larger than can be explained through stellar rotation alone, hence that an age spread of ∼ 35 Myr must be present within the cluster. However, the flux excesses from decretion disks, and in particular the presence of significant numbers of shell stars within the cluster, which are not represented in the stellar models, calls into question the need for a major age spread within the cluster 5 .
The same shell star phenomenon is also expected to happen for lower mass stars (down to ∼ 1.5 M ), although potentially at lower rates. It would be more difficult to find these lower mass stars with disks, because, due to their lower temperatures, they will not ionise their disks, i.e., the disks will not display emission lines. However, such disks can still produce shell absorption lines. In addition, the disks should be detectable through their mid-IR colours (e.g., Kenyon & Hartmann 1995), a promising avenue for future work when higher spatial resolution mid-IR observatories come online.

Rotation rates of Be and shell stars
In Fig. 11, we show the sin measurements obtained for the Be and shell stars in our sample and compare them to the results obtained for normal MSTO stars. As explained in Sec. 5.1, we discarded the S results for Be and shell stars, given the lack of any adequate models in the Allende Prieto et al. (2018) spectral library. To enable 5 We point the interested reader to the review of Bastian & Lardo (2018) for a discussion of the evidence against significant age spreads within massive clusters. a better comparison to the normal MSTO stars, we therefore show in Fig. 11 only results based on individual line fits. Fig. 11 shows that on average, Be stars rotate faster than normal main sequence stars, as expected. The tail towards low sin values observed for the Be stars is likely produced by stars observed at low inclination . Given that Be stars observed at low inclination should not show shell features, it is reassuring that no shell stars with sin < 160 km s −1 are observed.
It is interesting to note that the sin distributions of Be and shell stars shown in Fig. 11 do not extend to higher values compared to normal MSTO stars. This would have been expected if all near-critically rotating MSTO stars appear as Be stars. Our results support the idea that fast rotation alone is not a sufficient condition for a star to show a decretion disk, either because the disk is a transient feature or because additional mechanisms are needed to trigger disk formation (e.g., Baade & Rivinius 2020). Furthermore, a comparison of the histograms shown in Fig. 11 to the critical rotation velocity predicted for the MSTO of NGC 1850, crit ≈ 400 km s −1 (cf. Sec. 5.3) suggests that the Be stars are not rotating close to the expected break-up velocity. Recently, Dufton et al. (2022) found an average ratio of / crit = 0.68 for a sample of (more massive) Be stars in the 30 Doradus region, while Rivinius et al. (2013) found an average value of ∼ 0.7 for the Galactic field when compiling various literature studies. While the distributions shown in Fig. 11 also extend to / crit ∼ 0.0.7, they suggest a somewhat lower average ratio in NGC 1850.

CONCLUSIONS
We have presented an analysis of multi-epoch MUSE spectroscopy of the stellar populations within the ∼100 Myr LMC star cluster, NGC 1850. Our full sample consists of more than 4000 stars, and each star has been observed for a number of typically 16 epochs. For the present work, we have combined the single-epoch spectra on a star-by-star basis and, following cuts based on S/N and cluster membership probabilities, ended up with a sample of 2 184 stars with MUSE spectra that we analysed in order to understand the distribution of stellar rotation across the stellar population of NGC 1850. The set of combined MUSE spectra, along with the HST photometry, is made publicly available 6 .
The main results can be summarised as follows: • There is a clear correlation between the colour of stars on the MSTO and their observed sin value. This confirms previous suggestions that the extended MSTO phenomenon is driven by the stellar rotation distribution (Bastian & de Mink 2009).
• The two branches of the split main sequence have markedly different sin distributions, with the blue arm made up primarily of slow rotators while the red arm consists mainly of rapid rotators. Similar bimodal distributions have also been found for B stars in the Galactic field (e.g., Sun et al. 2021) or the 30 Doradus region of the LMC (Dufton et al. 2013). Our result confirms previous results found for other clusters, based on much smaller stellar samples and confirms predictions of stellar models that include rotation (e.g., D'Antona et al. 2015).
• Over a relatively large range in magnitude (17 < F438W < 20), we find that the fast rotating stars have a median sin of ∼ 200 km s −1 and that the sin distribution is barely populated beyond ∼ 300 km s −1 . Comparison to stellar models suggests that not many main sequence stars are rotating close to their break up velocities. However, the evolution of the latter for stars approaching the MSTO suggests that a large fraction of fast rotators could be critically rotating at the end of their main sequence lifetimes. Interestingly though, the MSTO of NGC 1850 shows a lack of stars with sin values close to the predicted critical value of 400 km s −1 .
• Similar to previous analyses (e.g. Milone et al. 2018), we find that the Be star fraction is a strong function of magnitude, increasing towards the MSTO, and going to zero on the nominal main sequence. In light of the predicted drop in critical velocity depicted in Fig. 6, this appears in qualitative agreement with Be model predictions based on single star evolution (see also Hastings et al. 2020). However, under certain assumptions, a similar trend can also be explained in a binary driven scenario (see the discussion in Hastings et al. 2021). A promising avenue for future research is the comparison of rotation velocities of main sequence stars across LMC clusters of different ages, as they provide snapshots of very similar stars (in terms of mass, metallicity, and environment) at different timestamps of their main sequence evolution.
• Within the Be star population, we find 23% of the stars to be "shell stars", i.e., Be stars that are seen nearly equator-on. By assuming an isotropic distribution of the Be stars' spin axes, we can translate this fraction to the Be star decretion disk opening angle, for which we find ∼ 15 degrees (half opening angle), in agreement to what has been found for Galactic Be stars.
• The shell stars are found almost exclusively on the red side of the MSTO, suggesting that they are self-extincted by their own disks. The photometric shifts are significantly larger than those predicted for "normal" Be stars owing to the growth and waning of their disks with time. Since decretion disks are not included in stellar models (even those that include rotation), the spread on the MSTO will never be matched by stellar models alone. This may result in the spurious inference of age spreads within clusters that host Be stars.
• This phenomenon is also likely to be present in older clusters (with ages up to ∼ 1.5−2.0 Gyr), although the stars with disks will be more difficult to identify (since they will not be hot enough to ionise their disks to be observed as late type A or F 'e' (emission) stars). Such stars should be identifiable, however, through their mid-IR colours, via shell-absorption lines or filled-in Balmer lines (Slettebak 1982).
Studying LMC clusters of different ages again appears as a promising avenue for future research, in particular given that the Balmer lines are much easier accessible in clusters other than NGC 1850, where the strong nebulosity complicates matters.