© The Authors 2024

1. 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).

The 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 and 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${}_{-2}^{+4}{M}_{\odot }$) and that also the upper mass limit from theoretical models should be shifted downwards 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 well-constrained 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 ∼6600 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.

2. 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$\stackrel{″}{.}$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¹. 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 target-phase-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).

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$\stackrel{″}{.}$49 × 0$\stackrel{″}{.}$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⁻¹ 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⁻¹ 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.538 GHz. We obtained σ_rms ∼ 1.9 mJy. We also made images at 1.3 and 10 km s⁻¹ 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 ∼22 km s⁻¹). 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 Appendices A–C.

For five RSGs (F01, F02, F03, F04, and F13), the CO(2–1) line emission was detected, with a spatial extent ≲1″ (see Appendix 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$\stackrel{″}{.}$75 centred on the peak of the continuum emission. Figure 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 (Figs. C.3–C.7) 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 symmetric². For each source, the ISM contamination and alternative fits are discussed in Appendix C.

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 Appendix A). The vertical dashed blue lines indicate the local standard of rest velocity, vLSR, as deduced from the ALMA CO data (see Appendix A). The dotted and dashed-dotted blue lines indicate the vLSR value as determined by 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 Appendix 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 and 29 mJy beam⁻¹ km s⁻¹, with corresponding peak signal-to-noise ratio between ∼12 and 34 (see Appendix 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⁻¹ km s⁻¹) 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.

3. 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⁻⁶ 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⁻⁷ ≤ Ṁ_SED ≤ 8.7 × 10⁻⁷ 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⁻⁶ 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.

3.1. 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$\stackrel{″}{.}$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³. 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

$$\begin{array}{c}\hfill T_{\mathrm{gas}}(r)=T_{\mathrm{eff}}{\left(\frac{R_{\star }}{r}\right)}^{ϵ}\end{array}$$(1)

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$\stackrel{″}{.}$5), but slightly lower than 1″ in diameter (see Fig. C.1). Given a distance of 6600 pc to RSGC1 and stellar radii between 7.55 × 10¹³ and 1.04 × 10¹⁴ cm, 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$\stackrel{″}{.}$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⁻⁶ 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 Appendix 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

$$\begin{array}{c}\hfill v_{\mathrm{gas}}(r)\simeq v_{\infty }{(1-\frac{R_0}{r})}^{\beta }.\end{array}$$(2)

