ALMA detection of CO rotational line emission in red supergiant stars of the massive young star cluster RSGC1 Determination of a new mass-loss rate prescription for red supergiants

,


Introduction
The evolution of massive stars up to the point of supernovae (SNe) remains poorly understood.The steepness of the initial mass function and their short lifetimes (∼15 Myr) make such stars rare, whilst the brevity of their post main-sequence (MS) evolution makes the direct progenitors of SNe rarer still.The pre-SN mass-loss behaviour is the key property that determines the appearance of the SN, since it dictates the extent to which the envelope is stripped prior to explosion.It also determines the nature of the end-state, that is complete disruption, neutron star, black hole, or total implosion with no supernova (e.g.Heger et al. 2003).
Send offprint requests to: Leen.Decin@kuleuven.beThe most common of the core-collapse SNe are of type IIP, which are observed to have red supergiants (RSGs) as their direct progenitors.Smartt (2009) noted that the range of initial masses of these SN progenitors inferred from pre-explosion photometry, 8≤M /M ⊙ ≤17, is at odds with conventional theory, which predicts that the upper mass limit should be closer to ∼30 M ⊙ ; referred to as the 'red supergiant problem' (e.g.Ekström et al. 2012).A potential explanation for this discrepancy is that the missing RSGs (i.e.those with an initial mass between 17 -30 M ⊙ ) collapse to form black holes with no observable SNe.Later on, Davies & Beasor (2020) cautioned that this observational cutoff (of 17 M ⊙ ) is more likely to be higher and is fraught with large uncertainties (19 +4 −2 M ⊙ ) and that also the upper mass limit from theoretical models should be shifted down-wards to ∼25 -27 M ⊙ .One of the main uncertainties in both the data analysis and stellar evolution predictions is our relatively poor knowledge of RSG mass-loss rates.
Mass loss during the RSG phase can affect the progenitors of SNe in two ways.Firstly, increased mass loss can strip the star of a substantial fraction of the envelope, causing the star to evolve back to the blue before SN (Georgy 2012), and possibly depleting the stellar envelope of hydrogen (hence changing what would have been a Type-II SN into a Type-I SN).Secondly, the mass ejected can enshroud the star in dust, increasing the visual extinction by several magnitudes (e.g. de Wit et al. 2008;Beasor & Davies 2018), and causing the observer to underestimate the pre-SN luminosity of the star, or perhaps preventing the progenitor from being detectable at all (Walmswell & Eldridge 2012).Hence, accurate knowledge of RSG mass-loss rates is crucial to our understanding of stellar evolution and SN progenitors.
Despite this, the mass-loss rates ( Ṁ ) of RSGs are relatively poorly known, in comparison to the winds of hot massive stars which have been studied extensively (e.g.Sundqvist et al. 2011;Hawcroft et al. 2021;Rubio-Díez et al. 2022).The mass-loss rate recipes most often used in evolutionary models are those from de Jager et al. (1988) and Nieuwenhuijzen & de Jager (1990), which are somewhat antiquated as they are basically scaled up from red giants, and which can only predict Ṁ of field RSGs to within ±1 dex (see Fig. 2 below).These recipes assume that Ṁ scales with mass, luminosity and temperature, but do not take into account how Ṁ may change as the opacity of the circumstellar material builds up over time.What is required is a study of the mass-loss rates of samples of RSGs with uniform initial abundances and masses, where the evolutionary phase is the only variable.
There are several established methods for measuring Ṁ in cool stars.Arguably the best is to monitor the spectrum of a companion star as it passes behind the primary's wind (e.g.Kudritzki & Reimers 1978), but unfortunately there are very few such systems.Until recently, the only way to study large numbers of RSGs was to observe and model the infrared excess arising from the circumstellar dust (e.g.van Loon et al. 1999;Bonanos et al. 2010).However, most of the material is in molecular gas, and so a large (∼×200-500), uncertain correction factor must be applied to convert the measured dust mass-loss rate into a gas mass-loss rate (hereafter referred to as Ṁ SED ), while there is no information on the outflow speed or radial density profile (required to get Ṁ ).Alternatively, OH masers can be used to measure the gas wind speed.The modelling of the spectral energy distribution (SED) for an assumed gas-to-dust ratio yields a predicted wind speed that can be compared to the observed one.The scaling of the SED models to account for the difference between expected and measured expansion velocity, then yields the gas-to-dust ratios and mass-loss rates of the sample under study (see, e.g., Goldman et al. 2017).A much better way is to observe the gas using CO molecular line transitions to derive the expansion velocity and gas mass-loss rate directly from the CO line profiles -hereafter referred to as Ṁ CO (e.g., Knapp & Morris 1985;Loup et al. 1993;Josselin et al. 1998;Decin et al. 2006;Ramstedt et al. 2008;Danilovich et al. 2015).The faintness of these lines has meant that, until now, such observations were only possible for nearby bright RSGs.One exception for the detection of CO from an extragalactic RSG is IRAS 05280−6910 situated in the Large Magellanic Cloud, but the spectral resolution of the data acquired with the Herschel Space Observatory did not allow the expansion velocity to be measured (Matsuura et al. 2016).Now, with the immense gain in sensitivity provided by the Atacama Large Millimeter/submillimeter Array (ALMA), we can directly measure the expansion velocity and gas mass-loss rates of homogeneous samples of RSGs with wellconstrained distance and stellar parameters.
With this study, we aim to provide the first measurements of the gas mass-loss rates ( Ṁ CO ) of RSGs as a function of the specific RSG age.To do this, we have identified a sample of RSGs which have roughly the same masses and identical initial chemical compositions.Such samples are uniquely found in star clusters.The cluster RSGC1 at a distance, d, of ∼6 600 pc contains 14 RSGs and one post-RSG, all with initial masses ∼25 M ⊙ (Davies et al. 2008;Beasor et al. 2020).The cluster is effectively coeval, since even a large cluster age spread (0.5 Myr) would be short compared to the cluster age of 12 Myr (Davies et al. 2008).This means that the range of initial masses of the RSGs in RSGC1 must be narrow, and since all stars with this initial mass must follow virtually the same evolutionary track on the Hertzsprung-Russell (HR) diagram, the differences in the stars' luminosities are entirely due to how evolved they are.By measuring the mass-loss rates of these stars, we will be able to determine not only accurate mass-loss rates for a sample of RSGs, but also how mass-loss behaviour changes as the star evolves.In a next step, using these measurements in conjunction with model estimates of the RSG lifetimes, one can then integrate over all stars in the cluster to find the total mass lost during the RSG phase, a key property when estimating the fate of a star and the type of SN it will produce.

