Red Supergiants in M31 and M33. II. The Mass-loss Rate

Mass loss is an important activity for red supergiants (RSGs) and can influence their evolution and final fate. Previous estimations of mass-loss rates (MLRs) of RSGs exhibit significant dispersion due to differences in method and the incompleteness of samples. With the improved quality and depth of surveys including the UKIRT/WFCAM observations in the near-infrared, and LGGS and PS1 in the optical, a rather complete sample of RSGs is identified in M31 and M33 according to their brightness and colors. For about 2000 objects in either galaxy from this largest ever sample, the MLR is derived by fitting the observational optical-to-mid-infrared spectral energy distribution with the DUSTY code of a 1D dust radiative transfer model. The average MLR of RSGs is found to be around 2.0 × 10−5 M ⊙ yr−1 with a gas-to-dust ratio of 100, which yields a total contribution to the interstellar dust from RSGs of about 1.1 × 10−3 M ⊙ yr−1 in M31 and 6.0 × 10−4 M ⊙ yr−1 in M33, a non-negligible source in comparison with evolved low-mass stars. The MLRs are divided into three types by the dust species, i.e., amorphous silicate, amorphous carbon, and optically thin, and the relations between MLR and stellar parameters, infrared flux, and colors are discussed and compared with previous works for the silicate and carbon dust groups.


Introduction
Mass loss is a common phenomenon of red supergiants (RSGs), evidenced by the observed infrared excess first in NML Cyg and VY CMa (Johnson 1968;Hyland et al. 1969) and lately in numerous RSGs. The destiny of an RSG depends on mass-loss rate (MLR) in addition to initial mass and chemical composition. Some RSGs may stay in the RSG stage and eventually explode as hydrogen-rich Type II-P supernovae (SNe), while others may evolve backward to the blue end of the H-R diagram and spend some short time as yellow supergiant stars, blue supergiant stars, or Wolf-Rayet stars before exploding as SNe (Smartt 2009;Humphreys 2010;Ekström et al. 2013;Meynet et al. 2015;Beasor & Davies 2018). In either case, the strong mass loss during the RSG phase would greatly influence its ultimate fate. Moreover, the mass loss of RSGs significantly contributes to interstellar dust, particularly in young galaxies at high redshift, where the potential dust producers-asymptotic giant branch stars (AGBs) are not yet evolved (Massey et al. 2005;Levesque 2010).
There have been various types of methods to estimate the MLR (expressed as mass loss per year, M e yr −1 ) of RSGs from their infrared emission, which increases with the amount of dust. One type fits the spectral energy distribution (SED) by a dust radiative transfer model (see, e.g., Gordon et al. 2018). Alternatively, spectral analysis may be used to derive the MLR with the shell expansion velocities determined from the spectral line, e.g., an OH double-peak maser (Goldman et al. 2017). The classical Reimers law (Reimers 1975;Kudritzki & Reimers 1978) was derived from high-resolution spectroscopy from which the outflow velocities of stellar winds and the optical depths of circumstellar lines were determined. There are some empirical relations between the MLR and stellar parameters such as luminosity or effective temperature, which are usually based on the results from the SED fitting.
RSGs in the Milky Way are the nearest, though they suffer heavy interstellar extinction in the Galactic plane. Still, some of them are close enough to be observed in multiple bands or even spatially resolved. The calculation leads to high MLR for some famous RSGs. Gordon et al. (2018) estimated the MLR of VX Sgr, NML Cyg, and S Per from infrared imaging to be 2-3 × 10 −5 M e yr −1 on average. Beasor & Davies (2018) studied the RSGs in the clusters NGC 7419 and χ Per, yielding an MLR from 10 −7 M e yr −1 to 10 −5 M e yr −1 and found a steep increase of MLR with stellar luminosity. For a sample of about 15 dustenshrouded RSGs in the Large Magellanic Cloud (LMC), van Loon et al. (2005) derived their MLR to lie mostly in the range 10 −4 to 10 −5 M e yr −1 by using the DUSTY code, a dust radiative transfer model (Ivezic & Elitzur 1997). For a large sample (∼30,000) of AGBs and RSGs in the LMC, Riebel et al. (2012) used the GRAMS grids and obtained a dust MLR mostly between 10 −7 and 10 −11 M e yr −1 , much lower than found by van Loon et al. (2005). Gordon et al. (2016) studied the RSGs in M31 and M33 and found that the massive RSGs have an average MLR of 10 −4 to 10 −5 M e yr −1 . Groenewegen & Sloan (2018) estimated an MLR between 10 −4 and 10 −10 M e yr −1 for RSGs and AGBs in the LMC, SMC, and dwarf spheroidal galaxies in the Local Group.
The large range of MLR in previous studies (from about 10 −11 to 10 −4 M e yr −1 ) originated from differences in the sample of RSGs, the dust properties, the models, and so on. Previous studies selected only a small number of RSGs, mostly bright in the optical and/or infrared, which would incline to cases of high MLR and cannot reflect the MLR of all RSGs unbiasedly. The differences in method also bring about difficulties in comparing the results. Ren et al. (2021, hereafter Paper I) constructed a sample of RSGs in M31 and M33 based mainly on the near-infrared color-color diagram and Gaia astrometric information, which is claimed to be complete. This provides us the opportunity to investigate the MLR of all the RSGs in an individual galaxy. Besides, both M31 and M33 are spiral galaxies with similar metallicity to our Milky Way, whose entire scenario cannot be viewed from within it; this would shed light on the MLR, the properties of dust from RSGs, and its contribution to the interstellar medium in our Galaxy. This work intends to calculate the MLR for a couple of thousand RSGs in M31 and M33 in this newly selected sample, the majority of which are identified as RSGs for the first time, and to discuss the properties of circumstellar dust and the contribution to interstellar dust from RSGs.
The paper is organized as follows. Section 2 describes the sample, and Section 3 presents the dust model we use. The last section is about the result and some discussion on MLRs of RSGs.