The terminal wind velocity, v_∞, is deduced from the ALMA CO(2–1) line profiles (see Table 1), uncertainty being ±3 km s⁻¹ (see Appendix 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₀ = 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), and Gottlieb 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⁻¹.

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⁻⁵.

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 and 42 × 10⁻⁶ M_⊙ yr⁻¹. The evolution of the four sources with mass-loss rate below ∼4 × 10⁻⁶ 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⁻⁷ M_⊙ yr⁻¹ is derived.

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 of Kee et al. (2021) is shown with a dashed line for 2 different values of the atmospheric turbulent velocity $v_{\mathrm{turb}}^{\mathrm{atm}}$. 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 and 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.

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 Appendix 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 Appendix 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 dust-to-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⁻¹, 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₀/f_CO(r) (with f₀ 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.

3.2. 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,

$$\begin{array}{c}\hfill log({\dot{M}}_{\mathrm{SED}}/M_{\odot }{\mathrm{yr}}^{-1})=(-26.4-0.23\times M_{\mathrm{ini}}/M_{\odot })+4.8log(L_{\mathrm{bol}}/L_{\odot }).\end{array}$$(3)

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 Appendix 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

$$\begin{array}{c}\hfill log({\dot{M}}_{\mathrm{SED}}/M_{\odot }{\mathrm{yr}}^{-1})=a+blog(L_{\mathrm{bol}}/L_{\odot })\end{array}$$(4)

to their data points (see Table E.1 and Appendix E). We repeat the analysis applying the same IDL routine FITEXY⁴ 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 of Beasor et al. (2020; copied in Table E.2), and lower for the slopes, although they agree within the 1-sigma uncertainties, the exception being NGC 2100. In particular, for RSGC1 we derive a = −32.98 ± 18.43 and b = 5.03 ± 3.44, as compared to Beasor et al. (2020) who listed a = −52.0 ± 51.2 and b = 8.8 ± 9.5 (see left panel in Fig. E.1 and discussion in Appendix E)⁵.

Fig. 3. ṀSED–luminosity relations for the four open clusters RSGC1, NGC 2100, NGC 7419, and χ Per. The coloured open symbols represent the (Lbol, Ṁ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 Appendix 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 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.

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⁴ 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:

$$\begin{array}{c}\hfill log({\dot{M}}_{\mathrm{SED}}/M_{\odot }{\mathrm{yr}}^{-1})=a_1+\mathrm{\Delta }{a}_i+blog(L_{\mathrm{bol}}/L_{\odot }),\end{array}$$(5)

where a₁ 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₁ (hence, Δa₁ = 0). That is, each cluster gets its own intercept (a_i = a₁ + Δa_i) and the slope is forced to be identical. The MC method yields $b=3.{34}_{-0.41}^{+0.48}$, $a_1=-23.{83}_{-2.64}^{+2.15}$ with corresponding Δa_i-values listed in the last column of Table 2; see right panel in Fig. 3⁶. 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 Appendix F), the errors on a₁, Δa₂, Δa₃, Δa₄ 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

$$\begin{array}{c}\hfill log({\dot{M}}_{\mathrm{SED}}/M_{\odot }{\mathrm{yr}}^{-1})=R+S{M}_{\mathrm{ini}}/M_{\odot }+blog(L_{\mathrm{bol}}/L_{\odot }).\end{array}$$(6)

The combined fit yields $R=-20.{63}_{-2.38}^{+1.93}$, $S=-0.{16}_{-0.04}^{+0.03}$, and $b=3.{47}_{-0.45}^{+0.57}$. 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 1-sigma uncertainties; see the standard deviations for the parameters of Eq. (3) derived in Appendix 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

$$\begin{array}{c}\hfill log({\dot{M}}_{\mathrm{SED}}/{10}^{-5}{M}_{\odot }{\mathrm{yr}}^{-1})=R+S{M}_{\mathrm{ini}}/10M_{\odot }+blog(L_{\mathrm{bol}}/{10}^5{L}_{\odot }),\end{array}$$(7)

with $R=1.{71}_{-0.44}^{+0.54}$, $S=-1.{63}_{-0.36}^{+0.30}$, and $b=3.{47}_{-0.45}^{+0.57}$; the uncertainties on the intercept now being of the order of 30%.

4. Discussion

4.1. 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⁻¹ 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, 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 CO-analysis 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 Sects. 4.2 and 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:

$$\begin{array}{cc}\hfill log({\dot{M}}_{\mathrm{SED}}/{10}^{-5}{M}_{\odot }{\mathrm{yr}}^{-1})& =R_{\mathrm{SED}}+S{M}_{\mathrm{ini}}/10M_{\odot }\hfill \\ \hfill & +blog(L_{\mathrm{bol}}/{10}^5{L}_{\odot })\hfill \\ \hfill log({\dot{M}}_{\mathrm{CO}}/{10}^{-5}{M}_{\odot }{\mathrm{yr}}^{-1})& =(R_{\mathrm{SED}}+\mathrm{\Delta }R)+S{M}_{\mathrm{ini}}/10M_{\odot }\hfill \\ \hfill & +blog(L_{\mathrm{bol}}/{10}^5{L}_{\odot }),\hfill \end{array}$$(8)

with best-fit values being $R_{\mathrm{SED}}=1.{77}_{-0.46}^{+0.58}$, $S=-1.{68}_{-0.40}^{+0.31}$, $b=3.{50}_{-0.46}^{+0.60}$, and $\mathrm{\Delta }R=0.{49}_{-0.41}^{+0.44}$; 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).

Fig. 4. Ṁ–luminosity relations parametrised in terms of log Lbol 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 Lbol values of Davies et al. (2008).

To assess the predictive power of Eq. (8), we have compared its predicted Ṁ-values with Ṁ_CO-values derived for some well-known Galactic M-type supergiants towards which several rotational CO lines have been observed (α Ori, μ Cep, VX Sgr, and VY CMa; see Appendix G). The analysis presented in Appendix 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 Ṁ_CO-values (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 and 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 Appendix 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. (2020) 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⁻¹, 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⁻⁵ 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.

4.2. 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⁵ L_⊙, a wind mass-loss rate ≤2 × 10⁻⁵ 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.

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 depletion⁷. 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) and Hamann et al. (1995). For temperatures below 10⁴ 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.

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 105 years; tracks go till carbon depletion indicated by the open square. The inset at the RSG, for a range in log Teff of 0.01 dex, shows how both tracks digress there.

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 mass-loss 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₂O, OH, and SiO maser emission⁸. 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 structure – for 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. Equation (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 Appendix 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 mass-loss 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. (2023).

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 mass-loss 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_{-8.5}^{+17}$%. 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⁻³ 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⁻⁴ M 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.

5. 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 conditions – mass, 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 mass-loss 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 mass-loss 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.

Acknowledgments

The authors acknowledge the referees for their constructive feedback. The authors thank Emma Beasor and Benjamin Davies for fruitful discussions. L.D. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreements No. 646758: AEROSOL with PI L. Decin) and from the FWO grant G099720N. L.D., H.S., and P.M. acknowledge support from the KU Leuven C1 excellence grant MAESTRO C16/17/007 (PI L. Decin). P.M. acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N. This paper makes use of the following ALMA data: ADS/JAO.ALMA2013.1.01200.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The authors acknowledge the UK Science and Technology Facilities Council (SFTC) IRIS for provision of high-performance computing facilities. This work was partly performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility. The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1.

