How stellar rotation shapes the colour−magnitude diagram of the massive intermediate-age star cluster NGC 1846

ABSTRACT

We present a detailed study of stellar rotation in the massive 1.5 Gyr old cluster NGC 1846 in the Large Magellanic Cloud. Similar to other clusters at this age, NGC 1846 shows an extended main-sequence turn-off (eMSTO), and previous photometric studies have suggested it could be bimodal. In this study, we use MUSE integral-field spectroscopy to measure the projected rotational velocities (vsin i) of around |$1400$| stars across the eMSTO and along the upper main sequence of NGC 1846. We measure vsin i values up to |$\sim 250\, {\rm km\, s^{-1}}$| and find a clear relation between the vsin i of a star and its location across the eMSTO. Closer inspection of the distribution of rotation rates reveals evidence for a bimodal distribution, with the fast rotators centred around |$v\sin i=140\, {\rm km\, s^{-1}}$| and the slow rotators centred around |$v\sin i=60\, {\rm km\, s^{-1}}$|⁠. We further observe a lack of fast rotating stars along the photometric binary sequence of NGC 1846, confirming results from the field that suggest that tidal interactions in binary systems can spin-down stars. However, we do not detect a significant difference in the binary fractions of the fast and slowly rotating sub-populations. Finally, we report on the serendipitous discovery of a planetary nebula associated with NGC 1846.

Hertzsprung–Russell and colour–magnitude diagrams, stars: rotation, galaxies: star clusters: individual: NGC 1846

1 INTRODUCTION

It is now well established that the resolved colour–magnitude diagrams (CMDs) of young and intermediate-age stellar clusters (i.e. <2 Gyr) display features that are not well reproduced through simple stellar isochrones, such as extended main sequence turn-offs (eMSTOs, e.g. Mackey & Broby Nielsen 2007) and split or dual main sequences (e.g. Milone et al. 2015). These features have been suggested to be due to distributions of stellar rotational velocities (e.g. Bastian & de Mink 2009; Brandt & Huang 2015; Gossage et al. 2019). Rotation changes the internal structure of the star, because the centrifugal support and extra mixing in the core region alter its hydrostatic equilibrium compared to that of a non-rotating star of the same mass and composition. These effects modify the evolutionary path of the rotating star in the CMD relative to the equivalent non-rotating star (e.g. Meynet & Maeder 2000). Furthermore, in stars rotating close to their break-up velocity, the spherical symmetry is broken, leading to effective temperature variations across its surface (von Zeipel 1924), so that the observed magnitudes also become a function of the inclination of the rotation axis towards the line of sight. Finally, the additional mixing can change the surface abundances of a star compared to the non-rotating case. However, the relative impact of these effects is still a matter of debate and varies across stellar evolutionary models. Further, we note that rotation is not the only possible process affecting the hydrostatic equilibrium of a star (see Johnston et al. 2019).

Several studies argued that the impact of stellar rotation on the CMDs would be minor (e.g. Girardi, Eggenberger & Miglio 2011; Platais et al. 2012) and the other factors such as age spreads or differential extinction would be responsible for the observed features. However, recent observations have directly linked the position of stars in the CMDs to their projected rotational velocity (vsin i) values (e.g. Dupree et al. 2017; Bastian et al. 2018; Kamann et al. 2018b; Marino et al. 2018).

In order to explain features such as dual main sequences, it has been suggested that the rotational distribution of stars in young clusters may also be bimodal, with peaks at near zero and near critical (i.e. break-up) velocity (e.g. D’Antona et al. 2015). Evidence for the presence of stars rotating close to their break-up velocities comes from the discovery of large populations of Be stars in young clusters (Bastian et al. 2017). While rotational velocities of field stars with similar masses and spectral types show an extended and perhaps bimodal distribution (see Royer, Zorec & Gómez 2007, for early A and late B stars), this is not extreme enough to match what is needed to explain the observed oddities in the CMDs of young massive clusters. In addition, Royer et al. (2007) found that the bimodality in the distribution of rotational velocities among Galactic field stars disappears for intermediate A-stars and later spectral types, where, aside from members of binary systems, only fast rotating stars are observed. This raises the question of what causes the large range of rotational velocities required to explain the eMSTOs of intermediate-age clusters. D’Antona et al. (2015) suggested that binary stars could play a role. Braking via tidal interactions appears to work efficiently in field binary stars, which Abt & Boonyarak (2004) found to be predominantly slow rotators. However, in some clusters, the slowly rotating stars appear to be more numerous than the fast rotating ones (e.g. Milone et al. 2015). As braking only works efficiently in hard binaries that are expected to survive in the clusters for a Hubble time, this would imply that the young and intermediate-age clusters have fundamentally different binary properties than the Galactic globular clusters, which have binary fractions of typically a few per cent only (e.g. Milone et al. 2012). Alternatively, a strong dependence of the properties of binaries on the masses of their constituent stars would be required, given the difference in stellar mass between the stars observed in young and intermediate-age clusters and in Galactic globulars.

Until now, the sample sizes of stars with measured vsin i values have been limited, meaning that it has not been possible to study statistically representative samples of stars in order to meaningfully compare with models that include physically motivated distributions of rotational velocities. Thanks to their large field of view and the excellent throughput, modern integral-field spectrographs such as MUSE enable us to significantly increase the spectroscopic sample sizes in star clusters compared to traditional (multi-object) spectrographs (e.g. Kamann et al. 2018a). In an initial study presented in Kamann et al. (2018b), we used adaptive optics enhanced VLT/MUSE observations of the intermediate age (∼1.5 Gyr) cluster NGC 419 in the Small Magellanic Cloud (SMC). We found a clear decrease in effective temperatures across the eMSTO, in agreement with expectations of stellar models that include rapid rotation. Due to the somewhat low S/N of the data, we were unable to measure the vsin i values of individual stars on the eMSTO in a reliable way. Instead, we co-added the spectra of stars on the blue and red sides of the eMSTO and carried out a differential analysis of the resulting spectra. We found, in agreement with rotating stellar models, that the blue spectrum had a warmer effective temperature and lower rotational velocity than the red spectrum. Here we expand this initial study to another intermediate-age cluster with significantly higher S/N data, allowing us to determine the stellar parameters and vsin i values of large numbers of individual stars across the eMSTO.

In this work, we focus on the massive cluster NGC 1846 in the Large Magellanic Cloud (LMC). The cluster has an age of about 1.5 Gyr (e.g. Mackey et al. 2008) and a mass of about 2 × 10⁵ M_⊙ (e.g. Goudfrooij et al. 2014). The cluster hosts an extended main sequence turn-off, which shows some evidence for being bi-modal (Mackey et al. 2008). However, there is no sign of a dual main sequence, which is typically only seen in the CMDs of younger clusters.

This paper is organized as follows. We present our data in Section 2 and their analysis in Section 3. The results are described in Section 4 and discussed in Section 5 before we conclude in Section 6. In the reduced data, we further detected a previously unknown planetary nebula. This discovery is presented in Appendix  B.

2 DATA

We observed NGC 1846 with the Multi-Unit Spectroscopic Explorer (MUSE), a panoramic integral-field spectrograph at the ESO Very Large Telescope (Bacon et al. 2010). The observations were carried out during four nights, 2018 October 01, 2019- March 11, 2019-rch 14, and 2019 August 22. Each observation consisted of three exposures of 880 s integration time each, targeting the central 1 arcmin × 1 arcmin of NGC 1846. To homogenize the image quality across the field of view, we applied a spatial dither pattern and derotator offsets of 90° between the exposures. All observations were performed with the GALACSI adaptive optics system (La Penna et al. 2016), in order to correct the ground-layer turbulence and enhance the spatial resolution of the data. The effective seeing in the SDSS r-band was measured as 0.5, 0.8, 0.8, and 0.7 arcsec, respectively, for the four nights.

The data reduction was performed with the standard MUSE pipeline (Weilbacher et al. 2012, 2014). It performs all the basic reduction steps – bias removal, slice tracing, wavelength calibration, and flat-fielding – for each individual integral-field unit (IFU). Afterwards, the data from all 24 IFUs are combined, flux-calibrated, and corrected for geometrical distortions. The final result is a data cube with a spatial sampling of 0.2 arcsec, and covering the wavelength range from 480 to |$930\, {\rm nm}$| with a spectral resolution of |$\sim 0.25\, {\rm nm}$| full width at half-maximum (FWHM, corresponding to |$R\sim 1700\!-\!3500$|⁠). We created one data cube from the entire set of available exposures as well as four cubes containing only the exposures taken during a single night. The latter was done in order to search for binary stars via radial velocity variations. In Fig. 1, we show a gri colour image created from the cube with maximum depth. The effective seeing of this cube was measured as 0. 8 arcsec at |$550\, {\rm nm}$| and 0.6 arcsec at |$850\, {\rm nm}$|⁠.

Figure 1.

gri colour image of NGC 1846 created from the MUSE data presented in this work. The image size is |$1\, {\rm arcmin}^2$|⁠, North is up and East is left.

3 ANALYSIS

3.1 Spectrum extraction

We used pampelmuse (Kamann, Wisotzki & Roth 2013) to extract stellar spectra from the reduced MUSE cubes of NGC 1846. The code deblends the spectra of the resolved stars based on a wavelength-dependent model of the point spread function (PSF) that is recovered from the integral-field data and the source coordinates provided in an astrometric reference catalogue. For the latter, we used the Hubble Space Telescope (HST) photometry presented in Martocchia et al. (2018). It consists of Advanced Camera for Surveys (ACS) magnitudes in three filters, F435W, F555W, and F814W. For the extraction, we relied on the sources with measured F814W magnitudes, as this passband has the largest overlap with the MUSE wavelength range.