Observations and data reduction
The 14 RSGs (F01-F14) and 1 post-RSG (F15) in RSGC1 (see Table 2 in Davies et al. 2008) were observed with ALMA on 2015 June 9 and 11 for proposal code 2013.1.01200.S.We requested observations at both band 9, centred on the CO v=0 J=6-5 rotational line transition and band 6, centred on the CO v=0 J=2-1 rotational line transition, but only the latter were obtained.Positions and other details are given in Table A.1.The observations have three spectral windows (spw); one 'line' spw with a width of 1.875 GHz and 3840 channels to cover the CO(2-1) transition, and two 'continuum' 2 GHz spw with 128 channels each centred at 228.5 GHz and 213 GHz.The recorded line channels are not independent and the minimum effective spectral resolution of 0.977 MHz is approximately double the channel width.The spectral resolution of the continuum data is 15.625 MHz.36 antennas were used with minimum and maximum baselines of 63-783 m, providing a maximum recoverable scale of ∼4.′′ 4 in a field-of-view of 26 ′′ .The total integration per target source was 4.8 min.Standard ALMA Cycle 2 observing and quality control procedures were used 1 .The flux scale was set relative to Titan (excluding its atmospheric lines).Compact quasars J1733-1304 and J1832-1035 were used for bandpass calibration and phase-referencing, respectively.
The data were calibrated using the ALMA Quality Assurance scripts implemented in CASA (the Common Software Applications package) (McMullin et al. 2007).The estimated accuracy of the flux scale as applied to the targets is ∼7%.The targetphase-reference separation is ∼3.7 • (depending on target).Inspection of the (small) slopes in the phase-reference phase solutions, along with the probable antenna position uncertainties in 2015 (Hunter et al. 2016), suggests an absolute astrometric accuracy ≳ 1/16 of the synthesised beam, depending on the target signal to noise ratio (S/N).Davies et al. (2008) and Nakashima & Deguchi (2006), respectively (see Table A.1).The noise in the spectrum is 1.2 -1.4 mJy (see text for derivation) and is shown as (red) error bar in each panel.The dotted black line in the upper left panel is an alternative fit to the CO(2-1) line profile of F01, as discussed in App. C.
We inspected the continuum spw for each target and excluded several channels covering the SiO v=3 J=5-4 line seen around 212.582 GHz for some sources.We imaged each target, achieving a noise from the full 1.7 GHz range of σ rms ∼0.05 mJy.The synthesised beam in all images is about (0. ′′ 49×0.′′ 37) at position angle ∼-65 • , depending on the frequency.Table A.1 shows that the continuum peaks are less than half the rms in a 3 km s −1 spectral channel, so we did not perform continuum subtraction.Image cubes were made for the spw covering the CO(2-1) line for each source at 3 km s −1 velocity resolution (approximately 4 input channels), adjusted to constant velocity in the Local Standard of Rest frame (v LSR ) with respect to the CO(2-1) rest frequency of 230.538GHz.We obtained σ rms ∼1.9 mJy.We also made images at 1.3 and 10 km s −1 resolution but these do not reveal any more detections or significant details.We imaged the SiO v=3 J=5-4 line which, where detected, covered 2-4 continuum channels (width ∼22km s −1 ).In these low spectral resolution continuum spw the per-channel σ rms is ∼0.6 mJy and no other lines were detected.The parameters for all detections are given in App.A-App.C.
For five RSGs (F01, F02, F03, F04, and F13), the CO(2-1) line emission was detected, with a spatial extent < ∼ 1 ′′ (see App. C).These are the first detections of spectrally and spatially resolved CO rotational line emission of sources in an open cluster, in this particular case CO emission arising from the stellar wind of red supergiants located in RSGC1.For each of those five RSGs, the CO(2-1) line profile was extracted for a circular aperture of 0. ′′ 75 centred on the peak of the continuum emission.) and line profiles, there is contamination by interstellar medium (ISM) emission at specific frequencies; see the dashed regions in Fig. 1.In general, there is much less (or no) ISM contamination visible on the high velocity side, although there might be a noise background.Therefore, the red part of the line profile can also be used to assess the ISM contamination at lower velocities, since the CO line profiles are expected to be symmetric2 .For each source, the ISM contamination and alternative fits are discussed in App. C.
We made total intensity (zeroth moment) maps for each CO line over the uncontaminated channels.The rms-noise in the total intensity maps ranges between 21 -29 mJy/beam km/s, with corresponding peak signal-to-noise ratio between ∼12 -34 (see Sect.C).We measured the azimuthally averaged flux in annuli 200-mas thick, taking the minimum of the rms or the median average deviation as the error (see Decin et al. 2018, see Fig. C.2).These gave flux distributions with full width half maximum sizes of 460, 400, 480, 600 and 460 mas for F01, F02, F03, F04, F13, respectively, uncertainty being ∼50 mas.This gives an indication of the relative sizes of the brightest ∼60% of the emission, rather than the true size, since the flux distribution is not necessarily Gaussian and may be irregular.It is more challenging to estimate the total size of the CO emission, since our observations are sensitivity-limited and provide a lower limit.We estimated where 3× the rms noise (3×25 mJy beam −1 km s −1 ) intersected the azimuthal average profiles.This gave diameters of 700, 530, 770, 650 and 900 mas for F01, F02, F03, F04, F13, respectively but the position uncertainty is proportional to the phase noise, ∼140 mas at S/N = 3, and the weaker sources, in particular, may be more extended.

Analysis and results
The five sources detected in CO(2-1) include (i) the three M-type RSGs with the highest Ṁ SED amongst the sample of red supergiants in RSGC1 analysed by Beasor et al. (2020) (F01, F02, and F03; with Ṁ SED > ∼ 4×10 −6 M ⊙ /yr); and (ii) the peculiar RSG F13, which is anomalously red compared to the other RSGs in the cluster (Davies et al. 2008), and for which the true luminosity is difficult to determine (Beasor et al. 2020).Ten sources (F05, F06, F07, F08, F09, F10, F11, F12, F14, F15) remain undetected at a CO(2-1) rms noise value of ∼1.7 mJy/beam.For five of those sources, Beasor et al. (2020) derived a mass-loss rate between 1.8×10 −7 ≤ Ṁ SED ≤8.7×10 −7 M ⊙ /yr, implying that a line sensitivity of at least a factor 5-20 better is required to detect those sources with a similar ALMA setup.The Beasor et al. Ṁ estimates were not yet available at the moment these observations were proposed and a general value of 1×10 −6 M ⊙ /yr was used for calculating the line sensitivities in the proposal.Moreover, the CO outer envelope radius was then estimated based on Mamon et al. (1988).However, both the Ṁ -estimate (for the undetected sources) and the CO outer envelope radius (for the detected sources) turned out to be smaller (see Sect. 3.1), implying an overestimate of the actual CO(2-1) line strength.
For the five sources for which CO(2-1) rotational emission was detected, we derive the properties of the red supergiant's wind from a radiative transfer analysis.The outcomes, in particular the wind mass-loss rates, Ṁ CO , are then compared to (literature) parameters retrieved from an analysis of the dust spectral features visible in the SED, Ṁ SED .