References

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⁻¹. 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⁻¹ (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⁻¹. 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⁻¹, but the ALMA CO and SiO maser data of Nakashima & Deguchi (2006) indicate a value around 117.5 km s⁻¹.

Appendix B: ALMA SiO and continuum images

For the five sources in RSGC1 detected in the continuum and in CO v = 0 J = 2-1 emission, bright SiO v = 3 J = 5-4 line emission (at 212.582 GHz) is present as well. No other sources showed either line or had significant continuum detections.

Given the weak continuum flux for the five detected sources, we have used the bright SiO emission for locating the source and its continuum emission. Fig. B.1 shows the SiO v = 3 J = 5-4 zeroth moment (i.e., total intensity) maps, and Fig. B.2 the continuum emission. The location of the maximum of the continuum emission and the maximum continuum flux are listed in Table A.1. The continuum rms noise is listed in Table A.1, and is ∼0.05 mJy/beam. The rms noise in the SiO zeroth moment maps is 28, 34, 29, 34 and 35 mJy/beam km/s for F01, F02, F03, F04, and F13, respectively.

Fig. B.1. ALMA band 6 SiO v = 3 J = 5-4 zeroth moment maps of five red supergiants in RSGC1: F01 (upper left), F02 (upper right), F03 (middle left), F04 (middle right), and F13 (bottom left). The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, with respect to the input coordinates for the ALMA observations (see columns 2 and 3 in Table A.1). The ALMA synthesised beam is shown as a white ellipse in the lower right corner of each panel. The contours (in orange) show the SiO emission at [3, 5, 10] times the rms value.

Fig. B.2. ALMA band 6 continuum maps of five red supergiants in RSGC1: F01 (upper left), F02 (upper right), F03 (middle left), F04 (middle right), and F13 (bottom left). Contours (in orange) are indicated at [50, 60, 90, and 99]% of the peak continuum emission (see Table A.1). The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, with respect to the peak of the SiO v = 3 J = 5-4 zeroth moment map. The ALMA synthesised beam is shown as a white ellipse in the lower right corner of each panel.

Appendix C: ALMA CO v = 0 J = 2-1 total intensity and channel maps

In this section, we first present the ALMA CO v = 0 J = 2-1 total intensity (zeroth moment) and channel maps for all sources that were detected in CO(2-1); i.e. F01, F02, F03, F04, and F13. Thereafter, we discuss for each source individually the spectral regions contaminated by emission from the ISM.

For all sources, the continuum was less than the spectral channel rms so no continuum subtraction was performed. Fig. C.1 shows the CO v = 0 J = 2-1 total intensity map of F01, F02, F03, F04, and F13. The rms in each image varies due to different amounts of ISM contamination as well as velocity widths averaged, with values being, respectively, 26, 21, 25, 23, and 29 mJy/beam km/s. The peak signal-to-noise ratio ranges between ∼12 – 34.

Fig. C.1. ALMA band 6 CO v = 0 J = 2-1 zeroth moment maps of five red supergiants in RSGC1: F01 (upper left), F02 (upper right), F03 (middle left), F04 (middle right), and F13 (bottom left). The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, with respect to the input coordinates for the ALMA observations (see columns 2 and 3 in Table A.1). The ALMA synthesised beam is shown as a white ellipse in the lower right corner of each panel. The contours (in orange) show the CO emission at [3, 5, 10] times the rms value. The cyan circle has a diameter of size 1″ and is centred on the peak of the continuum emission.

For each source, two channel maps are presented, first a zoom onto the central 2$\stackrel{″}{.}$2 region and then a larger map of size 8″ (Figs. C.3 – C.7). The combined use of these maps is needed to assess the ISM contamination.

Fig. C.2. Measurements of the azimuthally averaged flux in annuli of 200-mas thick. The grey line indicates the 3×25 = 75 mJy beam−1 km s−1 cutoff.

Fig. C.3. 12CO v = 0 J = 2-1 channel map of F01. The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, with respect to the peak of the continuum emission (indicated by the cyan cross in the upper panels). The velocity is given in the upper left corner of each panel (in units of km s−1). The cyan circle indicates the circular aperture of 0$\stackrel{″}{.}$75 (diameter) centred on the peak of the continuum emission that was used for extracting the CO(2-1) line profile. The ALMA beam is shown as a white ellipse in the bottom left corner of the bottom left panel. Upper plot: The contours (in orange) show the CO emission at [3, 5]×σrms (see Table A.1). The contours (in yellow) in the last panel show the continuum emission at [60, 90, and 99]% of the peak continuum emission.

Fig. C.4. 12CO v = 0 J = 2-1 channel map of F02. Same as Fig. C.3, but for F02.

Fig. C.5. 12CO v = 0 J = 2-1 channel map of F03. Same as Fig. C.3, but for F03.

Fig. C.6. 12CO v = 0 J = 2-1 channel map of F04. Same as Fig. C.3, but for F04.

Fig. C.7. 12CO v = 0 J = 2-1 channel map of F13. Same as Fig. C.3, but for F13.

For all sources, there is much less (or no) ISM contamination visible on the high velocity side, although there might be a noise background. In particular:

F01: At velocities lower than ∼105 km s⁻¹ the ISM emission starts contributing. In general, the angular extent of a spherically symmetric source with monotonically increasing velocity will be larger around zero rest velocity, and will decrease when going to redder and bluer velocities. For the source F01, we notice extra emission in the channels from 129.4–135.4 km s⁻¹ offset from the centre to the north-east, that is also visible as increase in flux density in Fig. 1. Either this is noise (the blob is at < 4×σ_rms) or it is a genuine source (as, for example, another cluster member, a companion, ...). When fitting the CO(2-1) line profile for retrieving the gas mass-loss rate, this extra emission is neglected. In the case that F01 would still contribute significant flux in the 129.4-135.3 km s⁻¹ channels, the terminal velocity would be underestimated. Under the rough assumption that the extra blob to the north-east would contribute ∼2.8 mJy in the 129.4 – 135.5 km s⁻¹ velocity range, the terminal velocity could reach up to ∼16 km s⁻¹. A good prediction of the (flux-corrected) red wing and peak flux around 135 km s⁻¹ (with peak flux then being around 5 mJy) is achieved for a mass-loss rate of 2.2×10⁻⁶ M_⊙/yr and v_LSR around 127 km s⁻¹, but the peak around 115 – 125 km s⁻¹ is then underpredicted by a factor of ∼2 (see dotted black line in upper left panel in Fig. 1). That value for Ṁ_CO is only ∼10% higher than the value listed in Table 1.

F02: The CO(2-1) spectrum of F02 is quite clean, with the exception of the ISM contamination at velocities lower than ∼105 km s⁻¹.

F03: For velocities lower than ∼114.4 km s⁻¹, there is a strong contribution from the ISM, which then vanishes at ∼99.4 km s⁻¹. The theoretical line profile for a terminal wind velocity of 13 km s⁻¹ and Ṁ_CO of 3.2×10⁻⁶ M_⊙/yr is slightly underpredicting the red wing of the ALMA CO(2-1) data (Fig. 1). A terminal wind velocity of 17 km s⁻¹ and Ṁ_CO of 3.8×10⁻⁶ M_⊙/yr gives a better prediction of the red wing for flux values below ∼4 mJy, but then the typical two-horn profile is predicted too wide. That difference of 15% in Ṁ_CO owing to an uncertainty in v_∞ is within the quoted uncertainty on the retrieved Ṁ_CO-values, see Table D.1.

F04: Similar to F03, the ISM is contaminating the channel maps and spectrum for velocities ≲114.4 km s⁻¹.

F13: F13 is the source for which the CO(2-1) emission is the least contaminated by the ISM. Only for velocities ≲99.4 km s⁻¹, the line profile might be contaminated.

Appendix D: Uncertainty estimate on CO(2-1) line profile predictions

The radiative transfer analysis of the CO(2-1) emission is based on a number of input parameters. To assess the impact of the model predictions to uncertainties on the different derived or assumed input parameters, we varied the input parameters independently and calculated the percentage impact on the line profile and the associated uncertainty on the mass-loss rate Ṁ_CO (see Table D.1). Uncertainties on the effective temperature (of ±130 K), luminosity and hence stellar radius R_⋆ (±215 R_⊙), and distance D (±890 pc) are taken from Davies et al. (2008), on the fractional carbon abundance and hence on [CO/H] (of ±1.4×10⁻⁵) from Davies et al. (2009), on the terminal wind velocity (of ±3 km/s) and the outer radius (of ±50 R_⋆) from our ALMA data, and on the coefficients for the wind acceleration profile β (±2) and for the temperature profile ϵ (±1) from the analysis presented in Decin et al. (2006) and De Beck et al. (2010). The canonical gas-to-dust ratio of 200 was taken representative for a Milky Way cluster (Beasor et al. 2020), but is quite ill-constrained. We therefore let it range between 100 – 1000. 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%). Uncertainties on all other parameters only have a minor effect. We note that the uncertainty estimates in Table D.1 are conservative values since only based on reproducing the peak flux, but not accounting for the full shape of the line profile which carries important diagnostics on, for example, the optical thickness of a specific line.