The extraction from the data cube at full depth yielded 5044 spectra. After removing spectra with low signal-to-noise (S/N < 5, measured per pixel and averaged over the MUSE wavelength range) and all stars for which the recovered F814W magnitude deviated significantly from the original one, we were left with a sample of 3638 spectra that we consider for further analysis. In Fig. 2, we show the location of the extracted spectra in the HST colour–magnitude diagram of NGC 1846. We are able to extract spectra down to a magnitude of m_F555W ∼ 23, approximately 2.5 mag below the eMSTO of NGC 1846. Spectra of eMSTO stars typically have S/N ∼ 25 in the MUSE data.

Figure 2.

HST colour–magnitude diagram of NGC 1846 (grey dots), produced from the photometry presented in Section 3.1. The stars for which useful spectra were extracted from the MUSE data are colour-coded according to the spectral signal-to-noise per pixel, averaged over the MUSE wavelength range. The green lines show MIST isochrone predictions for an age of log (t/yr) = 9.14, a metallicity of [Fe/H] = −0.37 (Z = 0.006), and different stellar rotation rates |$\, \Omega /\Omega _{\rm crit,\, ZAMS}$|⁠. See Section 3.2 for the definition of the critical rotation rate, |$\Omega _{\rm crit,\, ZAMS}$|⁠.

We repeated the extraction process for the four cubes created from the data taken during the individual nights. Applying the same quality cuts as before, we are left with 3331 spectra for night 2018 October 01, 2091 spectra for night 2019 March 11, 2162 spectra for night 2019 March 14, and 2374 for night 2019 August 22. The different numbers can be explained by the different seeing conditions during the individual nights.

3.2 Isochrone fitting

The isochrones we used to fit the aforementioned HST photometry of NGC 1846 are the same as those described in Gossage et al. (2019). Briefly, these models are an expanded version of those of MIST (Choi et al. 2016). The main difference is that these models contain a greater range of initial rotation rates, going from |$\, \Omega /\Omega _{\rm crit,\, ZAMS}=0.0$| to 0.9, in 0.1 dex steps, as well as a tailored stellar mass and metallicity range. The boundary conditions for these models are set by ATLAS12, while SYNTHE is used for the bolometric corrections (Kurucz 1970, 1993). These models are evolved until the end of core helium burning.

The rotation rate |$\, \Omega /\Omega _{\rm crit,\, ZAMS}$| is the equatorial angular velocity of the star, divided by its critical value as defined by Choi et al. (2016), at the zero-age main sequence (ZAMS). The critical value, |$\Omega _{\rm {crit,\, ZAMS}}$|⁠, is the velocity where centrifugal forces equal the surface gravity and is an intrinsic property of the star. Rotation is initiated at the ZAMS for these models, with the specified value (e.g. |$\, \Omega /\Omega _{\rm crit,\, ZAMS}=0.9$|⁠). The effects of gravity darkening (von Zeipel 1924) are included, following Espinosa Lara & Rieutord (2011). Unless a viewing angle is chosen, the models assume a surface averaged effect on the luminosity and temperature of the star due to gravity darkening. No viewing angle was assumed for the comparisons carried out in this work.

Furthermore, at cluster ages around |$1.5\,$|Gyr, turn-off stars have low enough mass to possess surface convection layers, leading to surface magnetic fields. These magnetic fields can interact with stellar winds, magnetically braking the rotation of the star. Our models do not include an explicit treatment of magnetic braking. Instead, the rotation rate is simply reduced according to a mass-dependent factor between 1.8 and |$1.3\, {\rm M_\odot }$|⁠, going from fully to non-rotating in that mass range, respectively. Realistic modelling of magnetic braking is being investigated (e.g. see recent implementations in Amard et al. 2016; Georgy et al. 2019). An implementation is planned for future generations of the MIST models.

The isochrone fits to the HST data were performed in the (m_F435W − m_F814W, m_F555W) CMD shown in Fig. 2. We started from the values provided in Goudfrooij et al. (2018) and adapted them until we achieved a good by-eye fit. Two parameter sets were found to provide a good fit. The first set is similar to the one determined by Goudfrooij et al. (2018), Z = 0.008 ([Fe/H] = −0.25) and A_V = 0.06, whereas the second set has lower metallicity, Z = 0.006 ([Fe/H] = −0.37), and higher extinction, A_V = 0.26. As cluster age, we assumed log (age/yr) = 9.14 in both cases. The following analyses are based on the isochrones with a metallicity of Z = 0.006. We verified that our results do not sensitively depend on the choice of isochrone (unless noted otherwise).

We show a comparison between the HST photometry and the isochrone predictions for different |$\, \Omega /\Omega _{\rm crit,\, ZAMS}$| ratios in Fig. 2. It confirms that the eMSTO of NGC 1846 can be modelled using a range of stellar rotation velocities. We observe a slight colour offset along the red giant branch (RGB) of the cluster, where the isochrone predictions are offset by Δ(m_F435W − m_F814W) ∼ 0.1−0.2 to the red from the data. This can be explained by an inconsistent treatment of convection in the stellar atmospheres (ATLAS; Kurucz 1970) and interiors (MESA; Paxton et al. 2011), which leads to a RGB effective temperature (T_eff) scale that is, in general, cooler than observed. Based on the results obtained by Choi et al. (2018), the next generation of MIST models will feature a self-consistent treatment of convection in MESA and ATLAS, leading to a larger convective mixing length parameter, α_MLT (see Choi et al. 2016, section 3.6.1 for a brief description of the mixing length theory in these models) and, thus, better agreement with observations.

We assigned each star with an available MUSE spectrum an effective temperature T_eff,iso and a surface gravity log g_iso via comparison with the final set of isochrones. To this aim, we determined the closest isochrone point to each star in (m_F435W − m_F814W, m_F555W) space and copied its values for the effective temperature and surface gravity. A small subsample of our stars is missing either the F435W or the F555W magnitude in the HST catalogue. In such cases, we assigned fixed values of |$T_{\rm eff, iso}=4\, 500\, {\rm K}$| and log g = 2.5. These values were used as initial guesses for determining the stellar parameters using full spectrum fitting (see Section 3.4 below).

3.3 Cross-correlation

We obtained initial radial velocities for the spectra extracted from each MUSE cube by cross-correlating them against synthetic PHOENIX spectra from the göttingen spectral library (Husser et al. 2013). For each spectrum, a matching template was selected using the stellar parameters from the isochrone comparison explained in Section 3.2.

The analysis was performed with a custom python programme inspired by the fxcor routine available in iraf. Prior to the cross-correlation, the continuum was subtracted from both the object and the template spectra and wavelength regions affected by telluric absorption were set to zero. Then both spectra were resampled to the same equidistant sampling in logarithmic wavelength space. The final step before the cross-correlation was a filtering in Fourier space, removing both the lowest and highest wavenumbers.

Each cross-correlation result was assessed using the method of Tonry & Davis (1979). In particular, we calculated the r′ parameter which scales with the ratio of the heights of the selected peak and an average noise peak in the asymmetric part of the cross-correlation function. As described in Kamann et al. (2014), r′ can be used to estimate the uncertainties of the radial velocity measurements. There is a tendency for this method to overestimate the true uncertainties. However, the radial velocities determined via cross-correlation were only used as initial guesses for the full-spectrum fitting described in Section 3.4 below. Hence we did not perform a calibration of their uncertainties.

We selected trustworthy results by making a simultaneous cut in r′ and S/N, namely r′ > 4 and S/N > 5. This resulted in a set of |$2\, 970$| radial velocities derived for the full-depth data set.

3.4 Full-spectrum fitting

The final step of the spectrum analysis consisted in a full-spectrum fit of each extracted spectrum. This was done using the spexxy tool which is described in Husser et al. (2016). It obtains a set of best-fitting stellar parameters for each spectrum by fitting it against the full göttingen spectral library of PHOENIX templates. spexxy comes with an interpolator that predicts the template spectrum at each point of the parameter space. It accounts for continuum mismatches between an observed spectrum and a template by a multiplicative polynomial that is applied to the template. Furthermore, spexxy is also able to account for the telluric absorption bands, using a dedicated library of atmospheric spectra that are fitted simultaneously with the PHOENIX templates. For more information on the individual features, we refer to Husser et al. (2016).

For each spectrum extracted from the data cube at full depth, we aimed to derive the effective temperature T_eff, the solar-scaled metallicity [M/H], the radial velocity v_los, and the line broadening along the line of sight, σ_los. Recall that these quantities are derived from a comparison between the template spectra and the entire observed spectrum, not individual lines. Hence, they are tied to the strong features in the MUSE wavelength range that are most sensitive to them, such as the H α, H β, and Paschen lines in case of T_eff, or the |${\rm Fe}\, \small {{I,II}}$| lines bluewards of |$600\, {\rm nm}$| as well as the Mgb and |${{\rm Ca}\, \small {II}}$| triplet lines in case of [M/H]. All of these features also contribute to measuring v_los and σ_los. However, owing to their narrower intrinsic widths, the discriminatory power of the metallic lines is enhanced compared to the hydrogen lines.