Radiative transfer analysis of the CO(2-1) emission
To retrieve the stellar wind parameters from the CO(2-1) line profiles, we used the non-local thermodynamic equilibrium (non-LTE) radiative transfer code GASTRONOOM (Decin et al. 2006) based on a multi-level approximate Newton-Raphson (ANR) operator.The molecular line data are as specified in Appendix A of Decin et al. (2010).When ray-tracing the modelled circumstellar envelope (CSE), we used a circular model beam with the same extraction aperture (of 0. ′′ 75 diameter) to allow direct comparisons.The modelled CSEs were divided into 150 shells, evenly spaced on a logarithmic scale from the stellar radius (R ⋆ ) out to the outer envelope radius.
The effective temperatures, T eff , and stellar luminosities, L bol , were taken from Davies et al. (2008); see Table 1 3 .This translates into stellar radii of ∼1100 -1500 R ⊙ (see Table 1).Since we only have one rotational CO line, we cannot put constraints on the radial distribution of the kinetic temperature for the individual sources.We therefore assume that the gas kinetic temperature follows a power-law radial profile 3 We note that those values are slightly different than the ones used in Beasor & Davies (2016, 2018); Beasor et al. (2020).
with ϵ = 0.6, since such a power law has been shown to be a good representation of the kinetic temperature in circumstellar envelopes (Decin et al. 2006, and references therein).
For each of the five detected sources, the radial extent of the CO emission zone is larger than the beam size (of ∼0.′′ 5), but slightly lower than 1 ′′ in diameter (see Fig. C.1).Given a distance of 6 600 pc to RSGC1 and stellar radii between 7.55×10 13 -1.04×10 14cm, this translates into a CO envelope radius of 250< R out < 650 R ⋆ .While the spectra of F01 and F13 are not spatially resolved for an extraction aperture of 0. ′′ 75, the spectra of F03 and F04 show the clear absorption depth reminiscent of spatially resolving the CSE.This is in line with the angular sizes estimated from the zeroth moment maps which imply that F03 and F04 have the largest angular extents (Sect.2).These considerations yield a CO envelope radius of ∼350 R ⋆ for F01 and F13, of ∼400 R ⋆ for F02, and of ∼600 R ⋆ for F03 and F04.Those values are smaller by a factor of ∼1.5 -2 compared to the value used for the proposal preparation, that was based on the r 1/2 -value -which marks the radius where the CO abundance drops to half of its initial value -of Mamon et al. (1988) for a wind mass-loss rate of 1×10 −6 M ⊙ /yr.
To calculate the CO excitation and hence level populations, we account for excitation by the stellar photons, the microwave background and the dust radiation field (Decin et al. 2006).For the latter, we assume silicates (Decin et al. 2006) that start condensing at a temperature of ∼1750 K, which translates into a dust condensation radius, R dust , of 3-4 R ⋆ .To get the local mean radiation field at each radial point in the grid, we calculate the dust radiation field for a canonical gas-to-dust ratio (r gd ) of 200, representative for a Milky Way cluster (Beasor et al. 2020).As we discuss in App.D, the specific dust-to-gas ratio has only a minimal effect on the CO excitation by the dust radiation field for these simulations.
The parametrised β-type accelerating wind is described by The terminal wind velocity, v ∞ , is deduced from the ALMA CO(2-1) line profiles (see Table 1), uncertainty being ±3 km s −1 (see App. A).As boundary for the velocity structure, we assume that the flow velocity of the gas is equal to the local sound velocity, v s , at R 0 = R dust .For the region between R ⋆ and R dust , β is assumed to be 1/2 (Decin et al. 2006); for the region beyond R dust we follow the general conclusions from Khouri et al. (2014), Decin et al. (2020), andGottlieb et al. (2022) that the value of β should be larger to represent the slowly accelerating flow in oxygen-rich winds.We here adopt β = 3.We also include a constant turbulent velocity v turb of 3 km s −1 .
To determine the fractional abundance of CO with regard to hydrogen, we assume all photospheric carbon to be locked in CO in the circumstellar envelope (CSE).For a value of A(C) = log(C/H) + 12 as derived by Davies et al. (2009) for RSGC1, this yields a CO fractional abundance of 8.9×10 −5 .This set of stellar and circumstellar parameters allows us to retrieve the mass-loss rate, Ṁ CO , from the ALMA CO(2-1) line profiles for each of the five detected sources; the derived values are listed in Table 1, and range between 1.8-42×10 −6 M ⊙ /yr.The evolution of the four sources with mass-loss rate below ∼4×10 −6 M ⊙ /yr (F01, F02, F03, and F04) is dominated by nuclear burning, while in the case of F13 the wind mass-loss rate is currently determining its RSG evolution (see Fig. 2).For the RSGC1 red supergiants that remained undetected in the CO(2-1)    line, an upper limit on the mass-loss rate of ∼7×10 −7 M ⊙ /yr is derived.

Stellar parameters CSE parameters Star
Given this set of input parameters, the largest uncertainty in Ṁ CO arises from uncertainties in the terminal wind velocity (∼35%), and then in the distance, the outer CO envelope size, and the CO fractional abundance (each ∼20%); for additional details, readers can refer to App.D. The combined effect of the errors in the individual input parameters on Ṁ CO is a factor ∼1.4; for additional details, readers can refer to App.D. This error on Ṁ CO is lower than the errors on Ṁ SED ; see Table 1.Measurements of dust emission are usually angularly unresolved, hence rely on SEDs, and depend on more factors known to vary such as composition and size of grains, the stellar contribution to the SED, etc. for which reason estimates of the dustto-gas ratio can differ by an order of magnitude.
We note that terminal wind velocities determined solely from the half line width at zero intensity might be lower limits since it was shown by Decin et al. (2018) that line widths are sensitivity limited, and hence that in some cases higher signal-to-noise data could indicate higher v ∞ values, and hence higher Ṁ CO values.However, if the terminal velocity increases, so does the width of the full line profile and of the width between the two horns in spatially resolved, optically thin line profiles.We therefore have used all these characteristics together to determine v ∞ .
An increase in terminal velocity by 3 km s −1 , would induce an increase in Ṁ CO by ∼40% (see Table D.1).
The ALMA data also constrain the CO outer envelope radius (see Table 1).In recent studies, Groenewegen (2017) and Saberi et al. (2019) improved the calculations on CO photodissociation in circumstellar envelopes made by Mamon et al. (1988).Using the ratio of the maximum flux density over 3 times the noise as a proxy for the fractional abundance of f 0 /f CO (r) (with f 0 being the initial CO abundance), we can use Eq.(3) of Saberi et al. (2019) and their values for r 1/2 and α (their Table B.1.) to compare the theoretical predicted CO envelope size with the ALMA observations.The observed ALMA CO envelope radius of the four RSGC1 sources F01, F02, F03, and F04 is only lower by a factor ∼1.1 -1.8 compared to what Saberi et al. (2019) predicted, which is a remarkable agreement given the fact that the values for r 1/2 and α were calculated for a standard Draine (1978) interstellar radiation field (ISRF).The ISRF in the massive young cluster RSGC1 might actually be higher given the fact that newly formed stars affect the surrounding materials strongly via their UV photons.However, unlike the case of 47 Tucanae (McDonald et al. 2015), an estimate of the strength and local variation of the interstellar radiation field in RSGC1 is currently lacking.Groenewegen (2017) has shown that an increase in ISRF by a factor ∼3 leads to a decrease in photodissociation radius by a factor ∼1.5.The notable exception is F13 for which the predicted outer radius is almost an order of magnitude larger than the observed one.F13 stands also out in other aspects; for further discussion, readers are referred to Sect. 4.