To calculate the combined effect of the errors in the individual input parameters on Ṁ_CO, one should perform a Monte Carlo Markov Chain (MCMC) analysis. However, this is not feasible since non-LTE radiative tranfser calculations are too computer-intensive to be embedded in a MCMC algorithm. We therefore need to resort to a simpler approach similar to, for example, Díaz-Luis et al. (2019) who have used an analytical prescription derived by Knapp & Morris (1985) to calculate the mass-loss rate given some individual parameters such as the expansion velocity, distance, integrated line intensity, and [CO/H] ratio (see, for example, Eq. 3 in Díaz-Luis et al. 2019). However, we here use our model grid and results presented in De Beck et al. (2010) since based on the same GASTRONOOM models, in particular we use Eq. 9 of De Beck et al. (2010). The CO envelope size, R_out, is not an individual input value of Eq. 9 of De Beck et al. (2010), but is inherently captured in Eq. 9 since the [CO/H] ratio determines the photodissociation radius of the circumstellar envelope. Contrary to Díaz-Luis et al. (2019), we do not add the individual uncertainties in (unweighted) quadrature, but use a proper error propagation based on Eq. 9 of De Beck et al. (2010) which weighs the sensitivity of each individual uncertainty. Hence, we derive that the total uncertainty on the estimated Ṁ_CO-values is a factor of ∼1.4. In that, the uncertainty in terminal wind velocity contributes some 50% of the total uncertainty factor and the distance ∼30%.