As initial guesses for the parameters, we used the results from the isochrone comparison described in Section 3.2 (for T_eff and [M/H]) and from the cross-correlation described in Section 3.3 (for v_los). In cases where the cross-correlation did not yield a reliable radial velocity, we used the systemic velocity of NGC 1846 (⁠|$239\, {\rm km\, s^{-1}}$|⁠, Mackey et al. 2013) as initial guess. A constant initial guess of |$30\, {\rm km\, s^{-1}}$| was used for σ_los. Because of the challenges involved in determining the surface gravity log g from low-resolution spectroscopy (see Husser et al. 2016), we fixed log g to the values derived from the isochrone comparison during the analysis. The analysis formally converged for 3189 out of the 3638 useful extracted spectra.

We adopted several criteria for excluding potentially unreliable results from the converged analyses. Unless otherwise noted, we only considered results stemming from spectra extracted with S/N > 10. In addition, we removed stars for which the analysis yielded |$\sigma _{\rm los}\lt 5\, {\rm km\, s^{-1}}$|⁠. We found that in this regime, the uncertainties returned for the σ_los and v_los measurements were sometimes unreasonably high, suggesting that the optimization got stuck at the boundary at σ_los = 0. Note that this does not imply that we can reliably measure a line broadening down to |$\sigma _{\rm los}=5\, {\rm km\, s^{-1}}$|⁠. Our sensitivity limit in this respect is higher (see discussion in Section 4.1).

In Figs 3 and 4, we show the results obtained for T_eff, [M/H], and the projected rotational velocity vsin i. Only stars that are likely members of NGC 1846 (cf. Section 3.5) and passed the aforementioned reliability cuts are shown. As mentioned earlier, instead of vsin i, spexxy parametrizes the line broadening via the standard deviation of the Gaussian broadening kernel, σ_los. To convert our results, we used dedicated simulations that are detailed in Appendix  A. We will discuss the results displayed in Figs 3 and 4 in detail in Section 4 below.

Figure 3.

The distributions of projected rotational velocity (left), effective temperature (centre), and scaled-solar metallicity (right), across the colour–magnitude diagram of NGC 1846. Only results derived from spectra extracted from the MUSE data with S/N > 10 are displayed. The results were obtained in full-spectrum fits as outlined in Section 3.4. In the centre panel, dashed green lines indicate the limits used to verticalize the main sequence.

Figure 4.

The same as Fig. 3, but focused on the eMSTO region of NGC 1846.

We further ran spexxy on the spectra extracted from the cubes created from the observations gathered during individual nights. The main purpose of this analysis was to obtain radial velocities for the different data sets, enabling us to identify stars with variable radial velocities. For this reason, we did not try to determine σ_los in the single-night data sets, but fixed it to |$0\, {\rm km\, s^{-1}}$|⁠. Otherwise, the analyses were performed in the same way as for the full-depth data.

To reliably separate physical radial velocity variations from those caused by the finite accuracy of our measurements, a proper knowledge of the measurement uncertainties is key. We calibrated the uncertainties from the observed epoch-to-epoch variations in a similar manner to Kamann et al. (2016). It is based on the idea that (in the absence of physical variations due to binary or variable stars) the normalized velocity offsets

$$\begin{eqnarray} \delta v_{\rm los} = \frac{v_{{\rm los},\, 1}-v_{{\rm los},\, 2}}{\sqrt{\epsilon _{{\rm v},\, 1}^2+\epsilon _{{\rm v},\, 2}^2}}\, , \end{eqnarray}$$

(1)

follow a Gaussian distribution with standard deviation unity. In equation (1), |$v_{{\rm los},\, 1},\, \epsilon _{{\rm v},\, 1}$| and |$v_{{\rm los},\, 2},\, \epsilon _{{\rm v},\, 2}$| are the velocity measurements and uncertainties obtained during two different epochs for the same star. Calibrated uncertainties are obtained by multiplying the original uncertainties with the actual standard deviation of the δv_los distribution.

In contrast to Kamann et al. (2016), where the calibration was performed by grouping stars with extracted spectra of comparable S/N, we grouped the stars according to their T_eff, log g, and M/H measurements. For each star, we inferred the correction factor for the formal uncertainties of its radial velocity measurements based on the distribution of the δv_los values of the 100 closest stars in this three-dimensional space. Variable stars were iteratively removed from the comparison sample via kappa-sigma clipping. We found the calibrated uncertainties to vary depending on the brightnesses and colours of the individual stars. Velocities of stars on the red giant branch were measured to an accuracy between |$1.0\, {\rm km\, s^{-1}}$| (around m_F555W ∼ 18) and |$2.0\, {\rm km\, s^{-1}}$| (around m_F555W ∼ 20) during each epoch. Along the main sequence, the typical uncertainties (per epoch) were larger, ranging from |$5.0\, {\rm km\, s^{-1}}$| at m_F555W ∼ 20 to |$20.0\, {\rm km\, s^{-1}}$| for the faintest stars in our sample. We present the results of our variability analysis in Section 4.4 below.

3.5 Membership determination

To identify non-cluster members across the |$1\, {\rm arcmin}^2$| MUSE field of view, we adopted a maximum likelihood approach that takes into account the measured radial velocities as well as the surface brightness profile of NGC 1846. Briefly, we modelled the line-of-sight velocity distribution in the MUSE field of view as the superposition of two Gaussians, one representing the cluster and the other representing the field stars. In addition, each star was assigned a membership prior, corresponding to the surface density predicted by a King (1962) profile at the position of the star relative to the central value of the profile. The structural parameters of the profile were taken from Goudfrooij et al. (2009).

We then determined the parameters of the Gaussian profiles and a background density of stars using the affine-invariant Markov chain Monte Carlo sampler emcee (Foreman-Mackey et al. 2013). Details on this process can be found in Kamann et al. (2018a). With a set of optimized parameters at hand, we were able to assign each star a posterior membership probability of being a cluster member. We removed 190 stars with a probability <0.5 of belonging to the cluster from our sample.

When comparing the CMDs of likely members shown in Fig. 3 with the one of all extracted stars shown in Fig. 2, one notices that the stars along the ‘bridge’ connecting the bottom of the red giant branch with the central main sequence (with 21 ≲ m_F555W ≲ 23 and 1.5 ≲ (m_F435W − m_F814W) ≲ 2.0) have been discarded as non-members. Indeed, their positions in the CMD suggest that they are red giant stars from the LMC field. Note that the position of a star in the CMD was not considered in the membership selection process as this could remove potentially interesting stars in unusual locations of the CMD (such as blue stragglers or sub-subgiants) from our sample.

A detailed description of the membership determination process will be included in a forthcoming paper on the internal kinematics of NGC 1846.

4 RESULTS

4.1 Stellar rotational velocities

In the left-hand panel of Fig. 4, it is already visible that along the main sequence and in particular around the eMSTO, the stars possess a large range of rotational velocities. We measure minimum values of vsin i that are consistent with no rotation (i.e. |$\lesssim 40\, {\rm km\, s^{-1}}$|⁠, see below and Appendix  A) and maximum values of |$\sim 250\, {\rm km\, s^{-1}}$|⁠. To illustrate how rotation alters the observed spectra, we compare in Fig. 5 the spectra of two eMSTO stars that were fitted with low and high vsin i values, respectively. Both stars have magnitudes of m_F555W = 20.5, but different colours, m_F435W − m_F814W = 0.75 for the slowly rotating star |$\#31598$| and m_F435W − m_F814W = 0.90 for the fast rotating star |$\#33406$|⁠. In the spectral analysis, we found consistent metallicities and a temperature offset of |$350\, {\rm K}$|⁠, with the fast rotating star being cooler, which is in line with the correlation between photometric colour and effective temperature discussed in Section 4.2. A significant broadening of the strong |${{\rm Ca}\, \small {{II}}}$| triplet lines in the spectrum of the fast rotating star is clearly visible in Fig. 5. This shows that despite the moderate spectral resolution of MUSE (⁠|$R\sim 3\, 000$| at the wavelength of the |${{\rm Ca}\, \small {{II}}}$| triplet, corresponding to an intrinsic FWHM of |$100\, {\rm km\, s^{-1}}$|⁠), spectral differences caused by stellar rotation are detectable. However, as the simulations described in Appendix  A show, our measurements are limited to the regime |$v\sin i\gtrsim 40\, {\rm km\, s^{-1}}$|⁠. For stars with projected rotational velocities below this threshold, we can only give upper limits (i.e. confirm that they are slow rotators), whereas the exact values become sensitive to systematic errors.

Figure 5.

Comparison of the |${{\rm Ca}\,\small {{II}}}$| triplet region in the MUSE spectra of two eMSTO stars with very different measured rotational velocities, as indicated in the plot. Note that the spectra have been smoothed with a median filter of 3 pixels width for clarity.

Fig. 4 further shows that the average vsin i increases towards the red end of the eMSTO. This is in agreement with previous studies of clusters with similar ages to NGC 1846 that found a correlation between the position of a star at the eMSTO and its vsin i (e.g. Bastian et al. 2018; Kamann et al. 2018b). A comparison to the isochrone tracks depicted in Fig. 2 also suggests a good qualitative agreement between our measurements and the model predictions, with both showing an increase in the average vsin i from the blue to the red end of the eMSTO.