Comparison to SED retrievals
In various works, Beasor et al. have derived the mass-loss rate for red supergiants in four open clusters, including RSGC1, NGC 7410, χ Per, and NGC 2100 (Beasor & Davies 2016, 2018;Beasor et al. 2020).The first 3 clusters reside in the Milky Way, while the latter one is an LMC cluster.For each of those clusters, the age and initial mass M ini of the RSGs was determined.This yielded values of 12±2 Myr and 25±2 M ⊙ for RSGC1.For the other 3 clusters, the respective values are 21±1 Myr and 10±1 M ⊙ (NGC 2100), 20±1 Myr and 11±1 M ⊙ (NGC 7419), and 21±1 Myr and 11±1 M ⊙ (χ Per).The aim of their works was to derive a new mass-loss rate prescription for red supergiants that can be used in stellar evolution codes.Beasor et al. (2020) derived a general Ṁ -luminosity relation that is dependent on the initial mass, By keeping M ini constrained, Beasor & Davies (2016, 2018) showed that the Ṁ -luminosity relation has a tighter correlation with a smaller value for the dispersion.The standard deviation for the slope is given to be 4.8±0.6; the standard deviations on the other numerical values were not listed by Beasor et al. (2020), but have been determined in App.E. The slope of Eq. ( 3) is steeper than any other mass-loss rate relations derived for red supergiants (see Fig. 2).We here re-examine the results of Beasor et al. (2020) and compare the SED mass-loss rate to the ones derived from the ALMA CO(2-1) line in the current work.
In a first step, Beasor et al. (2020) has determined a Ṁluminosity relation for all clusters in their sample by fitting the relation to their data points (see Table E.1 and App.E).We repeat the analysis applying the same IDL routine FITEXY4 and using the most conservative error estimates for both L bol and Ṁ SED (see Table E.1).The fit to Eq. ( 4) is shown in Fig. 3; the values derived for a and b are listed in Table 2, together with the Pearson correlation coefficient.Similar to Beasor et al. (2020), the standard deviation on the intercept a and the slope b is large for the cluster RSGC1.For all clusters, our values for the intercept are systematically higher than the values listed in Table 4  The values of the slope b obtained from each cluster are compatible with one another.This had led Beasor et al. (2020) to fix b to 4.8 ± 0.6 (corresponding to the weighted mean of the slope obtained for each cluster) in order to decrease the number of degrees of freedom of the fit from 8 to 5. We follow a similar strategy to limit the number of degrees of freedom, however we do not fix b from the previous exercise.We rather fit a common value of b for all four clusters simultaneously with each individual intercepts (a i ) using a multivariate linear fitting method implemented with the MPFIT Levenberg-Marquardt least-squares solver (Markwardt 2009).Uncertainties were estimated through Monte Carlo (MC) by generating 10 4 statistically equivalent data sets randomly drawn from normal distributions centred on the best fit solutions and repeating the fit for each artificial data set.The 68%-confidence interval on the best-fit parameters are given by the 0.16 and 0.84 percentiles of the distributions of obtained parameters.Fitted at face values, the errors on the individual intercepts a i are tightly correlated.To prevent this correlation, we adopt a slightly different functional form: where a 1 is the intercept for the cluster RSGC1, and ∆a i indicates the difference in intercept of each cluster i = 1,2,3,4 with regard to the intercept a 1 (hence, ∆a 1 = 0).That is, each cluster gets its own intercept (a i = a 1 + ∆a i ) and the slope is forced to be identical.The MC method yields b = 3.34 +0.48 −0.41 , a 1 = −23.83+2.15  −2.64 with corresponding ∆a i -values listed in the last column of Table 2; see right panel in Fig. 3 6 .While mathematically equivalent to directly fitting the intercept a i instead of their difference ∆a i to a reference intercept (here arbitrarily chosen to be that of RSGC1, but see App.F), the errors on a 1 , ∆a 2 , ∆a 3 , ∆a 4 are now uncorrelated, which allows for a better sense of whether the intercepts from each cluster vary from one another from a direct comparison of the ∆a i values and their errors.
To parametrise the mass-loss rate in terms of both the luminosity and the initial mass, we perform a MC fit to all four clusters together using a parametrisation similar to Eq. 3, that is  4) where the values of both a and b are free parameters in each fit.In the right panel, we force the slope b to be identical for all clusters and fit Eq. ( 5).The best-fit values are listed in the legend and in Table 2.The shaded areas provide a visualisation of the uncertainties of the fits.They were constructed through MC as the 68%-locus of all best-fit relations derived from a fit to mock data sets used to assess the errors on the best-fit parameters.Notes.The second column lists the linear Pearson correlation coefficient.The third and fourth column list, respectively, the intercept a and slope b with their standard deviation derived by fitting Eq. ( 4) to the data.The fifth column lists the fit to Eq. 5 which yields a best-fit slope b = 3.35 +0.43   −0.37   and intercept a1 = −23.96+2.00 −2.32 . (6) The combined fit yields R = −20.63+1.93 −2.38 , S = −0.16+0.03 −0.04 , and b = 3.47 +0.57−0.45 .This new mass-loss rate relation has a higher constant R and a shallower dependence on the initial mass and the luminosity as compared to Beasor et al. (2020); see Eq. 3 and Fig. 2 (full and dotted grey lines).Except for the mass-dependent intercept, these values do not agree within their respective 1sigma uncertainties; see the standard deviations for the parameters of Eq. (3) derived in App.E.
However, the uncertainties associated with the intercept R are substantial, spanning two orders of magnitude for the associated mass-loss rate.But it is crucial to acknowledge that these uncertainties pertaining to the intercept do not reflect the uncertainties on the mass-loss rates within the range of values for the initial mass and luminosities under study.By normalising the mass-loss rate, initial mass and luminosities to representative values, the resulting parameter uncertainties become more readily applicable (e.g.van Loon et al. 2005).This yields the following analytical relation with R = 1.71 +0.54  −0.44 , S = −1.63+0.30 −0.36 , and b = 3.47 +0.57−0.45 ; the uncertainties on the intercept now being of the order of 30%.