It might seem surprising that this total uncertainty on Ṁ_CO is similar to the highest uncertainty raised by an individual parameter (of ∼42%) listed in Table D.1. The reason thereof is that Table D.1 uses the sensitivity of the intensity peak to determine the associated change in Ṁ_CO by using a full non-LTE radiative transfer routine, while Eq. 9 of De Beck et al. (2010) is an analytic relation based on the integrated line intensity with various weightings (in terms of exponential factors) for the input parameters. Those different approaches hence rule out a one-to-one comparison of both uncertainty values.

Appendix E: Ṁ-luminosity relations from Beasor et al. (2020)

In a first step, Beasor et al. (2020) has determined a Ṁ-luminosity relation for all clusters in their sample by fitting the relation given in Eq. (4) to their data points taken from Beasor & Davies (2016), Beasor & Davies (2018) and Beasor et al. (2020) (see Table E.1). They therefore have used a linear least-square approximation that accounts for both errors on the ordinate and coordinate values. However, plotting the Ṁ-luminosity relation using the values for the intercept and slope as listed in Table 4 of Beasor et al. (2020) (repeated here in Table E.2) reveals a problem with those values (see left panel in Fig. E.1): (1) A typo occurred in Table 4 of Beasor et al. (2020): the intercept a for RSGC1 should not be −52, but −53 (Beasor, priv. comm.), which indeed yields a much better fit to the data (see dotted line in left panel of Fig. E.1); and (2) comparing Fig. 3 and Fig. E.1 with Fig. 2 in Beasor et al. (2020) shows a problem with the data points plotted in Beasor et al. (2020), for example the maximum Ṁ_SED value for NGC 7419 is 1.8×10⁻⁵ M_⊙/yr (see Table E.1) but is plotted around 7×10⁻⁵ M_⊙/yr in Fig. 2 of Beasor et al. (2020); similar problems occur for other data points. For that reason we continue with the data listed in Table E.1; a corrigendum is currently being prepared by E. Beasor et al.

Fig. E.1. Ṁ-luminosity relations as determined by Beasor et al. (2020). The coloured symbols represent the (Lbol, Ṁ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(Lbol/L⊙) for all clusters in the sample. In the left panel, the dashed lines show the Ṁ-Lbol 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 Ṁ-Lbol 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 Table 1 of Beasor 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).