Thanks to the deep MUSE observations, we are able to investigate the behaviour of vsin i not only for eMSTO stars, but also for main-sequence stars below the turn-off. In order to investigate in-depth the behaviour of vsin i along the main sequence and its turn-off, we followed Bastian et al. (2018) and determined the blue and red edges of the upper main sequence. This was done by obtaining the 5th and 95th percentiles of the distribution in m_F435W − m_F814W colour as a function of m_F555W magnitude for the interval 19.8 < m_F555W < 22.0. Outliers, such as blue stragglers or field stars, were removed by restricting the colour distributions to 0.5 < m_F435W − m_F814W < 1.1.¹ The resulting edges are shown in the central panels of Figs 3 and 4 as dashed green lines. Afterwards, we assigned each star on the main sequence a pseudo-colour |$\Delta _{F435W,\, F814W}$| by measuring its normalized horizontal distance to the blue edge. We further identified as eMSTO stars all stars that passed the aforementioned selection criteria and fulfilled the condition m_F555W < 20.6. This resulted in a sample size of 1411 main-sequence stars, out of which 474 are eMSTO stars.

In Fig. 6, we analyse the variation of vsin i as a function of the |$\Delta _{\rm F435W,\, F814W}$| pseudo-colour. We split up the stars in between the two fiducial lines shown in the central panels of Figs 3 and 4 into five magnitude bins, which are depicted in order of increasing magnitude in Fig. 6. The results for the brighter two bins, which mainly contain eMSTO stars, confirm a correlation between the vsin i of a star and its (pseudo-) colour, in the sense that stars spin faster, on average, the redder they are on the eMSTO. This trend continues down the main sequence. Only for the faintest bin shown in Fig. 6, no correlation between colour and projected rotational velocity is obvious. However, this bin also seems to contain significantly less fast rotating stars compared to the other bins. We suggest that this lack of fast rotators indicates that our observations have reached a regime in stellar mass where stars are efficiently spun down by magnetic braking (see discussion in Section 5.1).

Figure 6.

The distribution of projected rotational velocities, vsin i, along the main sequence of NGC 1846. The left-hand panels show the measured vsin i values as a function of the pseudo-colour |$\Delta _{F435W,\, F814W}$| (as defined in Section 4.1) for different magnitude bins as green circles, while values considered upper limits are shown as green Y-shaped symbols. Note that in order to visualize the densities of the distributions, all symbols are semitransparent. The median uncertainty of the vsin i measurements is indicated at the top right of each panel. In the right-hand panels, we show histograms of the vsin i measurements for the same magnitude bins as in the left-hand panels. Black vertical lines indicate the locations of the peaks of Gaussian mixture models with n = 2 components.

Another observation to be made in Fig. 6 is the emergence of a separate branch of red (⁠|$\Delta _{F435W,\, F814W} \gt \sim 0.3$|⁠) but slowly rotating (⁠|$v\sin i\lt 100\, {\rm km\, s^{-1}}$|⁠) stars. This branch appears particularly prominent in the faintest three bins, i.e. at magnitudes m_F555W > 20.7. As we will discuss in Section 5.2 below, this sequence of stars is likely to be related to binary stars.

Looking at the histograms of the vsin i measurements that are shown in the right-hand panels of Fig. 6, a bimodal distribution of projected rotational velocities is suggested at least in the brighter four magnitude bins of the sample, with one maximum being located at around |$v\sin i = 60\, {\rm km\, s^{-1}}$| and the other one being located at around |$v\sin i = 140\, {\rm km\, s^{-1}}$|⁠. To verify the bimodality, we tried to represent the vsin i measurements in each bin with Gaussian mixture models containing n = 1 to n = 4 components, and used the Bayesian information criterion (BIC) to select the optimal number of components n for each bin. We found that the model with n = 2 was preferred in all five magnitude bins. The locations of the peaks suggested by this model are included in the right-hand panels of Fig. 6 as black vertical lines. In addition, we provide the means μ and dispersions σ of the Gaussian components in Table 1. The uncertainties provided in Table 1 were determined from 1 000 bootstrap samples created for each magnitude bin. We find that the slowly rotating stars peak at |$v\sin i\sim 60\, {\rm km\, s^{-1}}$|⁠, irrespective of the magnitude bin that is considered. On the other hand, some variations with magnitude are observed for the peak of the fast rotating stars. In particular, we observe a significant drop in the mean value from |$\mu _{\rm fast}=165.6\pm 3.8\, {\rm km\, s^{-1}}$| to |$\mu _{\rm fast}=122\pm 5.2\, {\rm km\, s^{-1}}$| when going from the central to the brightest magnitude bin (cf. Table 1).

Table 1.

Parameters of the Gaussian mixture models introduced in Section 4.1 to represent the vsin i histograms depicted in the right-hand panel of Fig. 6. The first two columns contain the lower and upper limits of the considered magnitude bins. The third column provides the average uncertainties of our vsin i measurements in each bin. Columns 4 and 6 contain the mean and standard deviation of the Gaussian component representing the slowly rotating stars, while the same properties for the component representing the fast rotating stars are listed in columns 5 and 7. Column 8 provides the fraction of stars in the fast rotating component.

mF555W, min	mF555W, max	|$\overline{\epsilon _{v \sin i}}$|	μslow	μfast	σlow	σfast	ffast
		(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)
19.8	20.4	12.1	58.3 ± 8.2	122.0 ± 5.2	24.4 ± 13.7	23.4 ± 11.0	0.55 ± 0.09
20.4	20.7	15.3	66.2 ± 3.3	145.3 ± 2.3	22.3 ± 9.1	21.4 ± 8.4	0.58 ± 0.04
20.7	21.0	20.9	67.0 ± 4.2	165.6 ± 3.8	29.0 ± 11.4	29.3 ± 12.3	0.51 ± 0.04
21.0	21.4	24.9	61.8 ± 9.5	150.9 ± 9.8	28.0 ± 17.6	41.4 ± 18.0	0.53 ± 0.10
21.4	22.0	28.8	66.1 ± 3.0	131.1 ± 6.1	24.5 ± 8.6	36.1 ± 16.4	0.37 ± 0.05

The Gaussian mixture models further provided us with estimates of the fractions of slow and fast rotators in each magnitude bin. As can be verified from the last column of Table 1, we find consistent fractions of about |$55{{\ \rm per\ cent}}$| of fast rotating stars. Only for the faintest bin, the fraction is significantly lower, which we attribute to the aforementioned spin-down of stars at this mass range caused by magnetic braking.

Finally, we note that at first glance, we do not measure any significant rotational velocities for stars in later evolutionary stages (cf. Fig. 3), such as sub-giants or red giant stars. This can be explained from the strong braking of the surface such stars experience while expanding after leaving the main sequence. However, as explained in Section 5.3 below, the large stellar sample allows us to constrain the time-scale on which the stars are being braked.

4.2 Effective temperatures

We further looked at the effective temperatures measured for the stars across the eMSTO of NGC 1846. We found evidence for a temperature gradient across the eMSTO of the similar-age cluster NGC 419 in Kamann et al. (2018b), as the average spectra created for the stars on both sides of the eMSTO showed a difference of |$\Delta T_{\rm eff}=270\, {\rm K}$|⁠. Our analysis of the data for NGC 1846 also reveals a well-defined correlation between effective temperature and the m_F435W − m_F814W colour, that is depicted in Fig. 7. The end-to-end variation we observe is |$\sim 600\, {\rm K}$|⁠, i.e. higher than our measurement in NGC 419. However, when we split the sample in half and determine the average effective temperatures for the blue and the red subsample (similar to our analysis for NGC 419), we obtain a difference |$\sim 300\, {\rm K}$|⁠, i.e. comparable to what we found for NGC 419.

Figure 7.

The distribution of effective temperatures measured across the extended main-sequence turn-off of NGC 1846 as a function of photometric colour. The MUSE results are shown as green diamonds whereas the coloured lines indicate the predictions from MIST isochrones with an age of log (age/yr) = 9.14 and various rotation rates.

We also show in Fig. 7 the temperature trends predicted by the MIST isochrones presented in Section 3.2, using different stellar rotation rates. While the slope of the colour–temperature relation is recovered remarkably well, the comparison to the MUSE data reveals a vertical offset of |$\sim 200\, {\rm K}$|⁠. The origin of this offset could be a mismatch of the selected isochrone parameters. We found that when we used isochrones with a higher metallicity (Z = 0.008 instead of Z = 0.006, cf. Section 3.2) and lower extinction (A_V = 0.06 instead of A_V = 0.26), the predicted effective temperatures underestimated the measured ones by |$\sim 200\, {\rm K}$|⁠. On the other hand, effective temperatures determined from low- or medium-resolution spectra can also be affected by systematic effects comparable to the observed offset (e.g. Kirby, Guhathakurta & Sneden 2008).

4.3 Metallicities

The full-spectrum fits yielded a median² metallicity of 〈[M/H]〉 = −0.55 and a scatter of |$\pm 0.15\, {\rm dex}$|⁠, which we consider as the formal uncertainty of individual metallicity measurements. This uncertainty, however, does not account for any systematic effects in the measurements. For example, the synthetic PHOENIX spectra used in the full-spectrum fits have scaled-solar abundances. Therefore, any deviations from the assumed abundance pattern, such as an enhanced or depleted [α/Fe] ratio, are likely to result in global shifts of the fitted [M/H] distribution. This could be an explanation for the difference between our median measured metallicity and the one obtained from the isochrone fits, Z = 0.006, corresponding to [M/H] = −0.37 (cf. Section 3.2). For an extensive discussion of systematic effects involved in deriving metallicities from MUSE spectra, we refer to Husser et al. (2016). We note that our spectroscopic measurement is in good agreement with the value of [Fe/H] = −0.49 ± 0.03 determined by Grocholski et al. (2006) using the |${{\rm Ca}\, \small {{II}}}$| triplet. As far as we are aware, no high-resolution measurements exist for this cluster. However, our results are within the metallicity range found by Mucciarelli et al. (2008) for a sample of other intermediate-age LMC clusters.