CO-based mass-loss relation for RSGs
The CO gas mass-loss rate values of those red supergiants in common with Beasor et al. (2020) are systematically lower, on average by a factor of ∼2 (see Table 1) -although we here must provide the caveat of dealing with small number statistics.CO mass-loss rates, Ṁ CO , are not afflicted with uncertain dust extinction corrections, dust-to-gas conversion ratios, and unknown expansion velocities as is the case for Ṁ SED .The main reason for the difference in Ṁ SED and Ṁ CO for the 3 RSGC1 sources analysed in both this study and by Beasor et al. (2020) is the expansion velocity which was assumed to be 25±5 km s −1 by Beasor et al. (2020), but which is lower for all RSGs in which CO(2-1) was detected.Correcting for the terminal velocities as deduced from the ALMA CO(2-1) data, Ṁ CO and Ṁ SED agree very well for the three sources in common to both studies (F01, F02, and F03) with average difference being a factor ∼0.9 and maximum percentage difference of 32% (see also last two columns in Table 1).
A point of concern on the reliability of the mass-loss rates could be binarity.The binary fraction of unevolved massive stars is thought to be above 70% (Sana et al. 2012;Moe & Di Stefano 2017;Patrick et al. 2022).For a predicted merger fraction of 20-30%, and a binary interaction fraction of 40-50%, the total RSG binary fraction is estimated around 20% (Patrick et al. 2019(Patrick et al. , 2020;;Neugent et al. 2020;Sana 2022).As discussed by Decin (2021) and Gottlieb et al. (2022), a binary system with small orbital distance can develop an equatorial density enhancement (EDE) promoting the formation of dust grains.Hence, it is expected that for those systems mass-loss rate estimates based on dust spectral features in the SED might yield too high an Ṁ -estimates, and that mass-loss rate estimates based on a COanalysis should be preferred (Decin et al. 2019).The sensitivity of the CO(2-1) integrated line flux to binary-induced morphologies has been shown to be less than a factor of ∼2, for both spiral structures (see Fig. 16 in Homan et al. 2015) and equatorial density enhancements (Decin et al. 2019), the latter morphologies being much better traced in other molecular diagnostics, such as the SiO v=0 J=8-7 transition (Kervella et al. 2016;Decin et al. 2020).
Despite of (i) the caveat of dealing with small-number statistics, and (ii) the fact that we are dealing with the 3 RSGs in RSGC1 that have the largest Ṁ SED -values from the study of Beasor et al. (2020), we can try to improve on the mass-loss rate prescription derived in Eq. ( 6).Similar to the analysis of Beasor et al. (2020), we exclude F13 from the regression analysis since it stands out in the sample of 14 RSGs in RSGC1: F13 is unusually red, has an Ṁ CO being an order of magnitude larger than the other four RSGs detected in CO(2-1), and is the only source for which the observed CO envelope size is much lower than any prediction for CO photodissociation in a circumstellar envelope (see Sect. 3.1).This tentatively suggests that another mass-loss mechanism is active (see Sect. 4.2 and Sect.4.3).The limited and biased sample led us decide to fit both the (L bol , Ṁ CO ) and (L bol , Ṁ SED ) measurements together with a four parameter set of equations: with best-fit values being R SED = 1.77 +0.58 −0.46 , S = −1.68+0.31 −0.40 , b = 3.50 +0.60  −0.46 , and ∆R = 0.49 +0.44 −0.41 ; see Fig. 4. The Ṁ CO measurements have an intercept that differs by 0.5 ± 0.4 compared to the Ṁ SED measurements; the likelihood that such a difference occurs by chance is only ∼12%.This new Ṁ -luminosity relation derived for M-type supergiants is plotted as full red line in Fig. 2, and is clearly different from the Ṁ -relation derived by Beasor et al. (2020) (full grey line in Fig. 2).
To assess the predictive power of Eq. ( 8), we have compared its predicted Ṁ -values with Ṁ CO -values derived for some wellknown Galactic M-type supergiants towards which several rotational CO lines have been observed (α Ori, µ Cep, VX Sgr, and VY CMa; see App.G).The analysis presented in App.G proves that Eq. ( 8) predicts the gas mass-loss rate for M-type red supergiants with effective temperature between ∼3200 -3800 K with good accuracy, the average difference only being ∼30% (see Table G.1).
Moreover, we can use the derived Ṁ CO values to derive the gas-to-dust ratio in the winds of individual sources.This ratio can be determined by comparing Ṁ SED -values with Ṁ COvalues (adjusted for differences in wind speed and distance, as used by different authors).Our analysis reveals a gas-to-dust ratio of 235 ±71 for the RSGC1 sources.For the Galactic sources, the gas-to-dust ratio ranges between ∼200 -550.The exception is the extreme RSG VY CMa for which a low ratio of ∼20 is derived, but there are indications that that ratio is changing throughout its mass-loss history (see App. G).
Remarkably, the new Ṁ CO values for the M-type red supergiants F01, F02, F03, and F04 are smaller than predicted by all empirical mass-loss rate prescriptions derived prior to 2020 and used in stellar evolution codes (see right panel in Fig. 2).Again here, (part of) the explanation is based on the fact that all empirical Ṁ -relations shown in Fig. 2 are based on an SED analysis that are prone to uncertain gas-to-dust ratios and expansion velocities, and that samples of stars will be flawed by a large fraction of stars that experience binary interaction.In addition, those previous studies were often biased towards samples with high mass-loss rate objects for which the infrared excess is easier to detect and model, as acknowledged by those authors; see, for example, the discussion in van Loon et al. (2005).Moreover, distances -and hence corresponding luminosities -are very uncertain for Galactic samples.This latter caveat was avoided in the studies by Beasor et al. who focused on open clusters, which eventually led in 2020 to the RSG mass-loss prescription given in Eq. (3) (Beasor et al. 2020) and shown as full grey line in Fig. 2. The only theoretical mass-loss rate prescription for RSGs is derived by Kee et al. (2021).That relation can fit the RSGC1 data under condition of an atmospheric turbulent velocity of ∼15±1 km s −1 , which is lower than the values quoted in Table 2 of Kee et al. (2021) for stars with mass around 25 M ⊙ .
In that regard, it is also worth turning our attention to F13 and to the mass-loss rate prescriptions from van Loon et al. ( 2005) and Goldman et al. (2017) which are derived for dusty RSGs, some of which displaying clear OH maser action.The relation of Goldman et al. (2017) yields a mass-loss rate that is a factor ∼10 higher than what we derive, although we need to remark that the pulsation period estimated from the period-luminosity relation given by De Beck et al. ( 2010) is very uncertain.The van Loon et al. (2005) relation predicts a mass-loss rate for F13 of 5.8×10 −5 M ⊙ /yr, only a factor 1.4 larger than our derived Ṁ CO (see Fig. 2).van Loon et al. (2005) discussed that their recipe overestimates mass-loss rates for Galactic 'optical' RSGs, on average by a factor ∼2.8 (Mauron & Josselin 2011) with standard deviation on that value being ∼2.5, hence in line with our result for F13.8).Symbols, colours and shades have the same meaning as in Fig. 3, with the exception of the addition of the Ṁ CO measurements derived in Sect.3.1 using the L bol values of Davies et al. (2008).

Phase change in mass loss
The amount of mass lost during the RSG phase, its speed, and how soon before core collapse the material is removed can have a dramatic effect on the resulting supernova light curve and spectrum (Smith et al. 2009).For luminosities above 10 5 L ⊙ , a wind mass-loss rate ≤2×10 −5 M ⊙ /yr implies that the nuclear burning rate exceeds the wind mass-loss rate (see right panel in Fig. 2), and hence that core-He burning is dominating the evolution of F01, F02, F03, and F04 (see Fig. 2).Accounting for the uncertainty in L bol (Table 1) and Ṁ CO , only F13 is above the boundary where the wind mass-loss rate dominates the star's evolution.This outcome is in line with the results from van Loon et al. (1999) and Javadi et al. (2013), who suggested that red supergiants come in two flavours, those dominated by nuclear burning which takes about 75% of the RSG lifetime, and those dominated by intense mass loss taking ∼25% of the RSG lifetimes.
The measured wind speed of the four RSGC1 sources F01, F02, F03, and F04, whose evolution is dominated by nuclear burning (see Fig. 2), is 11±3 km/s.However, in the case of source F13, the wind speed is significantly higher, at ∼22 km/s.The measured wind speeds conform to the wind speed-luminosity relation (v ∞ ∝ ZL 0.4 ) established by Goldman et al. (2017) for OH/IR stars in the Large Magellanic Cloud (LMC).However, they are noticeably lower compared to the Galactic samples (see Fig. 17 in Goldman et al. 2017).This finding aligns with the outcomes reported by Davies et al. (2009), who identified a subsolar iron content in the RSGC1 sources and indicated a metallicity between Z = 0.008 and Z = 0.02.This preliminary alignment with the relation proposed by Goldman et al. (2017) for lower metallicities tentatively corroborates the findings of Goldman et al. (2017) of expansion velocities being consistent with the predictions of dust-driven wind theory (van Loon 2000).
One advantage of our study is the ability to observe essentially the 'same' star at different stages of post-main sequence evolution within a single cluster.When considering the changes in mass-loss rate and wind speed, this leads to the hypothesis that RSGs may undergo one or multiple phase changes in mass loss during their RSG evolution.It is conjectured that there could be intense and potentially eruptive mass loss occurring for a shorter period of the RSG lifetime, which disrupts the otherwise more tranquil mass-loss process that takes place over a larger portion of the RSG branch.For further discussion on this topic, we refer to Sect.4.3.