In a next step, Beasor et al. (2020) chose to fix the gradient in order to reduce the number of degrees of freedom in the fit. Accounting for the standard deviation on the derived b-values, one gets b = 4.8±0.6. The latter value was then used for deriving the intercept values in their more general mass-dependent Ṁ-luminosity relation given in Eq. (3) (which is Eq. (4) in Beasor et al. 2020). But for RGSC1, the Ṁ-luminosity relation from their Eq. (4) seems at odd with the empirical (L_bol, Ṁ_SED) values (see right panel in Fig. E.1), with the straight-line fit being systematically lower than the empirical values. We therefore have repeated their analysis for a slope of b = 4.8±0.6. Fixing the slope, only the intercept has to be determined; the derived values with their standard deviation are listed in the fourth column of Table E.2. From the analysis of the four clusters, we can determine the mass-dependent intercept, yielding the following mass-dependent Ṁ_SED-luminosity relation

$$\begin{array}{cc}\hfill log({\dot{M}}_{\mathrm{SED}}/M_{\odot }{\mathrm{yr}}^{-1})& =(-26.5-0.22\times M_{\mathrm{ini}}/M_{\odot })+4.8log(L_{\mathrm{bol}}/L_{\odot })\hfill \end{array}$$(E.1)

with standard deviations on the best-fit parameters being −26.50±1.33, −0.22±0.10, and 4.8±0.6; which is similar to Eq. (3) within some round-off errors. Eq. (E.1) implies a slight shift upwards for RSGC1 as compared to Eq. (3) (see red dashed-triple dotted line in the right panel of Fig. E.1), but that shift is not visible on the scale of the plot for the 3 other clusters.

Appendix F: Sensitivity to intercept a_i

Fig. F.1 shows the fit to Eq. (5), but instead of using RSGC1 for determining a₁, we use NGC 2100 for determining the intercept a₂ and the difference in intercept of the other 3 clusters, Δa_i (i = 1,3,4). We fit both the (L_bol, Ṁ_SED) and (L_bol, Ṁ_CO) together; corresponding values are listed in Table F.1.

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 (Lbol, ṀSED) measurements of the four clusters (shown as open symbols with corresponding colour), and the dashed purple line the best-fit to the (Lbol, ṀCO) values of four stars in RSGC1 (filled purple downward triangles).

Appendix G: Other Galactic red supergiants

In this section, we derive Ṁ_CO for some well-known Galactic M-type red supergiants following the same modelling strategy as in Sect. 3.1. The aim is to compare the retrieved Ṁ_CO with the values predicted using Eq. (8). We therefore perform a radiative transfer analysis of the CO rotational line transitions of the red supergiants published by De Beck et al. (2010). In particular, we analyse the James Clerk Maxwell Telescope (JCMT) data of α Ori, μ Cep, and VX Sgr. The M2 red supergiant VY CMa is excluded from the current analysis, since the CO rotational line data have already been analysed by Decin et al. (2006) and we will resort to their outcomes. We also exclude the hypergiant NML Cyg, since its estimated initial mass of ∼40 M_⊙ (Schuster et al. 2009) puts its above the mass range of 8≲M_⋆≲30 M_⊙ of interest to this study. The JCMT data of α Ori, μ Cep, and VX Sgr are plotted in Fig. G.1. The intensity is plotted as main beam temperature T_mb, in units of K, which can be converted to flux density F_ν or brightness I_ν by following formula

$$\begin{array}{c}\hfill T_{\mathrm{mb}}=\frac{{\lambda }^2{F}_{\nu }}{2k\mathrm{\Omega }}=1.36{\left(\frac{{\lambda }_{\mathrm{cm}}}{\theta ″}\right)}^2{I}_{\nu },\end{array}$$(G.1)

Fig. G.1. CO rotational line transitions observed towards α Ori (top panels), μ Cep (middle panels), and VX Sgr (bottom panels) with the JCMT. The CO(2-1), CO(3-2), and CO(4-3) data are plotted as red histograms. Synthetic line profiles are overplotted as black solid lines. The vertical dashed grey lines indicate the LSR velocity. Accounting for the absolute flux uncertainty of ∼30%, the JCMT data have been scaled by a factor, indicated in the upper right corner of each panel.

with λ the wavelength, k the Boltzmann constant, Ω the beam solid angle, λ_cm the wavelength in units of cm, θ the beam size in units of arcsec, and I_ν the brightness in units of mJy/beam. The beam size of the JCMT CO(2-1), CO(3-2), and CO(4-3) data are 19$\stackrel{″}{.}$7, 13$\stackrel{″}{.}$2, and 10$\stackrel{″}{.}$8, respectively.