We further checked if any correlation exists between [M/H] and vsin i, as changes in the measured abundances with rotational velocity could help to understand how rotational mixing changes the surface structure of the stars. However, no significant correlation with vsin i was found. As mentioned above, our metallicity measurements are mainly sensitive to elements that possess strong lines in the MUSE wavelength range, such as calcium or magnesium. As it is not expected that rotational mixing has a significant impact on those elements, it is not surprising that we do not find a link between rotational velocity and surface metallicity.

4.4 Radial velocity variations

We used the radial velocities determined from the single-night spectra to see if any relation can be drawn between binarity and the measured rotational velocities. To this aim, we investigated for each main-sequence star if the night-to-night scatter in its measured radial velocity, i.e. the χ² value obtained under the assumption of a constant velocity, was consistent with the measurement uncertainties or whether it provided evidence for intrinsic variations. This was done using the method developed by Giesers et al. (2019), which assigns each star a probability p that it is radial-velocity variable. The assignment is based upon a comparison of the observed cumulative distribution of χ² values to the same distribution expected in the absence of binary stars. For each given χ² value, the latter distribution gives a prediction of the number of stars expected above this value due to random variations, whereas the former distribution yields the actual number of stars observed above this value. A comparison of the two numbers then yields the probability p that a star with the given χ² is variable.

In Fig. 8, we show the histograms of p values obtained after dividing our sample of stars into two equally sized subsamples, according to their measured vsin i relative to the median for the full sample.³ No differences are evident between the two subsamples, a two-sided Kolmogorov–Smirnov (K–S) test yields a K–S statistic of D = 0.037, corresponding to a probability of p = 0.80 that the two subsamples are drawn from the same parent distribution. When we follow Giesers et al. (2019) and consider each star with p > 0.5 as intrinsically variable, we find consistent variability fractions of |$5.4\pm 1.4{{\ \rm per\ cent}}$| and |$7.3\pm 1.5{{\ \rm per\ cent}}$| for the fast and slowly rotating subsamples, respectively.

Figure 8.

The probability distributions for radial-velocity variability of eMSTO stars, split according to whether the stars rotate more slowly (blue solid histogram) or faster (red dashed histogram) than the overall sample mean.

Note that in order to convert these numbers into a binary fraction of NGC 1846, one would need to correct for the selection effects of our study (see e.g. Sana et al. 2013; Giesers et al. 2019). The median uncertainties of our single-epoch velocity measurements are |$13.6\, {\rm km\, s^{-1}}$| and |$9.7\, {\rm km\, s^{-1}}$| for the fast and slowly rotating subsamples, respectively, and thus higher than the velocity dispersion of the cluster. On the other hand, we find a consistent binary fraction of |$6.8\pm 2.6{{\ \rm per\ cent}}$| among the RGB stars, for which our median single-epoch velocity measurement uncertainty is |$1.9\, {\rm km\, s^{-1}}$|⁠. Nevertheless, our temporal coverage is limited, consisting of four epochs spread by less than a year. Hence it is likely that we missed a fraction of the present binary stars. We further note that our observations target the cluster centre, where an overdensity of binaries is expected because of mass segregation. However, to first order both populations (slow and fast rotators) should be affected in the same way by these effects. Hence it appears unlikely that the comparable binary fractions we find are a result of the aforementioned observational biases, i.e. we see no reason to doubt that the binary fractions of the two populations are most likely similar.

Binarity is not the only possible source for radial velocity variations, as they can also be caused by stellar pulsations. In fact, at the age of NGC 1846, the upper main sequence and the eMSTO crosses the instability strip, so that the presence of variable stars – mainly δ Scuti stars – is to be expected. Salinas et al. (2016) investigated the role that such stars could play in shaping the CMDs of young to intermediate-age clusters and found that a large fraction of pulsating stars would be required for them to have a detectable effect. In their study of variability in NGC 1846, Salinas et al. (2018) predicted a total number of ∼60 δ Scuti stars in the cluster. Unfortunately, their sample of actually detected variable stars is restricted to radii larger than the half-light radius of NGC 1846, so that there is no overlap with the MUSE sample. However, it is unlikely that we misclassify a large number of δ Scuti stars as binary candidates. δ Scuti stars typically pulsate with periods |$\ll 1\, {\rm d}$| and with peak-to-peak amplitudes |$\lt 10\, {\rm km\, s^{-1}}$|⁠. Such pulsation would be hard to detect with our current data. We note, however, that strongly pulsating δ Scuti stars which are most likely to affect our binary statistics are typically found to be slow rotators (Breger 2000).

4.5 Radial distributions

Finally, we looked at the radial distributions of the samples of fast and slow rotators identified in Section 4.4. NGC 1846 is among the clusters for which Goudfrooij et al. (2011) reported differences in the concentrations of stars across the eMSTO. This result was interpreted as evidence for an age spread among the cluster stars. However, different radial concentrations could also originate from the evolution of the cluster. For example, if any of the subpopulations had a larger binary fraction, it could appear centrally concentrated because of mass segregation.

We calculated cumulative distributions for the two subpopulations, adopting the cluster centre found by Werchan & Zaritsky (2011). The distributions are shown in Fig. 9, which illustrates that we do not detect any differences in the concentrations of slow and fast rotators. Instead, when we perform a two-sided K–S test on the distributions shown in Fig. 9, we obtain a K–S statistic of D = 0.061 and a p-value of 0.20. Hence we cannot reject the null hypothesis that the fast and slow rotators share the same radial distributions. We note that the radial range covered by our data is a factor of ∼2.5× smaller compared to Goudfrooij et al. (2011) and we cannot exclude the possibility that different concentrations exist in the cluster outskirts.

Figure 9.

Cumulative fractions of slowly and fast rotating eMSTO stars in NGC 1846 as a function of distance to the cluster centre. The grey vertical line indicates the half-light radius of the cluster listed in Usher et al. (2017).

5 DISCUSSION

5.1 The distribution of rotational velocities

5.1.1 Rotational velocities across the eMSTO

Our analysis in Section 4.1 revealed a large range of projected rotational velocities amongst the eMSTO stars in NGC 1846, with values up to |$v\sin i\sim 250\, {\rm km\, s^{-1}}$|⁠. At the age and metallicity of NGC 1846, evolutionary models that account for stellar rotation, such as MIST (cf. Sect 3.2) or SYCLIST (e.g. Georgy et al. 2014), predict maximum surface velocities of |$\sim 300\, {\rm km\, s^{-1}}$| for turn-off stars that had zero-age main sequence (ZAMS) rotational velocities close to critical. The reason why our vsin i distribution does not extend to such high values is likely that vsin i measurements for stars rotating close to critical tend to underestimate the true values, because the equatorial regions, which are rotating the fastest, are gravitationally darkened (Townsend, Owocki & Howarth 2004). Additionally, part of the difference could also be caused by uncertainties in the evolutionary models. Hence, we do not interpret the difference between the maximum measured and predicted rotational velocities as evidence for a lack of stars with near-critical rotation rates in NGC 1846.

The stars we observe along the eMSTO in NGC 1846 have spectral types from early F to late A. For field stars of such spectral types in the Milky Way, Royer et al. (2007) observed distributions of vsin i extending up to |$300\, {\rm km\, s^{-1}}$|⁠, albeit only very few stars were observed with |$v\sin i\gt 250\, {\rm km\, s^{-1}}$|⁠. However, a one-to-one comparison is difficult because the field star sample covers a wide range of ages and the evolution of rotation appears to be complex for stars with stellar masses around |$2\, {\rm M_\odot }$| (Zorec & Royer 2012). The gradual decline of the surface velocities for the last |$\sim 30{{\ \rm per\ cent}}$| of the main sequence lifetimes that was observed by Zorec & Royer (2012) in this mass range may explain the lack of stars with vsin i values close to |$300\, {\rm km\, s^{-1}}$| in our sample. This scenario is also consistent with our observations where the fast rotators in the brightest part of the eMSTO are shifted towards lower velocities than the fainter populations (cf. Table 1), suggesting that they have been braked during the end of their main sequence lifetimes.

Our finding of a clear link between the position of a star at the eMSTO of NGC 1846 and its vsin i agrees with previous results in intermediate-age clusters (Bastian et al. 2018; Kamann et al. 2018b; Sun et al. 2019a). It highlights that stellar rotation is indeed the dominant mechanism for creating the observed colour spreads along the upper main sequences of such clusters, confirming the scenario initially put forward by Bastian & de Mink (2009).

We note that the correlation between vsin i and the colour of a star is non-linear. As can be verified from the left-hand panels of Fig. 6, a step increase in vsin i is observed near the blue edges of the (pseudo-)colour distributions, implying that up to |$\sim 150\, {\rm km\, s^{-1}}$|⁠, the effect of vsin i on the observed colour is small. Only for higher projected rotational velocities, a stronger change in colour arises, evident from the flattening of the distributions shown in the left-hand panels of Fig. 6 above this value. This behaviour is in agreement with previous observations and the predictions of rotating stellar models (e.g. see fig. 5 of Bastian et al. 2018).

5.1.2 The lower mass limit of fast rotating stars