Implications for stellar evolution
Accurate stellar mass-loss predictions are fundamental for stellar evolutionary models, in particular for the prediction of the nature of the end-products.The outcome of this study impacts the formation frequency of core-collapse supernovae of type IIP, black holes and neutron stars, and hence the frequency of gravitational wave events that can be detected with current (and to be developed) detectors.For massive stars with initial mass < ∼ 30 M ⊙ , winds during the main-sequence phase will only remove ≤0.8 M ⊙ (Beasor et al. 2021).Hence, the only evolutionary phase during which these massive stars can potentially loose a significant amount of mass is during the cool red supergiant phase, which for a star of ∼25 M ⊙ lasts ∼10 5.5 yr (Meynet et al. 2015).
The MESA evolutionary code (Paxton et al. 2019) was used to compute the evolution of a 25 M ⊙ star until core carbon depletion7 .Following Brott et al. (2011), overshooting was modelled by extending the convective region by 0.335 pressure scale heights, while winds during the main sequence and the Wolf-Rayet phase are modelled using the prescriptions of Vink et al. (2001) andHamann et al. (1995).For temperatures below 10 4 K we switch to either the prescription of Nieuwenhuijzen & de Jager (1990) or the newly derived one from Eq. ( 8).Composition and opacities are determined from solar metallicity and metal fractions as given by Asplund et al. (2009).
Both simulations evolve to very different endpoints, with the one using the Nieuwenhuijzen & de Jager (1990) wind prescription managing to strip its outer stellar envelope and evolve to become a hot Wolf-Rayet star.Using the new prescription of Eq. ( 8), only 1.93 M ⊙ of the hydrogen-rich stellar envelope is lost (out of a total envelope mass of 11.97 M ⊙ ) implying that such massive stars would explode as RSG upon core-collapse (see Fig. 5).Following Smith et al. (2009), this leads to the suggestion that if F01, F02, F03, or F04 were to explode in their current RSG phase, this would produce a Type II SN with limited level of interaction with the circumstellar material, without enough inertia to substantially decelerate the blast wave and with no substantial narrow Hα emission from the post-shock gas.Mass-loss rates of order ∼11.97/1.93= 6.20 times higher would be needed to successfully strip the H-rich stellar envelope.In that sense, F13 might be an important outlier mass-loss wise, since its Ṁ CO is more than an order of magnitude larger than for the other four RSGs.The retrieved Ṁ CO value of F13 does not fit Eq. ( 8) which is ∼10 times above the prescription.Both the strong CO(2-1) line of F13 and its high near-infrared extinction (Davies et al. 2008) indicate that F13 is surrounded by a lot of circumstellar material produced by a wind with massloss rate higher than the other four RSGC1 sources.It might be that the luminosity of F13 is underestimated given the high Ṁ CO and the fact that it is anomalously red compared to the other RSGs in RSGC1.For fitting the relation Eq. ( 8), log(L bol /L ⊙ ) should be 5.74, which would imply that it would be brighter than any known RSG and is far above the Humphreys-Davidson limit.Another possibility is an evolutionary phase change.Davies et al. (2008) has shown that F13 is spatially coincident with H 2 O, OH, and SiO maser emission8 .Such maser emission is often associated with evolved stars that are long-period variables (periods of 300 -500 days) and have winds with a very high mass-loss rate.This induces the suggestion that another mass-loss rate mode has become active in F13.Among a coeval sample of RSGs one may expect to see enhanced mass loss, and hence masers, in those objects furthest along their evolution, that have high L ⋆ /M ⋆ ratios -close to the (modified) Eddington limit -so that only small changes in the atmospheric structurefor example caused by pulsations or changes in the high opacity due to variations in the hydrogen ionisation beneath the stellar surface -makes these stars unstable to more furious episodic mass ejections.Eq. ( 8) would then represent the more quiescent mass-loss process, while F13 is an example RSG that undergoes stronger, potentially eruptive, mass loss.
In that sense, F13 shares its status as extreme RSG with VY CMa (see App. G), the only Galactic RSG for which the Ṁ CO prediction using Eq. ( 8) is a factor ∼6.6 too low.VY CMa is one of a small class of evolved massive stars characterised by extensive asymmetric ejections and multiple high massloss events lasting several hundred years (Decin et al. 2006;O'Gorman et al. 2012;Kamiński 2019;Humphreys et al. 2021) attributed to large-scale surface and magnetic activity.Both F13 and VY CMa are indicative of a stronger, potentially eruptive, mass-loss process that breaks from the prescription given in Eq. ( 8); they demonstrate that one should be careful about applying any Ṁ -prescription, including Eq. ( 8), more globally in stellar evolution models as they will not reproduce the behaviour of those extreme RSGs; see also the recent discussion by Massey et al. (2022).
If we consider F to be the fraction of time spent on RSG phase with stronger mass loss, and B the enhancement of the mass-loss rate during that stage as compared to the more quiescent mass-loss rate (here B = 10.43),we can estimate that the enhancement of a lifetime averaged mass-loss rate would be equal to (1 − F ) + B × F .For an enhancement of 6.20 and B = 10.43,we need F = 0.55, so 55% of the lifetime on the strong massloss rate RSG phase to successfully strip the H-rich stellar envelope.However, this derived fraction of F = 55% is much higher than what one derives from using number statistics to determine the fraction of time spent in the strong mass-loss rate stage.I.e.out of 14 RSGs in RSGC1, only one of them is observed in the strong mass-loss rate stage.To assess the likelihood of observing such a ratio (1 out of 14), a binomial distribution can be employed.By conducting a Bayesian analysis with a flat prior for F , the posterior distribution of F can be determined, which corresponds to a beta-distribution with shape parameters α = 2 (representing the number of stars in high mass-loss rate phase + 1) and β = 14 (representing the number of stars in the quiescent RSG mass-loss rate stage + 1).Consequently, a 90% credible interval for F is obtained as 10.9 +17 −8.5 %.This interval represents the median value along with the range between the 5th and 95th percentiles.
This tension between the outcome of number statistics and the fraction of time needed to strip the H-envelope during the RSG phase when considering both the quiescent and (potentially eruptive) high mass-loss rate phase induces the suggestion that the RSGs in RSGC1 will not be able to strip their entire H-envelope and will explode as RSG upon core-collapse.For the full H-rich envelope to be stripped, mass-loss rates much stronger than observed for F13 would need to be invoked.A potential mechanism thereof has been explored by Heger et al. (1997), who suggested that sufficiently evolved stars become dynamically unstable and exhibit large-amplitude pulsations with periods of the order of the Kelvin-Helmholtz time scale, which eventually become strong enough to dynamically eject shells of matter from the stellar surface, implying the loss of (most of) the H-rich envelope.Simulations for episodic mass ejections from common-envelope binaries yield a similar outcome of dynamical unstable mass ejections with periods in the range of a few years to a few decades, leading to a time-averaged mass-loss rate of the order of 10 −3 M ⊙ /yr (Clayton et al. 2017).The probability of observing RSGs in such a stage of large amplitude pulsation and associated strong mass loss ( > ∼ 10 −4 Ṁ /yr) is not very large; however such events have marked consequence on the appearance of the supernova explosion.
Similarly, F15 is an important outlier, as being the only post-RSG in the RSGC1 sample of Davies et al. (2008).F15 matches the picture of being a post-RSG star, with its luminosity of log(L bol /L ⊙ ) ∼ 5.36 matching that of the brightest RSGs in RSGC1.Even though no mass-loss constraint could be made, F15 could be indicative of the mass-loss process operating in the RSG shutting off as the star becomes hotter (T eff ∼ 6850 K) and transitions to a radiative envelope.If F15 is indeed a post-RSG, it indicates that the RSG mass loss managed to strip the H-rich stellar envelope to the point that it would evolve to the blue.

Conclusions
The ALMA detection of CO(2-1) emission towards RSGs residing in the open cluster RSGC1 provides us with a powerful diagnostic to derive the gas mass-loss rates of those RSGs.Of importance is that the RSG cluster stars are co-eval, which allows for stars to be studied with the same initial conditionsmass, metallicity, local environment, etc.Since the cluster stars all have roughly the same initial masses (of ∼25 M ⊙ , within a few tenths of a solar mass), the evolutionary path should be the same, allowing for the luminosity to be used as a proxy for evolution.Based on the CO(2-1) detections, we propose a new massloss rate relation for M-type RSGs with effective temperatures between ∼3200 -3800 K, that scales with luminosity and mass.The new Ṁ -luminosity relation proposed in Eq. ( 8) is validated against some other well-known Galactic RSGs towards which multiple CO rotational lines have been observed.
The gas mass-loss rates derived from CO diagnostics are systematically lower than the values retrieved from an SED analysis on which current stellar evolution codes are based ( Ṁ CO < Ṁ SED ).Implementing our new mass-loss rate relation will impact the frequency rate of type IIP SNe, neutron stars, and black holes.In particular, models suggest that the RSG mass loss would not allow single massive stars to evolve back to the blue and explode as a H-poor SN.However, the mass-loss rate of both the RSG F13 in RSGC1 and the well-known Galactic extreme RSG VY CMa are almost an order of magnitude higher than predicted by Eq. ( 8), which is indicative of a stronger massloss process different than captured by Eq. ( 8).Statistical reasoning implies that the RSGs in RSGC1 will not be able to strip their entire H-rich envelope and will explode as a RSG upon core collapse.
Only five RSGs in RSGC1 were detected during the current observation run.A completion of this RSGC1 study with ALMA should allow a more accurate mass-loss rate relation to be derived, that can be checked against lower mass-loss rate RSG stars with lower luminosities.Given the fact that ALMA is now 50 -100% more sensitive than in 2015, such deeper observations are now well feasible and would allow the observation of large samples that are, as much as possible, uniform and unbiased.In addition, we intend to return to observe other clusters with slightly different ages, such as RSGC2 (Davies et al. 2007), which will then allow us to repeat this study but for RSGs with slightly different masses.Ultimately, we will be able to provide the time-averaged mass-loss rates and total mass lost during the RSG phase as a function of initial mass, crucial inputs for the theory of stellar evolution, and SN progenitors.

Appendix A: ALMA observation
For each star, the input for the ALMA observations is listed in Table A.1.For the five sources detected in the ALMA band 6 continuum data and in CO(2-1) line emission, Table A.1 lists the coordinates of the peak of the continuum emission, the peak intensity, and the rms noise of the continuum and CO(2-1) observations.For each source, Davies et al. (2008) has estimated the local standard of rest velocity, v LSR , based on high spectral-resolution observations of the 2.293 µm CO band head (see Table A.1), uncertainty being ±4 km s −1 .These v LSR values were used as input for the ALMA observations.For sources F01, F02 , F04, and F13, Nakashima & Deguchi (2006) also measured v LSR using SiO maser data, with an uncertainty of ±2 km s −1 (see last column in Table A.1).Using the ALMA data, we have re-assessed the v LSR values based on the considerations that (i) the CO line profile is symmetric around the v LSR value (see also footnote 2), and (ii) for spatially resolved sources for which the CO optical thickness is not too high (typically, a maximum τ ν of ∼2), the CO profile is two-horn like (see Fig. 1).Those values are listed in Table 1 and in Table A.1.Given the spectral resolution of the data, the uncertainty on the v LSR values is ±3 km s −1 .The newly derived v LSR values for all sources agree with the values from Nakashima & Deguchi (2006).For the sources F03, F04, and F13, the v LSR values also agree with Davies et al. (2008), F02 is just within the uncertainty range of both studies, but the value is off for F01.In particular, in the case of F01 Davies et al. (2008) lists a value of 129.5 km s −1 , but the ALMA CO and SiO maser data of Nakashima & Deguchi (2006) indicate a value around 117.5 km s −1 .Notes.First part of the table lists the target identifier, the input coordinates for the ALMA observations, and the vLSR as derived by Davies et al. (2008).The second part of the table lists for the four sources detected in the ALMA band 6 continuum data and in CO(2-1) the coordinates of the peak of the continuum emission, the peak flux, and the rms noise of the continuum observations.The third part lists the CO(2-1) line rms noise and the vLSR deduced from the ALMA CO(2-1) data.The last part list the vLSR values deduced by Nakashima & Deguchi (2006) from the analysis of SiO masers. (a) The stellar IDs from Figer et al. (2006). (b) From Davies et al. (2008). (c) From Nakashima & Deguchi (2006).
Article number, page 13 of 29    3 of Beasor & Davies (2016), Table 2 of Beasor & Davies (2018), and Table 2 of Beasor et al. (2020).Notes.The second and third column list, respectively, the intercept a and slope b with their standard deviation as derived by Beasor et al. (2020).
The fourth column lists the intercept a and its standard deviation for a slope fixed to b = 4.8±0.6.b) 11 (e) 1570 (b) 6.5 22 3400 (l) 25 28 13 VY CMa 5.48 (g) 2800 (g) M2/4II (b) 1680 (f ) 25 (h) 1200 (h) 22.6 (g) 35 (g) 950 (g) 80 (g) 5 1.4 Notes.Listed are the stellar luminosity L⋆, the effective temperature T eff , the spectral type, the stellar radius R⋆, the initial mass Mini, the distance D, the local standard of rest velocity vLSR, the terminal wind velocity v∞, the radius of the CO envelope Rout, the gas mass-loss rate Ṁ CO, the mass-loss rate as predicted from Eq. ( 8), and the mass-loss rate predicted using the luminosity-Ṁ -relation of Beasor et al. (2020) (Eq.( 3)).  i) Estimates of µ Cep's distance vary between 390±140 and 1818±661 pc (Montargès et al. 2019, and references therein).We here take the average value of the distance estimate of Montargès et al. (2021) (based on physical considerations on the relative size of the MOLsphere, D = 641 +148 −144 pc) and of Davies & Beasor (2020) (based on the average parallax of neighbouring OB stars, under the assumption that the RSG is part of the same association; D = 940 +140 −40 pc).I.e., D = 790 pc.For a change of distance from 790 pc to 390 pc, the retrieved Ṁ CO changes from 3×10 −6 to 1×10 −6 M⊙/yr. (j) The wind of Betelgeuse has various components.O' Gorman et al. (2012) proposes that the S2 flow of α Ori extends out to a radius of 17 ′′ (or 800 R⋆), although the measured intensity distribution of CO emission as a function of projected radius extends to ∼8.5 ′′ (∼400 R⋆).Changing the radius of the CO envelope from 800 R⋆ to 400 R⋆ changes the retrieved Ṁ CO from 0.4×10 −6 M⊙/yr to 0.5×10 −6 M⊙/yr. (k)  Montargès et al. (2019) have used the NOEMA interferometer to get a channel map of the CO(2-1) emission of µ Cep at a spatial resolution of 0. ′′ 92×0.′′ 72, with maximum recoverable scale (MRS) being 8 ′′ .Emission is detected up to ∼3. ′′ 5 from the central star (or ∼600 R⋆). (l) For a mass-loss rate around 2×10 −5 M⊙/yr, the predicted CO photodissociation radius is ∼3500 R⋆ (Groenewegen & Saberi 2021) or a full extent of 17 ′′ .The only CO interferometric data currently available for VX Sgr have been obtained in the framework of the ALMA ATOMIUM large program.However, the MRS of those data is only ∼8-10 ′′ and hence cannot be used to estimate the CO envelope size.For that reason, we resort to the predictions of Groenewegen & Saberi (2021).