Similar to the analysis of the five red supergiants in RSGC1 (Sect. 3.1), we assume that the main wind component can be represented by a spherically symmetric envelope. As has been discussed in detail in the literature, high spatial and spectral resolution observations of those Galactic RSGs prove that those winds are composed of a smooth outflow and individual mass ejections in the form of gas and dust clumps (Decin et al. 2006; Montargès et al. 2019, 2021; Humphreys & Jones 2022). The smooth outflow contains by far most of the mass and its mass-loss rate is a factor of a few up to almost 2 orders of magnitude higher than that traced by the individual clumps (Montargès et al. 2019, 2021; Humphreys & Jones 2022). With the aim of (i) deriving the overall (average) mass-loss rate Ṁ_CO, and (ii) comparing Ṁ_CO with the predictions from Eq. (8), we only focus on modelling the dominant outflow. This implies that spikes, visible on top of the JCMT CO line profiles and testifying the presence of individual gas clumps, will not be modelled.

To retrieve Ṁ_CO, we follow the same modelling strategy as outlined in Sect. 3.1. The input stellar and CSE parameters are listed in Table G.1. For α Ori and μ Cep, the CO outer envelope size can be estimated from radio interferometric observations, while for VX Sgr we use the photodissociation radius following the results of Groenewegen & Saberi (2021). For each of the sources, the individual modelling could be optimised by fine-tuning some of the input parameters described in Sect. 3.1. However, we refrain from doing so since we aim to assess the predictive power of Eq. (8) in relation to our overall modelling strategy of retrieving Ṁ_CO. Given the complexity of the line profiles and the absolute flux calibration accuracy of ∼30% (Kemper et al. 2003), the fits shown in Fig. G.1 prove that for all sources the main CO component is well reproduced. In contrast to the analysis of the RSGC1 RSGs, the main uncertainty on the retrieved Ṁ_CO values arises from uncertainties on the distance; translating into an Ṁ_CO uncertainty of a factor of ∼2 – 3.

The gas mass-loss rates retrieved from this analysis, Ṁ_CO, are listed in Table G.1 and can be compared to the Ṁ-values predicted using Eq. (8) and using Eq. (3), the latter once following the results of Beasor et al. (2020). With the exception of the cool extreme supergiant VY CMa, Eq. (8) predicts the gas mass-loss rate for the other 3 red supergiants well — the average difference only being ∼30%. This result supports the use of Eq. (8) for the prediction of mass-loss rates for M-type red supergiants with effective temperature ≳3200 K. Eq. (3) predicts the gas mass-loss rates well for μ Cep and VX Sgr, but fails for α Ori and VY CMa.

The derived Ṁ_CO values can be compared to Ṁ_SED retrievals to determine the gas-to-dust ratio of the winds of those Galactic RSGs. Assuming a wind speed of 10 km/s and gas-to-dust ratio of 100, Verhoelst et al. (2009) derived that Ṁ_SED equals 6.3×10⁻⁸ M_⊙/yr (for a distance of 131 pc) in the case of α Ori, and 2.9×10⁻⁷ M_⊙/yr (for a distance of 830 pc) in the case of μ Cep. Similarly we use the values derived by Liu et al. (2017) and Shenoy et al. (2016) for VX Sgr and VY CMa, respectively. Scaling for the difference in adopted distance and wind speed, we derive a gas-to-dust ratio of 297, 544, 226, and 18 for α Ori, μ Cep, VX Sgr, and VY CMa, respectively. The gas-to-dust ratio derived for VY CMa might seem very low. However, it should be noted that Decin et al. (2006) derived that VY CMa underwent a phase of high mass loss (∼3.2×10⁻⁴ M_⊙/yr; see Table G.1) some 1000 yr ago. This phase took some 100 yr and was preceded by a low mass-loss phase (∼1×10⁻⁶ M_⊙/yr) taking some 800 yr. The current mass-loss rate (of ∼8×10⁻⁵ M_⊙/yr) was used to derive the gas-to-dust value, which obviously will change during other mass-loss phases. Moreover, the assumption of a spherical symmetry for the circumstellar envelope of VY CMa is obviously a simplification on the basis of which the dust mass-loss rate has been derived by Shenoy et al. (2016); see discussion in that paper.

All Tables

All Figures

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 Appendix A). The vertical dashed blue lines indicate the local standard of rest velocity, vLSR, as deduced from the ALMA CO data (see Appendix A). The dotted and dashed-dotted blue lines indicate the vLSR value as determined by 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 Appendix C.

In the text

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 of Kee et al. (2021) is shown with a dashed line for 2 different values of the atmospheric turbulent velocity $v_{\mathrm{turb}}^{\mathrm{atm}}$. 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 and 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.

In the text