Age spreads have been proposed as an alternative explanation for eMSTOs, although this scenario appears at odds with the observed correlation between cluster age and eMSTO width (Niederhofer et al. 2015). Recently, Goudfrooij et al. (2018) proposed that both stellar rotation and age spreads play a role in shaping the CMDs of intermediate-age clusters. The idea is based on the presence of a kink in the main sequences of young massive clusters, which Goudfrooij et al. (2018) attribute to the disappearance of rotational effects in stars fainter than the magnitude of the kink. For NGC 1846, Goudfrooij et al. (2018) predict rotational effects to become negligible at m_F555W ∼ 21, i.e. at a magnitude where the eMSTO is already established. However, as visible in Fig. 4, we detect fast rotating stars well below this limit, all the way down to the onset of the eMSTO at m_F555W ∼ 21.5. Hence, we conclude that there is no need to invoke age spreads to explain the eMSTO morphology of NGC 1846.

Regardless, the presence of fast rotating stars is expected to cease below a certain magnitude, as stars of lower masses are efficiently spun down by magnetic braking. Interestingly, we find a lack of fast rotating stars at the faintest magnitudes covered by the MUSE data, as can be verified from the histograms shown in Fig. 6 and the values listed in the last column of Table 1. Our faintest bin, where the lack of fast rotating stars is most obvious, includes mainly stars fainter than m_F555W = 21.5. This magnitude roughly coincides with the disappearance of the eMSTO, providing further evidence for the importance of stellar rotation in explaining its presence. At m_F555W = 21.5, the MIST isochrones predict a current stellar mass of |$1.36\, M_{\rm \odot }$|⁠.

5.1.3 Bimodal rotational distributions

Our analysis provides strong evidence for a bimodal distribution of projected rotational velocities in NGC 1846. This finding appears in agreement with the results of Mackey et al. (2008), who found photometric evidence for a bimodal eMSTO in this cluster. To our knowledge, this is the first detection of a vsin i bifurcation in an intermediate-age cluster. However, bifurcated vsin i distributions have been observed in young star clusters, e.g. amongst O- and B-stars in 30 Doradus (see Dufton et al. 2013; Ramírez-Agudelo et al. 2013). For stars of later spectral types, the populations of slow and fast rotators have been linked to the two branches of the split main sequences observed in their young host clusters (Marino et al. 2018; Sun et al. 2019b). Interestingly, these split main sequences are observed for a range in stellar mass that is comparable to the present MSTO mass of NGC 1846. This makes it likely that NGC 1846 showed a bimodal main sequence at younger ages.

The individual peaks in the vsin i distributions shown in Fig. 6 appear remarkably narrow, given that one expects any distribution of intrinsic rotational velocities to appear broadened by the different inclinations of rotation axes towards the line of sight. Hence, the narrowness of the peaks in the histograms of Fig. 6 could indicate a non-isotropic distribution of inclination angles. To investigate if our data are consistent with an isotropic distribution of inclination angles, we created mock data sets, starting by assigning intrinsic rotational velocities of |$60\, {\rm km\, s^{-1}}$| and |$160\, {\rm km\, s^{-1}}$| to |$45{{\ \rm per\ cent}}$| and |$55{{\ \rm per\ cent}}$| of the stars, respectively. These numbers were chosen to resemble the parameters obtained for the brighter bins shown in Fig. 6 (cf. Table 1). Then, we randomly chose an inclination for each star, assuming an isotropic distribution. Afterwards, we added noise to the sample, by randomly assigning each mock star the uncertainty of an observed star and shifting its vsin i accordingly. Finally, we used a two-component Gaussian mixture model and measured the widths of both components. This process was repeated 1 000 times and the median widths were compared to those listed for the observed components in Table 1.

We found that the widths of the mock distributions were narrower or comparable to those of the observed ones. Hence, our analysis is consistent with a random distribution of inclination angles. However, in that case, the intrinsic distributions of rotational velocities must be rather narrow. For example, for the second brightest magnitude bin the widths of the observed components are |$\sigma _{\rm low}=22.3\pm 9.1\, {\rm km\, s^{-1}}$| and |$\sigma _{\rm high}=21.4\pm 8.4\, {\rm km\, s^{-1}}$| (see Table 1), whereas the mock components have widths of |$\sigma _{\rm low,\, mock}=24.2\pm 1.8{\rm km\, s^{-1}}$| and |$\sigma _{\rm high}=23.0\pm 2.1{\rm km\, s^{-1}}$|⁠. This yields upper limits for the widths of the intrinsic rotational velocity distributions of the slowly and fast rotating stars of |$22.0\, {\rm km\, s^{-1}}$| and |$21.2\, {\rm km\, s^{-1}}$|⁠, respectively. We conclude that we cannot differentiate between a highly peaked rotational distribution or a non-uniform inclination angle distribution. In principle, one could also use the isochrone model predictions to try and discriminate between the two scenarios (see de Juan Ovelar et al. 2019). However, as noted in Section 3.2, the isochrones used in this work are insensitive to inclination angle variations. Therefore, we leave this aspect to a future publication focusing on data model comparisons.

5.2 The impact of binary stars

The question emerges of what causes the bimodality in the observed vsin i values. Interestingly, while Royer et al. (2007) found bimodal vsin i distributions among more massive field stars of spectral types early A or late B (see also Dufton et al. 2013), the peak at low rotational velocities disappeared for stars with spectral types comparable to those of the eMSTO stars in NGC 1846. D’Antona et al. (2015) proposed that binary stars could be responsible for the emergence of the sequence of slowly rotating stars in clusters, under the assumption that all stars were born as fast rotators. Binaries were excluded from the samples of Royer et al. (2007), and Abt & Boonyarak (2004) found lower rotational velocities compared to single stars for F-, A-, and B-stars that were in binary systems. This is interpreted as evidence for losses of angular momentum due to tidal interactions.

When going to fainter magnitudes in Fig. 6, we note the emergence of a sequence of red stars (⁠|$\Delta _{\rm F435W,\, F814W} \gt 0.3$|⁠) that do not follow the general trend of an increase in vsin i towards the red edge of the eMSTO and instead have low rotational velocities (⁠|$v\sin i\lt 100\, {\rm km\, s^{-1}}$|⁠). Their position relative to the main sequence suggests that those stars are photometric binary stars. This would also explain why the sequence disappears for the brightest bins depicted in Fig. 6, where the main sequence appears vertical so that the light contribution of a companion will not offset a star from it. Hence, our analysis suggests that the process of braking via tidal binary interactions also operates in star clusters.

But can binaries be responsible for all the slowly rotating stars in NGC 1846? Abt & Boonyarak (2004) found that the trend towards lower rotational velocities was restricted to close binaries, with orbital periods <500 d, suggesting that tidal effects are restricted to binaries in short orbits (see also de Mink et al. 2013). In a cluster like NGC 1846 with a central velocity dispersion |$\sim 5\, {\rm km\, s^{-1}}$|⁠, such system would be considered as hard binaries (i.e. binaries with a binding energy above the average kinematic energy of a cluster star) that are extremely unlikely to be destroyed in interactions with other stars. Therefore, essentially all slowly rotating stars should still be in binaries if braking via tidal interactions was the main reason for their presence. Given that we found |$\sim 45{{\ \rm per\ cent}}$| of the eMSTO stars to be slow rotators, this requires a fraction of hard binaries that seems too high to be consistent with the low binary fractions (of typically a few per cent, e.g. Milone et al. 2012) observed in Galactic globular clusters of comparable mass and density as NGC 1846. We note, however, that all estimates of the binary fraction in Galactic globular clusters are based on stars with stellar masses considerably below the stars we probe in NGC 1846. Therefore, a strong stellar mass dependence of the properties of binary systems may serve as an explanation of this discrepancy.

Our analysis of the radial velocity variations yielded binary fraction estimates of about |$6{{\ \rm per\ cent}}$| (cf. Section 4.4), comparable to estimates in Galactic globular clusters of similar masses to NGC 1846. In addition, we do not find evidence for an enhanced binary fraction among the slowly rotating stars. At first glance, this result seems at odds with the finding that photometrically detected binary stars are preferentially slow rotators. However, as discussed by Giesers et al. (2019), the MUSE observations are not sensitive to the (approximately) equal-mass binaries that are detected via photometry. On the other hand, stars with relatively faint companions will only be detected via their radial velocity variations. Therefore, our two approaches to detect binary stars cannot be readily compared to each other. Nevertheless, we would have expected an enhanced fraction of binary candidates among the slowly rotating stars in our analysis of Section 4.4 if all of them were to be in hard binary systems (even when considering that our sample of slowly rotating stars likely contains a small fraction of fast rotating stars observed at low-inclination angles). So it is likely that additional ingredients are required to explain the presence of the slowly rotating subpopulation of NGC 1846, such as a broad distribution of natal rotational velocities.

We note that Ramírez-Agudelo et al. (2013, 2015) found qualitatively similar results, i.e. a bimodal distribution of rotational velocities, for both single and binary O-stars in the young cluster 30 Doradus. However, some differences were present between the distributions. In particular, the distribution of the binaries was lacking a tail towards high vsin i values that was present in the distribution of the single stars. This difference can be explained under the assumption that these fast rotators have been spun up in binary interactions (see de Mink et al. 2013).

5.3 The spin-down of stars along the sub-giant branch

As already noted in Section 4.1, the lack of any fast rotating stars in advanced evolutionary stages is expected, because stars leaving the main sequence are efficiently being braked during their expansion on the sub-giant branch (SGB). Wu et al. (2016) previously used the morphology of the SGB to provide indirect evidence for a deceleration of fast rotating stars in the massive intermediate-age cluster NGC 419. The decrease in SGB width when moving from the blue end to the red end was interpreted as a signature of fast rotating stars being spun down. Thanks to the large MUSE sample, we are in a situation where we can study the spin-down directly, and compare the observed trend to model predictions.