Preliminary Sample
The initial sample of RSGs in M31 and M33 is taken from Paper I. Massey et al. (2006Massey et al. ( , 2007 selected 437 and 749 RSG candidates in M33 and M31 respectively by their optical spectra and colors. In comparison with the number (∼1400) of RSGs in the SMC identified by Yang et al. (2018Yang et al. ( , 2019, the optically selected sample in M31 and M33 is far from complete. Paper I uses the deep UKIRT/WFCAM photometry in the JHK bands and an innovative method to exclude foreground stars, i.e., to remove foreground dwarf stars from their obvious branch in the J − H/H − K diagram due to the effect of their high surface gravity on the H band. Afterwards, RSGs are recognized by their large brightness and low effective temperature, and in practice from their location in the J − K/K diagram that coincides well with the theoretical evolutionary track of a massive star. Consequently, Paper I identified 5498 and 3055 RSGs in M31 (along with M32 and M110) and M33 respectively, several times the numbers in the previous samples. After excluding sources outside the areas of M31 and M33, there are 5253 and 3001 RSGs left, which form the initial sample of this work. Details of the selection and discussion on the completeness of the sample can be found in Paper I.

Final Sample: Combination of Catalogs in Various Wave Bands
The initial sample is further reduced to be appropriate for estimating MLR. In complying with the traditional method, we choose to fit the SED by the dust radiative transfer model. In order to characterize the SED of an RSG, it is better to have the photometric results covering wider wave bands (in particular in the infrared) with higher accuracy. In principle, the optical bands can constrain the properties of the central star while the infrared bands demonstrate the radiation of the dust. The optical photometric data are taken from the Pan-STARRS1 large-scale survey in the g, r, i, z, and y bands (Chambers et al. 2016) and the dedicated Local Group Galaxies Survey (LGGS) in the U, B, V, R, and I bands (Massey et al. 2006(Massey et al. , 2007. In addition to the UKIRT/WFCAM photometry in the J, H, and K bands available for all the objects, the photometry at longer wavelengths is taken from the Spitzer/IRAC (IRAC1, IRAC2, IRAC3, and IRAC4) (Werner et al. 2004), Spitzer/MIPS24 (Rieke et al. 2004), and Wide-field Infrared Survey Explorer (WISE) (W1, W2, W3) (Wright et al. 2010) observations. The WISE/W4 photometry is not used because its flux is usually much higher than that of the Spitzer/MIPS24, which may be caused by the low resolution in this band. These data consist of an SED from optical-violet (366 nm) to mid-infrared (24 μm) that covers the silicate features around 9.7 and 18 μm. Though Gordon et al. (2018) reported the detection of emission at Herschel/70 and 160 μm for a few RSGs in the Milky Way, the photometry at such long wavelengths is unavailable for RSGs in M31 and M33 at a distance modulus as large as ∼24.5 mag. Fortunately, the very high luminosity of RSGs ensures that the circumstellar dust is hot or warm, so the lack of longer wavelength photometry has little influence on determining the MLR. On the photometric error, 0.1 mag is set for the PS1 force point-spread function magnitude (photometric measurements conducted at a specified location) and the LGGS and UKIRT magnitudes, while the Spitzer and WISE data are limited to be better than 3σ detection. Moreover, the source is required to have eligible photometry in at least three PS1 or LGGS bands, and an IRAC or a WISE band. The final sample consists of 1741 RSGs in M31 and 1983 RSGs in M33 from the initial sample that satisfy the above conditions.
The sample used to calculate the MLR is compared with the preliminary sample in the J − K/K diagram in Figure 1 for M31 and in Figure 2 for M33. The histogram in J − K and K exhibits no apparent systematic bias from the initial sample, and the MLR derived can then be regarded as representative of the whole sample.

Interstellar Extinction
The interstellar extinction is subtracted before estimating the MLR. For this purpose, the SFD98 (Schlegel et al. 1998) extinction map of the M31 and M33 regions is adopted. Though there are some arguments that the SFD98 map may underestimate or overestimate the interstellar extinction along some sightlines, it has some advantages. The SFD98 extinction is derived from the infrared emission that comes not only from the foreground Milky Way but also from M31 or M33 itself, so consequently the derived extinction includes both sources. With E(B − V ) from SFD98, the extinction in the U to K band is calculated with the CCM89 law (Cardelli et al. 1989) at R V = 3.1, and the extinction in the infrared bands is calculated with the law of Wang & Chen (2019) and Xue et al. (2016). These laws are representative of the average Galactic extinction, but may also be plausible for M31 and M33 (see, e.g., Bianchi et al. 1996) because they are spiral galaxies like the Milky Way; needless to say, the extinction law of M31 and M33 is not yet well determined. The distribution of the adopted A V for the RSGs is shown in Figure 3. It can be seen that A V to the RSGs in M31 ranges from about 0.3 to 1.5 mag, while it is much smaller in M33, i.e., from 0.1 to 0.2 mag, and both are higher than the foreground Galactic extinction, which is 0.17 mag and 0.11 mag for M31 and M33 respectively. In comparison, Massey & Evans (2016) adopt a constant A V = 1, which overestimates the extinction on average. It can be imagined that the overestimation of interstellar extinction would decrease the MLR. Quantitatively, we examine the influence of interstellar extinction by running the procedure for two cases with constant A V = 0 and A V = 1.0 respectively. It is found that the mean MLR would increase by 1% at A V = 0 and decrease by 5% at A V = 1.0.

Effective Temperature and Luminosity
The effective temperature (in kelvin) is calculated by using the following equation from Neugent et al. (2020) with the intrinsic color index (J − K ) 0 : The luminosity of RSGs is then calculated from the K-band magnitude and its bolometric correction with the distance moduli of 24.32 for M31 and 24.50 for M33 (Paper I), where the bolometric correction is also taken from Neugent et al. We choose the DUSTY code for its flexibility of input parameters, and its correctness has been confirmed by many applications. Because our sample of RSGs is large, we follow a way similar to GRAMS whereby the model grids are constructed in advance, covering a wide range of dust parameters and stellar models, then the model that best matches the observational SED is chosen. In comparison with GRAMS, our grids have higher resolution in the parameters space so that a better fitting can be obtained.
For a dust shell expanding with a constant velocity and consequently a density distribution of ρ ∝ r −2 , the MLR can be calculated from the dust and stellar parameters as follows (Beasor & Davies 2018): where a Q e / refers to the ratio of the dust radius a to its extinction efficiency Q e , ρ d is the dust density, v the expansion velocity, r 0 the inner radius of the dust shell, and τ V the optical depth of the dust shell at 550 nm. The parameter f refers to the gas-to-dust ratio.
The wind velocity of RSGs is difficult to measure. By studying AGBs and RSGs, Groenewegen et al. (2002), Ramstedt et al.  carbon-rich AGB stars respectively, which means there is little difference between them. In this work, we choose the formula of Goldman et al. (2017) for all dust species (silicate or carbonaceous dust) independent of the RSG. As the metallicity in M31 and M33 is similar to that in the Milky Way, the solar metallicity is used to calculate the wind speed, which leads to = By substituting the luminosity derived above, the resultant expansion velocities vary from 4 km s −1 to 15 km s −1 , which  agrees with the previous works. Since the MLR is proportional to the wind speed for a static wind, it would be easy to convert for other velocities.
As the DUSTY model (Ivezic & Elitzur 1997) generates the output SED at a given luminosity of 10 4 L e , the MLR can be scaled as / is a constant for the given dust species and size, and τ V and r 0 are the parameters to be derived from the model fitting. Finally, the MLR depends on three parameters for the given dust properties: optical depth τ V , inner radius r 0 , and luminosity factor F.

The Models Library
The model grids are built in a wide range of dust and stellar parameters as follows.
1. The dust species is either silicate, for which the optical constants from Draine & Lee (1984) is adopted with ρ = 3.5 g cm −3 , or amorphous carbon for which the optical constant from Rouleau & Martin (1991; RM91) is adopted with ρ = 2.26 g cm −3 . RSGs are massive stars and originally oxygen-rich, complying with cosmic abundance since no deep dredge-up as in AGB stars occurs. The silicate spectral features are already detected in many RSGs, e.g., in RS Per by Gordon et al. (2018) and in a group of RSGs by Verhoelst et al. (2009  If a constant radius is taken as is done in some previous works, this parameter will change. For a constant radius of 0.1 μm, this ratio changes to 0.139 for amorphous silicate and 0.052 for amorphous carbon, which means the MLR would decrease by 91% and 80%, respectively, i.e., by an order of magnitude. Definitely this size effect should be taken into account when comparing the MLR from various works. The change in this ratio is displayed in Figure 4 for the popular values. 4. The dust temperature at the inner radius has four options: 1200, 1050, 900, and 750 K. The condensation temperature for silicate and amorphous carbon depends on the pressure and other factors and thus varies. Gail & Sedlmayr (1984, 1999 found this temperature to be around 1200 K in the circumstellar envelope around latetype star. The dust shell may be further away from the condensation line so that a lower temperature is possible. Since far-infrared data are not available, a temperature lower than 750 K is not considered, which may underestimate the MLR, but this should not be significant for luminous RSGs. 5. The optical depth has 100 options from 0.001 to 10 in approximately equal logarithmic intervals. 6. The stellar atmosphere model comes from the MARCS code (Gustafsson et al. 2008). We select the RSG models at the given mass of 15 M e , metallicity of where Z is the metallicity mass fraction), and surface gravity of = g log 0, which are typical for an RSG in M31 or M33. Meanwhile, the effective temperature has 10 options: 3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000, 4250, and 4500 K. The influence of metallicity and surface gravity will be discussed later. The observational SED needs some models with an even lower effective temperature, which is not available for massive stars, so we take a 1 M e red-giant-like model at 2800, 2900, 3000, and 3200 K as the substitute. This is practically feasible because the SED resolution is too low to distinguish the spectra of RSGs and red giant branch stars at this low temperature. Moreover, only a very small proportion (1%) of sources have such a low temperature.

Fitting
The best model in the library is chosen by calculating the smallest χ 2 . After the flux is normalized to the K-band flux, which is available for all sources and least affected by the dust shell among the infrared bands, χ 2 is calculated as where F(M, λ) and F(M, K ) refer to the model flux, and F(O, λ) and F(O, K ) to the observational flux, in the λ and K bands respectively, and N|f (O, λ)| is the number of observational points. The χ 2 distribution is shown in Figure 5 after dropping 12 sources with χ 2 > 0.5. The unsuccessful cases are mainly caused by photometric error. In addition the wrong crossidentification may be responsible for a few cases, and variation in the light from RSGs may bring about the disparity between optical and infrared fluxes. Some RSGs have optically thin shells with τ V < 0.1, for which no discrimination can be made on dust species. For such sources, the amorphous silicate is assigned according to the cosmic abundance, and they constitute an "optically thin" group to be treated separately. The sources with optical depth τ V 0.1 are divided into two groups by dust species, either amorphous silicate or carbon determined by which yields the smaller χ 2 . Figure 6 displays examples of four typical cases of circumstellar dust, i.e., amorphous silicate, amorphous carbon, optically thin, and one unsuccessful case. The parameter T eff derived from the model fitting is compared with that from the color index J − K (Equation (1)) in Figure 7. Because the model grids provide only some discrete values, many sources have the same T eff , thus the mean and the standard deviation are marked. The agreement is very good at T eff > 3400 K, while T eff becomes smaller than that from J − K at T eff < 3400 K. This discrepancy may be attributed to replacing the massive star model by a low-mass star model at low T eff , but might also arise because Equation (1) for transforming the color index to T eff is not valid at low T eff .

Proportion of RSGs with Silicate and Carbon Dust
The results of fitting are divided into three types according to the species of circumstellar dust: (1) silicate, (2) carbon, and (3) optically thin with τ V < 0.1. In general, the silicate emission feature around 9.7 μm leads to a relatively steep increase from the Spitzer/IRAC3 band at ∼5.8 μm to the IRAC4 band at ∼8.0 μm, while the carbonaceous dust produces a plateau-like shape from 3.5 to 8 μm. For the third type, the optical depth of the circumstellar shell is very small such that the infrared SED can be almost described by the stellar model; consequently there is no significant difference in taking the amorphous silicate or carbon dust. Because RSGs are supposed to be O-rich in origin, we assign Type 3 objects to have silicate dust as Type 1.
The results for all the objects are listed in Tables 1 and 2. The number of objects in each type is presented in Table 3. Among the 3712 RSGs (1733 + 1979 for M31 and M33 respectively), 1665 sources (671 + 994) are best fitted by silicate dust and 581 sources (329 + 252) are attributed to amorphous carbon dust. This result means that 16% of RSGs have carbon-rich circumstellar dust on average, with 19% in M31 and 13% in M33. These percentages represent the   statistical study of the proportion of carbon-rich dust around RSGs for the first time. It may be argued that the SED of some assigned carbon-rich RSGs can also be fitted by silicate dust. This is true if the priority is given not to the smaller χ 2 but to the silicate dust. On other hand, our fitting is not skewed to carbon dust. Infrared spectroscopy from facilities such as the James Webb Space Telescope would be very helpful to ascertain the chemical nature. Indeed, there have been other attempts to identify the carbon dust around RSGs from photometry. For example, Yang et al. (2018) attributed the difference in the mid-infrared colors related to the Spitzer and WISE bands to the PAH emission for the RSGs in the LMC. Basically, PAHs are an indicator of a carbonaceous environment. Interestingly, they found about 15% of RSGs showing evidence of PAHs; this proportion coincides with the present work. Similarly, Yang et al. (2020) studied the SMC RSGs and also found the mid-infrared colors to be an indicator of the PAH emission, but only 4% being so, a much smaller proportion than in the LMC.
In comparison, the ratio of carbon-rich to oxygen-rich AGB stars is found to be 0.53 by Humphreys (2010) or 0.62 in the inner 4 kpc for NGC 6822, which translates into an iron abundance of [Fe/H] = −1.29 by Sibbons et al. (2012). It seems our derived proportion of carbon-rich RSGs is not very exaggerated if compared with the AGB stars. However, the metallicity of RSGs in NGC 6822 is found to be [Z] ∼ −0.52 by Patrick et al. (2015), which apparently differs from [Fe/ H] = −1.29 derived from the C-/M-AGB ratio (from the number ratio of carbon-rich to oxygen-rich AGB stars). RSGs and AGBs may be different populations due to their distinction in mass. Although low-and intermediate-mass stars experience a deep convection process that can dredge up the carbon produced in helium-burning to the surface to become carbonrich, the mechanism for an RSG to become carbon-rich is lacking, though convection exists at the surface of RSGs as well. The existence of carbon-rich RSGs is still puzzling.
The proportion of RSGs with carbon dust is found to increase with brightness, as shown in Figure 8. Interestingly the brightest RSG in both M31 and M33 has a carbon dust shell. This may be accidental because the number of bright RSGs has no statistical significance, thus the proportion for K < 13 is unreliable. But the trend of increasing carbon dust with brightness continues for fainter RSGs. In other words, the more massive RSGs may possess a carbon-rich dust shell more easily. In the color-color diagrams, Figures 9 and 10, the carbon dust is distinguishable from the silicate dust. The RSGs with carbon dust are on average redder in K − IRAC2. Previous works have noticed that this color index becomes negative at low temperature for RSGs (see Yang et al. 2019), which is believed due to the CO absorption in the IRAC2 band. Yet an RSG with carbon dust has a positive value, indicating normal continuum emission of dust, which may imply that the CO band absorption can be neglected. In contrast, an RSG with carbon dust is bluer in K − IRAC4. This can be understood if the carbon-rich dust brings about some PAH emission in the IRAC4 band, though the dust species used in the model fitting is amorphous carbon. The discrimination of carbon and silicate    dust in these two color-color diagrams generally agrees with our model fitting of the dust species. Nevertheless, there does exist some confusion in both diagrams, and no clear borderline can be recognized. Definitely, the SED of some RSGs cannot be explained by silicate dust, and an additional dust component must be supplemented. At present, carbonaceous dust seems to the only applicable one. The percentage of carbon dust derived in this work still has some uncertainty, but this comes from the model fitting. In later analysis, we will divide the sample into carbon and silicate RSGs according to the model-fitting results.

Mass-loss Rate
The distribution of MLR is presented in Figure 11. The statistical properties are listed in Table 4, where all the values are calculated with a gas-to-dust ratio of 100 and the column "all" includes the optically thin sources. The whole catalog is listed in Tables 1 and 2. The MLR ranges from 10 −9 to 10 −3 M e yr −1 , with a median value of ∼3.1 × 10 −6 M e yr −1 , which is an order of magnitude smaller than the mean value of ∼2.0 × 10 −5 M e yr −1 , caused by numerous sources with very low MLR, i.e., the optically thin sources. In general, the average MLR in M31 is comparable to M33. On the other hand, the summed MLR of the sample RSGs in M33 (3.97 × 10 −2 M e yr −1 ) is about 13% higher than that in M31 (3.51 × 10 −2 M e yr −1 ). If the RSGs are divided into two groups according to the species of circumstellar dust, some difference appears. The MLR is higher in the case of silicate dust than for carbon dust on average, and the histogram shows clearly that the carbon dust sources lack very high MLR. In addition, the range of MLR in the case of carbon dust is more concentrated. The reason may come from the difference in the dust properties. As illustrated in Section 3.1, the parameter a Q e / is 1.573 for silicate and 0.263 for amorphous carbon, and the MLR is proportional to it. This implies that the MLR differs by a factor of 6 between silicate and  carbon circumstellar shells even if they have the same optical depth and inner dust radius. The density difference, i.e., that ρ silicate = 3.5 g cm −3 is 1.5 times ρ AC = 2.26 g cm −3 , yields even higher MLR in the case of silicate dust. Apparently the higher-MLR RSGs with silicate dust need not be redder than those with carbon dust, which agrees with the observations.
Our results can be compared with the MLR of RSGs in previous works mentioned in the Introduction. All the previous values of MLR lie apparently within the range of the present sample from 10 −9 to 10 −3 M e yr −1 . The MLR of a few bright RSGs in the Milky Way from Gordon et al. (2018)     the high end of our sample, which can be understood as due to either their high luminosity or large optical depth. The large sample of AGB and RSG stars in the LMC from Riebel et al.
(2012) exhibits a range of dust MLR between 10 −7 and 10 −11 M e yr −1 , which coincides with our results as well if the gas-todust ratio of 400 is adopted. Generally our results support previous works with a much more complete sample that covers various cases of MLR of RSGs.

Relation of MLR with Stellar Luminosity
The relation of MLR with stellar parameters is widely calibrated, and some of these are mentioned in the Introduction. With the measured J − K, the effective temperature T eff and the bolometric correction to the K magnitude are calculated in Section 2.4. A linear relation is fitted between MLR and L L log( ) /  . Excluding the optically thin RSGs, the relations derived for silicate and amorphous carbon dust are  Table 5 shows that the correlation is moderate. The results are displayed in Figure 12, and compared with previous results in Figure 13. The slope, i.e., the power-law index γ in the relation µ g M L  , is the same in various cases, i.e., the same in M31 and M33, and the same for silicate and carbon dust. This index also agrees with those from van Loon et al. (2005) andde Jager et al. (1988) for dusty RSGs in the LMC. It should be mentioned that van Loon et al. (2005) andde Jager et al. (1988) actually derive the relations of MLR with not only luminosity but also effective temperature, and the -M L  relation displayed in Figure 12 is calculated with a mean effective temperature of our sample of 3921 K. In addition, this index is around 1.0, which means that MLR is approximately proportional to luminosity. On the other hand, the index is apparently smaller than that from Beasor & Davies (2018) for RSGs in Galactic clusters. The discrepancy of Beasor & Davies (2018) with others may be caused by their sample being different, and that they adopted a gas-to-dust ratio of 200.
Unlike the common slope, the intercept, i.e., the coefficient of proportion in the µ g M L  relation, is not the same. As discussed in Section 3.1, the amount of MLR depends on multiple factors, including the dust species, the dust size, the stellar expansion velocity, and luminosity. The effects of dust species and size are discussed in previous sections, and can be as large as one order of magnitude. The expansion velocity and luminosity should have less effect, probably changing the MLR by a factor of a few. Taking all the factors into account, the MLR may differ by two orders of magnitude for a given observational SED, therefore the comparison of the vertical shift is meaningful only under the same model conditions.
Taking T eff into account, the relationships become   Figure 14 together with the value of R 2 , which implies an improved relation in comparison with the luminosity alone. In general, the fitting is reasonably good with the residual -< M M log log 1.0 predict | |   , but it deviates somewhat from the measurements, particularly for the silicate dust. This may be understood because MLR depends on other parameters in addition to L and T eff , such as stellar mass, radius or surface gravity (Reimers 1975;Beasor et al. 2020), and pulsation period (Goldman et al. 2017). But these parameters are not conveniently available from observation.

Relation of MLR with Infrared Brightness
After examining the infrared-band brightness, a very good relation is found with the mid-infrared flux in the Spitzer/ IRAC4 band at 8 μm and the MIPS24 band at 24 μm. Linear fitting yields the following relations with the absolute magnitude, for which the distance modulus is 24.32 for M31 and 24.50 for M33 (Paper I).


The fitting lines are shown in Figures 15 and 16. The correlation coefficients are all large (<−0.7) as listed in Table 5, which implies a tight relation of the mid-infrared flux with the MLR. This relation is expected from the warm temperature of the circumstellar dust, which emits mainly in the mid-infrared. On other hand, the relation with near-infrared flux is weak and not presented.

Relation of MLR with Infrared Color Indexes
Examining the infrared colors, it is found that the MLR is sensitive to IRAC1 − IRAC4: As seen in Figure 17, the scattering around the fitting line is significant, and the correlation of MLR with the color index is not so tight as with the mid-infrared flux. The correlation coefficients in Table 5 also manifest that the correlation is very weak for silicate dust and moderate for carbon dust. Previously, Groenewegen & Sloan (2018) and Matsuura et al. (2009) derived the relation of MLR with IRAC1 − IRAC4 as well for AGB stars, but they take an exponential function as illustrated in Figure 17 by the dasheddotted lines, which hardly delineates the tendency for the RSG sample. It seems no good relation exists for the MLR with the infrared colors.

Total MLR
The initial sample from Paper I should be complete since the lower limit of the RSG brightness in M31 and M33 is higher than the observationally complete magnitude. But the sample for which we derived the MLR is incomplete after removing those without appropriate photometric SED. Nevertheless, the comparison with the initial sample shown in Figure 3 finds no bias to the K magnitude or the J − K colors. Thus we simply multiply the summed MLR from our sample by the ratio of the numbers in the initial and final samples, i.e., 5253/1733 = 3.03 for M31 and 3001/1979 = 1.52 for M33. As listed in Table 4, the summed MLR of all the modeled sample RSGs in M31 is ∼0.035 M e yr −1 , and that in M33 is ∼0.040 M e yr −1 . After multiplying by this ratio, the total MLR of RSGs is then 0.11 M e yr −1 in M31 and 0.060 M e yr −1 in M33. Dividing by the gas-to-dust ratio, i.e., 100 in this work, the rate of dust production by RSGs is   10 −3 M e yr −1 in M31 and 6 × 10 −4 M e yr −1 in M33. At present, we have no idea of the dust mass contributed by AGB stars in M31 or M33. Taking the Milky Way galaxy as a reference, the rate of dust production by evolved low-mass stars is a few 10 −3 M e yr −1 according to Draine (2009), or about 10 −3 M e yr −1 according to Jones (2001). Though the uncertainty of these values may be a factor of a few, it seems that RSGs can be a non-negligible source of interstellar dust. A precise estimation of the proportion of RSG dust can be made only after a systematic calculation of the MLR of AGB stars in M31 and M33 in our further work.

Influence of Metallicity and Surface Gravity
Metallicity and surface gravity change the stellar spectrum as well as effective temperature, and they are assumed to be constant in modeling. The metallicity of both M31 (Blair et al. 1982) and M33 (U et al. 2009) is found to exhibit some gradient with the distance from the center. Considering that RSGs are located in the spiral arm or ring area, the case of [Z] = 0.25 with the MARCS RSGs model is examined. Although Goldman et al. (2017) derived a relation of the wind speed proportional to metallicity, a scale of [Z] = 0.25 would yield an unreasonably high expansion velocity. Moreover, they suggest that the mass-loss rate is independent of metallicity if the latter is between a half and twice solar. Thus the velocity is still calculated with [Z] = 0. The resultant MLR is listed in Table 6. It can been seen that a higher metallicity will bring about a higher MLR, which is possibly caused by a lower bump in the J, H, K bands.
Quantitatively, the mean MLR shows very little difference when [Z] changes from 0 to 0.25, with a value of 2.03 × 10 −5 M e yr −1 , meaning an increase of 1%. For the surface gravity, the case of =g log 0.5 is examined. As shown in Table 6, the mean MLR increases to 2.03 × 10 −5 M e yr −1 with =g log 0.5, i.e., by 1%.