Fig. 3. ṀSED–luminosity relations for the four open clusters RSGC1, NGC 2100, NGC 7419, and χ Per. The coloured open symbols represent the (Lbol, Ṁ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 Appendix 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 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.

In the text

Fig. 4. Ṁ–luminosity relations parametrised in terms of log Lbol 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 Lbol values of Davies et al. (2008).

In the text

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 105 years; tracks go till carbon depletion indicated by the open square. The inset at the RSG, for a range in log Teff of 0.01 dex, shows how both tracks digress there.

In the text

Fig. B.1. ALMA band 6 SiO v = 3 J = 5-4 zeroth moment maps of five red supergiants in RSGC1: F01 (upper left), F02 (upper right), F03 (middle left), F04 (middle right), and F13 (bottom left). The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, with respect to the input coordinates for the ALMA observations (see columns 2 and 3 in Table A.1). The ALMA synthesised beam is shown as a white ellipse in the lower right corner of each panel. The contours (in orange) show the SiO emission at [3, 5, 10] times the rms value.

In the text

Fig. B.2. ALMA band 6 continuum maps of five red supergiants in RSGC1: F01 (upper left), F02 (upper right), F03 (middle left), F04 (middle right), and F13 (bottom left). Contours (in orange) are indicated at [50, 60, 90, and 99]% of the peak continuum emission (see Table A.1). The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, with respect to the peak of the SiO v = 3 J = 5-4 zeroth moment map. The ALMA synthesised beam is shown as a white ellipse in the lower right corner of each panel.

In the text

Fig. C.1. ALMA band 6 CO v = 0 J = 2-1 zeroth moment maps of five red supergiants in RSGC1: F01 (upper left), F02 (upper right), F03 (middle left), F04 (middle right), and F13 (bottom left). The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, with respect to the input coordinates for the ALMA observations (see columns 2 and 3 in Table A.1). The ALMA synthesised beam is shown as a white ellipse in the lower right corner of each panel. The contours (in orange) show the CO emission at [3, 5, 10] times the rms value. The cyan circle has a diameter of size 1″ and is centred on the peak of the continuum emission.

In the text

	Fig. C.2. Measurements of the azimuthally averaged flux in annuli of 200-mas thick. The grey line indicates the 3×25 = 75 mJy beam−1 km s−1 cutoff.
In the text

Fig. C.3. 12CO v = 0 J = 2-1 channel map of F01. The ordinate and co-ordinate axis give the offset of the right ascension and declination, respectively, with respect to the peak of the continuum emission (indicated by the cyan cross in the upper panels). The velocity is given in the upper left corner of each panel (in units of km s−1). The cyan circle indicates the circular aperture of 0$\stackrel{″}{.}$75 (diameter) centred on the peak of the continuum emission that was used for extracting the CO(2-1) line profile. The ALMA beam is shown as a white ellipse in the bottom left corner of the bottom left panel. Upper plot: The contours (in orange) show the CO emission at [3, 5]×σrms (see Table A.1). The contours (in yellow) in the last panel show the continuum emission at [60, 90, and 99]% of the peak continuum emission.

In the text

	Fig. C.4. 12CO v = 0 J = 2-1 channel map of F02. Same as Fig. C.3, but for F02.
In the text

	Fig. C.5. 12CO v = 0 J = 2-1 channel map of F03. Same as Fig. C.3, but for F03.
In the text

	Fig. C.6. 12CO v = 0 J = 2-1 channel map of F04. Same as Fig. C.3, but for F04.
In the text

	Fig. C.7. 12CO v = 0 J = 2-1 channel map of F13. Same as Fig. C.3, but for F13.
In the text

Fig. E.1. Ṁ-luminosity relations as determined by Beasor et al. (2020). The coloured symbols represent the (Lbol, Ṁ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(Lbol/L⊙) for all clusters in the sample. In the left panel, the dashed lines show the Ṁ-Lbol 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 Ṁ-Lbol 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 Table 1 of Beasor 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).

In the text

	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 (Lbol, ṀSED) measurements of the four clusters (shown as open symbols with corresponding colour), and the dashed purple line the best-fit to the (Lbol, ṀCO) values of four stars in RSGC1 (filled purple downward triangles).
In the text

Fig. G.1. CO rotational line transitions observed towards α Ori (top panels), μ Cep (middle panels), and VX Sgr (bottom panels) with the JCMT. The CO(2-1), CO(3-2), and CO(4-3) data are plotted as red histograms. Synthetic line profiles are overplotted as black solid lines. The vertical dashed grey lines indicate the LSR velocity. Accounting for the absolute flux uncertainty of ∼30%, the JCMT data have been scaled by a factor, indicated in the upper right corner of each panel.

In the text