To investigate how fast the stars spin-down, we show in Fig. 10 the vsin i measurements obtained for sub-giants as a function of photometric colour and compare them to the trends predicted by the MIST models for different ZAMS rotation rates. The model predictions correspond to the average surface velocity of the star at the equator, and can be considered as an upper envelope (i.e. before correcting for the effects of inclination) for the measurable rotational velocity of a star with a given |$\, \Omega /\Omega _{\rm crit,\, ZAMS}$|⁠.

Figure 10.

Projected rotational velocities of stars along the sub-giant branch of NGC 1846 as a function of photometric colour. While the MUSE measurements are shown as green diamonds, the coloured lines show the trends in average equatorial velocity predicted by the isochrones for stars with different zero-age main sequence rotation rates.

At the beginning of the SGB evolution, we still observe a wide range of rotational velocities, reaching |$v\sin i \sim 180\, {\rm km\, s^{-1}}$|⁠. About halfway along the SGB, the rotation rates are restricted to |$v\sin i \lt 100\, {\rm km\, s^{-1}}$| and at the end of the SGB, no stars with |$v\sin i \gt 50\, {\rm km\, s^{-1}}$| are observed. Hence, the spin-down appears to extend across most of the SGB evolution, which lasts about 100 − 200 Myr for stars of this mass (⁠|$1.65\, {\rm M_\odot }$| as predicted by the MIST isochrones) and metallicity.

An interesting aspect visible in Fig. 10 is that the model prediction for a ZAMS rotation rate that is 0.5 times the critical one appears to provide an upper limit for our data. This appears to be at odds with our finding for the eMSTO, were our vsin i measurements extended up to higher fractions of the predicted critical ZAMS rotational velocity. The reason for this difference can be seen in Fig. 6. In the brightest main sequence bin, the peak at high vsin i is observed at significantly lower velocities, |$\mu _{\rm fast}=122\pm 5\, {\rm km\, s^{-1}}$|⁠, compared to |$\mu _{\rm fast}\sim 150\, {\rm km\, s^{-1}}$| for the fainter bins. We interpret this behaviour as evidence that the spin-down already starts when the stars are still on the eMSTO. This also explains why in Fig. 2 the SGB appears very narrow and no stars follow the isochrone predictions for close to critical ZAMS rotational velocities.

6 CONCLUSIONS

We analysed the properties of a large sample of eMSTO and main-sequence stars in the |$\sim 1.5\, {\rm Gyr}$| old massive cluster NGC 1846 located in the Large Magellanic Cloud. Thanks to the large field of view of MUSE and our dedicated analysis tools, we were able to extract a large sample of |$1\, 411$| useful spectra of individual stars at or below the extended main sequence turn-off from the data. Compared to previous works that studied the stellar rotational velocities in young and intermediate-age clusters using 10s of stars per target, this represents an increase by a factor of ∼30 in sample size. The exquisite efficiency of MUSE further allows us to study stars as faint as m_F555W = 22, i.e. >2 magnitudes below the eMSTO of NGC 1846. To our knowledge, this is the first time that genuine individual main-sequence stars in a massive intermediate-age cluster are studied spectroscopically.

Our analysis of the projected rotation velocities confirms the results from previous studies (e.g. Bastian et al. 2018; Kamann et al. 2018b) that vsin i correlates with the position of a star on the eMSTO. However, we also find evidence for a bimodal rotation rate across the eMSTO, with two clearly separated peaks at |$v\sin i\sim 60\, {\rm km\, s^{-1}}$| and |$150\, {\rm km\, s^{-1}}$|⁠. This finding appears to be in line with previous photometric studies of NGC 1846, who reported that the eMSTO also appears to be bimodal (Mackey et al. 2008; Milone et al. 2009). It is also consistent with the observation of split main sequences in young stellar clusters that are characterized by different rotation rates (e.g. Marino et al. 2018).

We investigated if the detection of two well-defined peaks in the vsin i distribution can be considered as evidence for an anisotropic distribution of inclination angles in NGC 1846, given that even for an intrinsically bimodal distribution of rotational velocities, one would expect the bimodality to be partially washed out by the different inclinations of the stars towards the line of sight. Finding evidence for non-isotropic distributions of inclination angles could provide crucial insights into cluster formation. In particular when linked to other properties like the orientation of binary systems or the rotation of a cluster as a whole, it would allow one to study how angular momentum cascades down from the scales of giant molecular clouds to stellar scales. However, our data are still consistent with an isotropic distribution of rotation angles, provided that the intrinsic distributions are considerably peaked, with widths of |$\sigma \lesssim 20\, {\rm km\, s^{-1}}$| (assuming a Gaussian shape).

We find evidence for the disappearance of fast rotating stars at both ends of the eMSTO. At the bright end, the spin-down of stars as they evolve along the sub-giant branch can be explained by their simultaneous expansion. Our data suggest that fast rotating eMSTO stars already start their spin-down at the end of their main sequence lives, earlier than predicted by the stellar evolutionary models. At the faint end, we observe a significant drop in the number of fast rotators for the faintest stars in our sample. We attribute this observation to the magnetic braking of stars in this magnitude range. The fact that this roughly coincides with the magnitude where the spread in photometric colour disappears suggests that stellar rotation alone can explain the morphology of the eMSTO in NGC 1846.

The role that binaries play in shaping the vsin i distribution of NGC 1846 remains unclear. On the one hand, we find that photometric binary stars are preferentially slow rotators, in agreement with what is found in the Galactic field. On the other hand, the total binary fraction of the slow rotators appears to be only marginally enhanced compared to the binary fraction of the fast rotators. This suggests that while stars can be spun down by tidal interactions in tight binary systems, other factors must be at play to explain the high fraction of slowly rotating stars (⁠|$\sim 45{{\ \rm per\ cent}}$|⁠) that we observe in NGC 1846.

Finally, we also report the discovery of a new planetary nebula (PN) in NGC 1846, which is most likely a cluster member (see Appendix  B). Another PN, Mo-17 (e.g. Mackey et al. 2013), was already known in NGC 1846. Given that PNe in star clusters appear to be very rare, the coexistence of two PNe in a single cluster is quite remarkable.

ACKNOWLEDGEMENTS

The authors are grateful for a detailed report from the anonymous referee that has led to improvements in the paper. Based on observations made with the ESO Telescopes at the La Silla Paranal Observatory under Programme ID 0102.D-0268(A). SK, NB, and CU gratefully acknowledge funding from a European Research Council consolidator grant (ERC-CoG-646928-Multi-Pop). SdM acknowledges funding by the European Union’s Horizon 2020 research and innovation program from the European Research Council (ERC) (Grant agreement No. 715063), and by the Netherlands Organization for Scientific Research (NWO) as part of the Vidi research program BinWaves with project number 639.042.728. CG thanks the Equal Opportunity Office of the University of Geneva. CL thanks the Swiss National Science Foundation for supporting this research through the Ambizione grant number PZ00P2_168065. Support for this work was provided by NASA through Hubble Fellowship grant HST-HF2-51387.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This research made use of astropy,⁴ a community-developed core python package for Astronomy (Astropy Collaboration et al. 2013, 2018). This work made use of PyAstronomy and scikit-learn (Pedregosa et al. 2011).

Footnotes

1

We changed the upper limit to the 85th percentile for m_F555W < 21 given the relatively large number of stars to the red of the main sequence in this magnitude range.

2

Calculated over all likely cluster members with S/N > 10

3

We note that the results remained unchanged when we split the sample at the 45th percentile, as suggested by the fractions of fast rotators in Table 1.

4

http://www.astropy.org

5

Using the implementation in https://github.com/sczesla/PyAstronomy.

APPENDIX A: CALIBRATION OF THE v sin i SCALE

Our approach to measuring vsin i is different from most previous studies, where the widths of individual lines in high-resolution spectra have been determined. In particular, in the full-spectrum fit performed by spexxy, the line broadening is parametrized as the standard deviation σ_LOS of a Gaussian broadening kernel.

To verify that we can determine vsin i even at the limited spectral resolution of the MUSE data and to calibrate our σ_LOS measurements on a vsin i scale, we designed dedicated simulations. One obvious approach would be to take synthetic spectra, broaden them to a certain projected rotation velocity, and analyse them in the same way as the observed spectra. However, we expect that the largest uncertainty in our analysis will be systematic mismatches between the synthetic PHOENIX templates and the observed spectra, which would not be accounted for by such a simulation. Therefore, we followed a different route and took an actual observed spectrum from the X-Shooter Spectral Library (XSL, Chen et al. 2014). We selected the spectrum of HD 18769 for our test as it offers the best compromise in having stellar parameters comparable to those expected for MSTO stars in NGC 1846 and having a low projected rotation velocity itself. The latter was estimated to |$v\sin i\sim 40\, {\rm km\, s^{-1}}$| from measuring the widths of metal lines in the XSL spectrum. The stellar parameters were determined by Allende Prieto & del Burgo (2016) as |$T_{\rm eff}=8057\, {\rm K}$|⁠, log g = 3.63, and [Fe/H] = +0.15. Hence, HD 18769 is about |$+0.5\, {\rm dex}$| more metal rich than the stars of NGC 1846.