Fig. 1 .
Fig. 1.CO(2-1) line profiles of 5 red supergiants in RSGC1.The ALMA data are plotted as red histograms.Synthetic line profiles (see Sect. 3) are overplotted as black solid lines.The grey dashed regions indicate frequency regions that are contaminated by the ISM, a strong noise background, or other genuine emission (see App. A).The vertical dashed blue lines indicate the local standard of rest velocity, vLSR, as deduced from the ALMA CO data (see App. A).The dotted and dashed-dotted blue lines indicate the vLSR value as determined byDavies et al. (2008) andNakashima & Deguchi (2006), respectively (see TableA.1).The noise in the spectrum is 1.2 -1.4 mJy (see text for derivation) and is shown as (red) error bar in each panel.The dotted black line in the upper left panel is an alternative fit to the CO(2-1) line profile of F01, as discussed in App. C.
Fig. 1 shows the observed line profiles.The noise in the spectrum (clear of ISM contamination) is given by the σ rms values from Table A.1 divided by the square root of the number of beams, resulting in a spectral noise of ∼1.2 -1.4 mJy.As visible both in the CO(2-1) channel maps (Fig. C.3-Fig.C.7

Notes.
Listed are the stellar luminosity L bol , the effective temperature T eff , the spectral type, the stellar radius R⋆, the local standard of rest velocity vLSR (see discussion in App.A), the dust condensation radius R dust , the terminal wind velocity v∞, and the mass-loss rate as deduced from the CO(2-1) line, Ṁ CO, and as deduced from an SED analysis byBeasor et al. (2020), Ṁ SED.The last two columns compare a density measure, Ṁ /v∞, based on the values deduced in this study (column 11) and used in Beasor et al. (2020) who assumed a terminal wind velocity of 25 km s −1 (column 12). (a) From Davies et al. (2008). (b) As deduced from the ALMA CO(2-1) lines.The uncertainties on Ṁ CO are discussed in App.C -D. (c) From Beasor et al. (2020).(d)F13 was classified as K2 byDavies et al. (2008).However, later on it became clear that the CO-spectral type correlation is flawed if a star has a strong stellar wind, and that under these circumstances also the spectral-type -T eff relation fromLevesque et al. (2006) does not hold.We therefore use the effective temperature and spectral type classification fromMessineo et al. (2021). (e) Beasor et al. (2020) have updated the luminosities of the RSGC1 sources published by Davies et al. (2008).In particular, Beasor et al. (2020) derived log L bol (F01) = 5.58 and log L bol (F03) = 5.33.Changing the luminosities to the values from Beasor et al. (2020) induces an increase in Ṁ CO of 10% for F01 and of 3% for F03.

Fig. 2 .
Fig. 2. Mass-loss rate as a function of luminosity.Various mass-loss rate relations derived for red supergiants are shown for a fixed stellar mass of 10 M⊙ (left panel) and of 25 M⊙ (right panel) for an assumed effective temperature of 3450 K(Reimers 1975;de Jager et al. 1988;Nieuwenhuijzen & de Jager 1990;Salasnich et al. 1999;van Loon et al. 2005;Beasor et al. 2020;Goldman et al. 2017;Kee et al. 2021).Empirical mass-loss rate relations are displayed with a solid line, the theoretical relation ofKee et al. (2021) is shown with a dashed line for 2 different values of the atmospheric turbulent velocity v atm turb .For the empirical relation of Goldman et al. (2017), we use the RSG period-luminosity relation as given in Eq. (C.17) of De Beck et al. (2010), which is valid for pulsation periods between ∼300 -800 days.The corrected Ṁ SED-relation based on the data of Beasor et al. (see Eq. (6)) is shown as dotted grey line; the new Ṁ CO-relation derived from the ALMA CO(2-1) measurements of the M-type RSG winds in RSGC1 (Eq.(8)) is shown as full red line.The rate at which hydrogen and helium are consumed by nuclear burning are shown as thick dashed-triple dotted lines; the single-scattering radiation pressure limit for an expansion velocity of 12 km s −1 is shown as dashed dark grey line.Stellar mass loss rules the evolution of RSG stars if the wind mass-loss rate exceeds the nuclear burning rate, as indicated by the light-blue region; the nuclear-burning dominated region is indicated by the light-orange region.The red triangles in the right panel indicate the (L, Ṁ )-values as derived in Sect.3.1 from the ALMA CO(2-1) line profiles of F01, F02, F03, F04, and F13.

Fig. 3 .
Fig. 3. Ṁ SED -luminosity relations for the four open clusters RSGC1, NGC 2100, NGC 7419, and χ Per.The coloured open symbols represent the (L bol , Ṁ SED)-values for the four clusters studied by Beasor & Davies (2016) and Beasor et al. (2020); error bars indicate their most conservative error estimates (see App. E).The dashed lines in the left panel show the individual fits to Eq. (4) where the values of both a and b are free parameters in each fit.In the right panel, we force the slope b to be identical for all clusters and fit Eq. (5).The best-fit values are listed in the legend and in Table2.The shaded areas provide a visualisation of the uncertainties of the fits.They were constructed through MC as the 68%-locus of all best-fit relations derived from a fit to mock data sets used to assess the errors on the best-fit parameters.

Fig. 4 .
Fig. 4. Ṁ -luminosity relations parametrised in terms of log L bol and Mini and fitted jointly for the four open clusters RSGC1, NGC 2100, NGC 7410, and χ Per.Left panel presents the best-fit to the Ṁ SED measurements of the four clusters (Eq.7) while the right panel includes in the fit the Ṁ CO-values of four stars in RSGC1 (purple downward triangles) using Eq.(8).Symbols, colours and shades have the same meaning as in Fig.3, with the exception of the addition of the Ṁ CO measurements derived in Sect.3.1 using the L bol values ofDavies et al. (2008).

Fig. 5 .
Fig.5.Evolutionary track as computed with the MESA evolutionary code for a star of initial mass 25 M⊙ applying the RSG mass-loss rate prescription from Nieuwenhuijzen & de Jager (1990) (blue) and from Eq. (8) (orange).To discern long-lived versus short thermal evolution phases, dots have been added every 10 5 years; tracks go till carbon depletion indicated by the open square.The inset at the RSG, for a range in log T eff of 0.01 dex, shows how both tracks digress there.
Values used for Fig. 3 and Fig. E.1, and reproduced from Table

Fig
Fig. E.1.Ṁ -luminosity relations as determined by Beasor et al. (2020).The coloured symbols represent the (L bol , Ṁ SED)-values for four open clusters as derived by Beasor & Davies (2016); Beasor et al. (2020); error bars indicate their most conservative error estimates.The straight lines in both panels show the individual straight-line fits to each relation log( Ṁ SED/M⊙yr −1 ) = a + b log(L bol /L⊙) for all clusters in the sample.In the left panel, the dashed lines show the Ṁ -L bol relation using the values for the offset and slope as given in Table 4 of Beasor et al. (2020); the full lines in the right panel show the fits to the Ṁ -L bol relation once the gradient is fixed to b = 4.8, following Eq. 4 of Beasor et al. (2020) and using the values for the initial mass as given in Table1ofBeasor et al. (2020).The red dotted line in the left panel shows the fit to RSGC1 for an intercept a = -53; see footnote 7. The red dashed-triple dotted lin in the right panel shows the fit to RSGC1 following Eq.(E.1).

Fig
Fig. F.1.Ṁ -luminosity relations for the four open clusters RSGC1, NGC 2100, NGC 7410, and χ Per.Dashed gold, red, black, and blue lines present the best-fit of Eq. (5) to the (L bol , Ṁ SED) measurements of the four clusters (shown as open symbols with corresponding colour), and the dashed purple line the best-fit to the (L bol , Ṁ CO) values of four stars in RSGC1 (filled purple downward triangles).

Table 1 .
Stellar and CSE parameters for the five RSGs in RSGC1 for which CO(2-1) emission was detected.

Table 2 .
Best-fit parameters for the Ṁ SED-luminosity relation for each cluster.

Table A .
1. Data for the stars observed within the ALMA proposal 2013.1.01200.S.

Table E .
2. Best-fit parameters for the Ṁ SED-luminosity relation for each cluster derived by fitting Eq. (4) to the data.
Article number, page 25 of 29 Table G.1.Stellar and CSE parameters for the red supergiants observed in various CO rotational lines by De Beck et al. (2010).bol /L ⊙ ) T eff Spectral R ⋆