We simulated mock MUSE spectra by broadening the wavelength-dependent MUSE line-spread function (LSF) to a certain vsin i value, according to the formula of Gray (2008).⁵ The resulting kernel was used to broaden the XSL spectrum, and the result was resampled to the wavelength range and spectral sampling of the MUSE data. Finally, noise was added to the spectrum.

Our analysis needs to account for the intrinsic rotation of HD 18769. Analysis of the XSL spectrum yields |$v\sin i_{\rm HD\, 18769}\approx 40\, {\rm km\, s^{-1}}$|⁠. This value needs to be added in quadrature to the vsin i values we used to generate the mock spectra in order to obtain their true rotational broadening. Note that this step assumes that the rotational broadening kernels have Gaussian shape, which is not strictly true. In addition, this step implies that we cannot use the simulations to test the recovery of values in the range |$v\sin i \lt 40\, {\rm km\, s^{-1}}$|⁠. To address these two shortcomings, we performed a second set of simulations, using synthetic spectra, which is described below.

In total, we generated |$1\, 000$| mock spectra. Each was randomly assigned a rotation velocity from the interval |$0 \lt v\sin i \lt 200\, {\rm km\, s^{-1}}$|⁠. The MUSE LSF was recovered from the dedicated data product (LSF_PROFILE) generated by the MUSE pipeline. The pipeline recovers the LSF for every slice of every IFU, so we randomly selected the IFU and slice number for each simulated spectrum in order to verify how robust our results are against LSF mismatches. Furthermore, we applied a random radial velocity to each spectrum (from the interval |$30\, {\rm km\, s^{-1}}\lt v_{\rm LOS}\lt 70\, {\rm km\, s^{-1}}$|⁠). The S/N of the final spectra was between 10 and 100.

When analysing the mock spectra, we selected the initial guesses for the stellar parameters randomly from the ranges around the literature values, i.e. log g = 3.63 ± 0.10, |$T_{\rm eff}=(8057\pm 100)\, {\rm K}$|⁠, and [Fe/H] = +0.15 ± 0.10. The initial guesses for the radial velocity were also drawn randomly from an interval around the true value, where the width of the interval scaled inversely with the S/N of the simulated spectra. In the analysis with spexxy, we allowed every parameter to vary, only log g was fixed to the initial guess.

We show the result of our simulation in the left-hand panel of Fig. A1, where the line broadening σ_LOS determined by spexxy is plotted as a function of the vsin i used to generate the spectra. We obtain a relation that is very close to being linear. The best-fitting linear relation that passes through the origin gives

$$\begin{eqnarray} \sigma _{\rm LOS} = 0.55\times v\sin i\, . \end{eqnarray}$$

(A1)

It is reassuring that no strong dependence on the S/N of the simulated spectra is observed. Only at the lowest values, S/N < 15−20, the recovered σ_LOS are slightly too low, leading to vsin i values underestimated by |$\sim 5\!-\!10\, {\rm km\, s^{-1}}$| when using equation (A1). At higher S/N, the scatter around the relation is |$\sim 10\, {\rm km\, s^{-1}}$|⁠, confirming that we can measure vsin i to good accuracy in the simulated range from |$40-200\, {\rm km\, s^{-1}}$|⁠. Furthermore, it is reassuring that no significant impact of the LSF or the initial stellar parameters on the results is observed.

Figure A1.

Relation between the input projected rotational velocity vsin i and the recovered line broadening σ_LOS, obtained using mock MUSE spectra generated from an observed spectrum of the star HD 18769 (left-hand panel) and from a synthetic PHOENIX spectrum (right-hand panel). In both panels, different symbols sizes and colours indicate the S/N of the mock spectra. The thin black dashed lines show a one-to-one relation while the thick red dashed lines show the best-fitting relation according to equation (A1). Note that the reason for the cut-off at low vsin i visible in the left-hand panel is the intrinsic rotation of HD 18769 (⁠|$v\sin i\approx 40\, {\rm km\, s^{-1}}$|⁠).

We performed a second set of simulations in which we replaced the observed X-Shooter spectrum by a synthetic PHOENIX spectrum from the göttingen spectral library (Husser et al. 2013), with stellar parameters resembling those of eMSTO stars in NGC 1846 (⁠|$T_{\rm eff}=7200\, {\rm K}$|⁠, log g = 4.0, [M/H] = −0.5). As mentioned above, the reason for performing this second set of simulations was the intrinsic rotation of HD 18769, which for example prevented us from studying the recovery of vsin i in the regime |$\lesssim 40\, {\rm km\, s^{-1}}$|⁠.

We set up and analysed the additional simulations in the same way as described above. The results of this process are shown in the right-hand panel of Fig. A1. In comparison to the previous results, shown in the left-hand panel of Fig. A1, we notice several characteristics. First, we obtain a very similar slope for the σ_LOS–vsin i relation, 0.52 instead of 0.55 (cf. equation A1). This illustrates that the ansatz of convolving the two rotation kernels as if they were Gaussians we made in setting up the first set of simulations does not lead to a deviant vsin i scale. Second, we find that for low simulated vsin i, the relation flattens and we reach a plateau at |$\sigma _{\rm LOS}\sim 15\, {\rm km\, s^{-1}}$|⁠. We interpret this as evidence that our analysis is insensitive towards determining values of |$v\sin i \lesssim 40\, {\rm km\, s^{-1}}$|⁠.

The aforementioned flattening also implies that our measurement uncertainties become asymmetric at low vsin i. In order to quantify this effect, we determined the 16th and 84th percentile of the distribution of recovered vsin i values for different values of input (i.e. true) vsin i (using a step size of |$5\, {\rm km\, s^{-1}}$| in input vsin i). We found that the uncertainties of our measured vsin i values become increasingly asymmetric below |$100\, {\rm km\, s^{-1}}$| and are comparable to or larger than the actual measurements for |$v\sin i\lesssim 40\, {\rm km\, s^{-1}}$|⁠, again highlighting the sensitivity limit of our analysis. We incorporated this uncertainty calibration into our actual data.

APPENDIX B: SERENDIPITOUS DISCOVERY OF A PLANETARY NEBULA IN NGC 1846

As illustrated by the discovery of a nova remnant in the Galactic globular cluster M22 (NGC 6656) by Göttgens et al. (2019), MUSE is very efficient in uncovering faint nebular emission inside the 1 arcmin × 1 arcmin field of view. Here, we report the discovery of a planetary nebula in NGC 1846. As illustrated in Fig. B1, close inspection of our data revealed nebular emission at the northern edge of the MUSE pointing shown in Fig. 1. Its spectrum, obtained via simple aperture extraction of the fluxes from all spaxels inside the green circle shown in Fig. B1, shows emission lines of hydrogen, oxygen, nitrogen, and sulphur (cf. right-hand panel of Fig. B1). Based on the strong |${{\rm O}\,\small {{III}}}$| lines, the weak |${{\rm S}\,\small {{II}}}$| lines, and the ratio of the two visible |${{\rm S}\,\small {{II}}}$| lines, we classify the object as a planetary nebula (PN).

Figure B1.

Discovery of a planetary nebula in NGC 1846. The left-hand panels show an on-/off-band comparison for the |$[{{\rm O}\,\small {{III}}}]5007$| line, with the nebula visible at the northern edge of the MUSE field of view in the on-band image. The right-hand panel displays the spectrum extracted from the green circle visible in the left-hand panels. Gaussian profiles have been fitted to the labelled emission lines. The grey-shaded areas mark the sodium laser gap and spectral regions affected by strong telluric emission.

Unfortunately, the MUSE field of view only covers a fraction of the PN, so its morphology and exact location remain unknown. We also inspected the HST images available for the location under question, but none of them clearly suggested the presence of the PN, only the F555W image appeared to show a marginal glow around this position. Based on the green circle included in Fig. B1, we estimate its central coordinates to |$\alpha =05^{\rm h}\, 07^{\rm m}\, 36^{\rm s}78$| and |$\delta =-67^{\circ }\, 27\, ^\prime \, 08{^{\prime\prime}_{.}} 8$|⁠. The catalogue of PNe in the Large Magellanic Cloud by Reid & Parker (2006) does not list any objects near that location. In particular, we stress that this PN is not to be confused with Mo-17, which was discovered ∼80 arcsec south of NGC 1846 by Morgan (1994) and is likely associated with the cluster (Mackey et al. 2013). Hence we consider the object as a PN that was previously unknown.

We further fitted Gaussian line profiles to the main emission lines visible in the right-hand panel of Fig. B1. Their centroids yielded a radial velocity of |$241.7\pm 7.1\, {\rm km\, s^{-1}}$|⁠, which is in good agreement with the systematic velocity of NGC 1846. This makes it very likely that the PN is part of NGC 1846, instead of being a foreground or background object. When inserting the radial velocity and the distance to the cluster centre in our membership determination method (cf. Section 3.5), we obtain a probability of |$96{{\ \rm per\ cent}}$| that the PN is a cluster member.

We also searched the HST photometry for a possible central star of the PN. However, within a distance of 5 arcsec to the estimated centre position, no obvious candidates that are bright enough and show very blue colours were detected (cf. Fig. B2). We note the presence of one blue star (with m_F336W = 23.4 and m_F435W = 25.1), which is consistent with being a central star of a PN on the cooling track. However, given its faintness, the associated nebula should have already dissolved. It is also possible that the central star is strongly extincted.

Figure B2.

Near-UV colour–magnitude diagram of NGC 1846, created from the photometry presented in Section 3.1. Stars available in the HST catalogue within 5 arcsec of the estimated centre of the planetary nebula are highlighted according to their distance to said position.

Author notes

Hubble Fellow