Evolved Massive Stars at Low Metallicity. VI. Mass-loss Rate of Red Supergiant Stars in the Large Magellanic Cloud

Mass loss is a crucial process that affects the observational properties, evolution path, and fate of highly evolved stars. However, the mechanism of mass loss is still unclear, and the mass-loss rate ( MLR ) of red supergiant stars ( RSGs ) requires further research and precise evaluation. To address this, we utilized an updated and complete sample of RSGs in the Large Magellanic Cloud ( LMC ) and employed the 2-DUST radiation transfer model and spectral energy distribution ﬁ tting approach to determine the dust-production rates ( DPRs ) and dust properties of the RSGs. We have ﬁ tted 4714 selected RSGs with over 100,000 theoretical templates of evolved stars. Our results show that the DPR range of RSGs in the LMC is 10 − 11 M e yr − 1 – 10 − 7 M e yr − 1 , and the total DPR of all RSGs is 1.14 × 10 − 6 M e yr − 1 . We ﬁ nd that 63.3% RSGs are oxygen-rich, and they account for 97.2% of the total DPR. The optically thin RSG, which comprise 30.6% of our sample, contribute only 0.1% of the total DPR, while carbon-rich RSGs ( 6.1% ) produce 2.7% of the total DPR. Overall, 208 RSGs contributed 76.6% of the total DPR. We have established a new relationship between the MLR and luminosity of RSGs in the LMC, which exhibits a positive trend and a clear turning point at » L L log 4.4  .


INTRODUCTION
Stars with initial masses between 7 and 30 M ⊙ follow a path of evolution that leads them to become red supergiants (RSGs) before ultimately ending their lives as core collapse supernovae (Smartt et al.Corresponding author: Jian Gao, Ming Yang, Bingqiu Chen jiangao@bnu.edu.cn,myang@nao.cas.cn,bchen@ynu.edu.cn2009).During the RSG phase, these stars have low effective temperatures ranging from 3,500 to 4,500 K, high luminosities from 4,000 to 400,000 L ⊙ , and large radii from 100 to 1,500 R ⊙ (Massey & Johnson 1998;Massey & Olsen 2003;Levesque 2010Levesque , 2017)).For stars with initial masses below 30 M ⊙ , their mass loss mainly occurs during the RSG phase (Beasor et al. 2020).The type of supernovae that a massive star will become is primarily determined by its initial mass, metallicity, as well as its mass-loss rate (MLR) during the RSG phase (Humphreys 2010;Ekström et al. 2013;Smith 2014;Meynet et al. 2015;Beasor & Davies 2018).Moreover, RSGs are important contributors to the interstellar medium.The mass loss and evolution of RSGs will eventually influence the chemical evolution of their surrounding mediums.However, the mass loss mechanism of RSGs, which is the result of a combination of factors including radiation pressure, stellar pulsation, stellar wind, convection, as well as the luminosity and metallicity of the stars , is not yet fully understood (Gehrz & Woolf 1971;MacGregor & Stencel 1992;Harper et al. 2001;Wilson et al. 2000;Yoon & Cantiello 2010;Mauron & Josselin 2011;Beasor & Davies 2016).
Precisely determining the MLRs or dust-production rates (DPRs) of RSGs is necessary to investigate their mass-loss mechanism.One commonly used method to obtain the MLRs of RSGs is to fit the observational optical to mid-infrared (MIR) spectral energy distributions (SEDs) of the RSGs with theoretical models.Using the DUSTY radiative-transfer code (Ivezic & Elitzur 1997), Wang et al. (2021, hereafter W+21) found an average MLR of ∼ 2.0 ×10 −5 M ⊙ yr −1 for 1,741 and 1,983 RSGs in M31 and M33.Yang et al. (2023, hereafter Y+23)found that the typical MLR of RSGs in the Small Magellanic Cloud (SMC) is about 10 −6 M ⊙ yr −1 .Gordon et al. (2018) estimated the MLR of some well-known RSGs to be ranging from ×10 −4 M ⊙ yr −1 to ×10 −6 M ⊙ yr −1 .Riebel et al. (2012, hereafter R+12) obtained DPRs based on 30,000 RSGs and asymptotic giant branch stars (AGBs) in the Large Magellanic Cloud (LMC) using the GRAMS model grid (Srinivasan et al. 2011;Sargent et al. 2011), which were mostly between 10 −7 and 10 −11 M ⊙ yr −1 .Boyer et al. (2012, hereafter B+12) calculated the DPRs of AGBs and RSGs in the LMC and SMC from their infrared excesses and obtained a total DPR of approximately 2.4 × 10 −7 M ⊙ yr −1 for 3,908 RSGs in the LMC.Srinivasan et al. (2016, hereafter S+16) conducted a study on the DPR of 1,410 RSGs in the SMC, which had a total DPR of 4.6 × 10 −8 M ⊙ yr −1 and an average DPR of 3.3 × 10 −11 M ⊙ yr −1 .The MLRs obtained from these works vary significantly, mainly due to the use of different dust models and possible bias of the adopted samples sizes (Srinivasan et al. 2016;Groenewegen & Sloan 2018).A complete and pure RSG sample is vital for us to accurately calculate the MLRs of RSGs.Ren et al. (2021) constructed comprehensive and uncontaminated RSG samples of Local Group galaxies including the LMC and SMC, which offers the oppotunity to investigate RSGs mass loss in the LMC and their overall contribution to interstellar medium.The following sections of this paper are organized as follows: Section 2 details the sample, Section 3 presents the dust model grid, Section 4 explains the fitting procedure, Section 5 shows the results of this research, and Section 6 summarizes the paper.

DATA
In the current work, we adopt the RSG sample from Ren et al. (2021), who utilized the nearinfrared photometric data from UKIRT/WFCAM (Irwin 2013) and the 2MASS point source catalog (Skrutskie et al. 2006) and astrometric data from Gaia EDR3 (Gaia Collaboration et al. 2021) to identify RSGs in 12 low-mass galaxies of the Local Group.In their sample, there are 4,823 RSGs in the LMC.
In order to estimate MLRs of the RSGs, SEDs with wide coverage of wavelength range is necessary.The infrared (IR) photometry is essential as it reveals information about the circumstellar dust, while the optical data provide atmospheric parameters of the star.We crossmatched the RSG sample with various photometric data, including Skymapper (u, v, g, i, r, and z bands; Keller et al. 2007;Bessell et al. 2011;Wolf et al. 2018), Gaia EDR3 (G BP , G, and G RP bands; Gaia Collaboration et al. 2021), Magellanic Clouds Photometric Survey (MCPS; U, B, V, and I bands; Zaritsky et al. 2004), Vista Magellanic Cloud Survey (VMC; Y, J, and K S bands; Cioni et al. 2011), InfraRed Survey Facility (IRSF; J, H, and K S bands; Kato et al. 2007), The Two Micron All Sky Survey (2MASS; J, H, and K S band; Skrutskie et al. 2006), AKARI (N3, S7, S11, L15, and L24 bands; Onaka et al. 2007;Murakami et al. 2007;Kato et al. 2012), Wide field Infrared Survey Explorer (WISE; W1, W2, W3, and W4 bands; Wright et al. 2010), and Spitzer ([3.6], [4.5], [5.8], [8.0], and [24] bands; Werner et al. 2004), as shown in Table 1.A search radius of 1 ′′ was used in the crossmatching.These photometric data provided a broad wavelength coverage ranging from ultraviolet (0.36 µm) to IR (24 µm).To ensure data quality, we also imposed additional constraints on the photometric data for each band as following: 1.For Skymapper, Gaia, VMC, 2MASS, IRSF, and AKARI data, we required that the magnitude errors must be less than 0.1 mag.
3. For Spitzer data, we require the SNR of each band to be no less than 3 and the close source flag close f lag = 0.
Meanwhile, it was not necessarily for each target in our sample to have observations in all filters and satisfy above conditions, which would significantly reduced the number of stars in our sample.However, to ensure relatively accurate estimation of MLRs, 2MASS JHK S measurements for each star must be good.We also required that each source must had good observations in WISE W 3 or Spitzer [8.0] band, which was important to constrain the dust emission information.After applying these criteria, 85 sources were excluded based on photometric quality constraints.Among these 85 sources, one did not meet the 2MASS criteria , while the remaining 84 lacked data meeting the selection criteria for W3 or [8.0].We performed cross-matching with Simbad (Wenger et al. 2000) and removed 24 sources that were identified as planetary nebula, post-AGB, or YSO etc.After applying these criteria, there were 4,714 RSGs in the final RSG sample.
We then corrected for the extinction effect of our sample.In the current work, we utilized the extinction map of LMC by Chen et al. (2022) and the extinction coefficients from the classical CCM89 model (Cardelli et al. 1989) with R V = 3.41 (Gordon et al. 2003).In our sample, there are 8 stars are not covered by the extinction map of Chen et al. (2022) and another 155 stars have reddening values E(B − V ) < 0 in the map, which could be due to minimal extinction and errors.We used the same method as Y+23 to correct the extinction of these 163 sources, setting R V = 3.1, E(B−V ) = 0.1 (Schlafly & Finkbeiner 2011;Freedman et al. 2020), adopting the extinction coefficient obtained by Wang & Chen (2019).The 2MASS J and K S bands color-magnitude diagram (CMD) of the background targets and final RSG sample after extinction correction are presented in Figure 1.
Various radiation transfer models are available for calculating the MLR, including one-dimensional model DUSTY (Ivezic & Elitzur 1997), two-dimensional model 2-DUST (Ueta & Meixner 2003), and three-dimensional model MCMax3D (Min et al. 2009), etc. Y+23 used the 1D spherically symmetric DUSTY radiative transfer model and identified a prominent turning point on the luminosity-MLR diagram for RSGs.We speculate that this turning point may be related to the mass loss mechanism of RSGs (e.g., Vink & Sabhahit 2023).To confirm the presence of this turning point and determine if it is influenced by factors such as models and samples, we decided to employ different radiative transfer models.We adopted 2-DUST to generate the theoretical SEDs.2-DUST is a flexible, axisymmetric 2D radiation transfer code with highly adaptable input parameters.We assumed a spherically symmetric circumstellar dust shell for our study.With input parameters such as the spectrum of the central star, dust shell geometry, and dust grain properties, 2-DUST solves the radiation transfer problem and provides DPR, luminosity, and other output parameters of the star (Ueta & Meixner 2003).The settings of our input parameters are outlined in Table 2.We also provided the parameter settings of GRAMS (Sargent et al. 2011;Srinivasan et al. 2011), a template library generated using the 2-DUST radiation transfer model, for comparison.The input parameters we used in 2-DUST can be divided into three parts: the spectrum of the central star, the shape of the dust shell, and the dust grain properties of the circumstellar dust.We described them in detail below.

Photosphere Model of the centeral star
For the spectrum of the central star, we adopted the MARCS stellar atmosphere model (Gustafsson et al. 2008).A blackbody spectrum may be easier to use, but it fails to reflect evolutionary effects such as the "H-bump" feature (Yang et al. 2021).Our aim was to generate an SED grid that covers the properties of the evolved stars, including both AGB stars and RSGs in the Magellanic Clouds.40 MARCS models are chosen as templates for RSGs, with temperatures ranging from 3,300 to 4,500 K (3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000, 4250, 4500 K), a mass of 15 M ⊙ , surface gravity values ranging from −0.5 to 0 dex, and metallicities ranging from −0.5 to 0 dex.For evolved stars in other mass ranges, we also used 78 MARCS models with effective temperatures ranging from 2,500 to 5,000 K (2500, 2600, 2700, 2800, 2900, 3200, 3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000, 4250, 4500, 4750, 5000 K), masses ranging from 0.5 to 5 M ⊙ (0.5, 1, 2, 5 M ⊙ ), surface gravity values ranging from −0.5 to 0 dex, and a metallicity is 0.5 dex.Among them, the MARCS models with a temperature higher than 4500 K are only used to generate oxygen rich (O-rich) templates.We extended the wavelength coverage of the MARCS spectra to 1,000 µm by adding blackbody in the infrared region ( ≥ 20µm).The MARCS models were then downgraded to a resolution of R = 100 and fed into 2-DUST.To reduce the computation time, we set the number of data points output by the 2-DUST to 142. Figure 2 shows the 2-DUST input and output spectra of an example model.5,7,10,12,15,20,25 3,5,7,10,12,15,20,25 1.5,3,4.5,7,12 3,7,11,15 Rout (R in) 1000 1000 1000 1000 vexp (cm s −1 ) 10 10 10 10 We assumed that the density distribution of the dust envelope followed a power-law relation of ρ(r) ∝ r −2 , and that the dust envelope was spherically symmetric.For templates with masses between 0.5 and 5 M ⊙ , we calculated the model with inner dust shell radii of R in = 3, 5, 7, 10, 12, 15 R ⋆ , and for those with mass of 15 M ⊙ , R in =5, 7, 10, 15, 20, 25 R ⋆ , for which R ⋆ was the radius of the central star.The ratio between the outer radius and inner radius was assumed to be R out /R in = 1000.

Circumstellar Dust Shell
The temperature of the inner dust shell was determined by its size in 2-DUST.Previous works estimated the condensation temperature of silicate dust to be around 400 to 1,500 K, which was influenced by the central star and the dust composition (Gail & Sedlmayr 1984, 1999;van Loon et al. 2005;Sargent et al. 2010;Gail et al. 2020).On the other hand, the condensation temperature of carbon dust is higher (Groenewegen et al. 1994;Nanni et al. 2013;Brunner et al. 2018), and the upper limit of carbon dust temperature in GRAMS is set to 1800 K (Srinivasan et al. 2011).For our case, we used 1,800 K as the upper limit for carbon-rich (C-rich) models and 1,400 K for O-rich models.
The calculation of DPR by 2-DUST requires wind speed, which is difficult to measure.Previous studies of evolved AGB and RSG stars obtained a wide velocity range (van Loon et al. 1999;Ramstedt et al. 2006;De Beck et al. 2010;Cox et al. 2012;Marshall et al. 2004).We assumed a constant expansion velocity (v exp ) of v exp = 10 km s −1 without considering the influence of luminosity and metallicity, etc. v exp and DPR are simply proportional, that different DPR can be obtained by changing the values of wind speeds after all grids are generated.We will discuss later (Section 5.2) the result of scaling v exp and v exp = 25 km s −1 .
Based on above configuration, we generated a total of more than 10,000 models.We selected models with a dust shell inner temperature lower than 1,400 K (O-rich) and 1,800 K (C-rich), which yielded 77,458 models, consisting of 35,159 C-rich models and 42,299 O-rich models.To calculate the luminosity of each target, we assumed a distance of 50 kpc for the LMC (Westerlund 1997;Lee et al. 2011).The luminosity was estimated in three ways: using the 2-DUST output, the integrated observational SED (L obs ), or the integrated best fitted model SED (L mod ).The outputs for all the model are listed in Table 3, this table is available in the machine-readable format, a portion is shown here for guidance regarding its form and content.f The mass of the entire dust shell. g The mass of the central star.

SED FITTING
We determined the best-fit model by minimizing the χ 2 value, χ 2 for each source was defined by, where i is the number of the model, F (M, λ) are the simulated photometric data derived by the convolution of the model flux and each individual filters5 , F (O, λ) are the observed photometric data, N |f (O, λ)| are the number of adopted observation bands, and different for each source.C are the weight used for data of different wavelengths.In order to increase the weight of the infrared data, the C is changing based on different wavelengths, e.g., C = 5 for λ ≥ 1.0 µm and C = 1 for λ < 1.0 µm.The distribution of the number of measurements, N |f (O, λ)|, ranging from 6 to 32 filters for our RSG sample, was shown in Figure 3.
Prior to fitting, all the sample RSGs are split into two groups.The first group, named "S1", is consists of RSGs lacking W3 data or photometric data with wavelength longer than the W3 band.The second group, named "S2", is comprised of RSGs containing measurements of wavelength no less than the W3 band (i.e. the WISE W 3 and W 4, Spitzer [24], AKARI S11, L15 or AKARI L24 band).Stars in the S1 sample can only be fitted with models where τ 9.7/11.3≤ 0.01, while those in the S2 sample can be fitted with all models.Photometric data in those long wavelength are essential in determining the dust species present in a circumstellar dust shell (specifically, Si-O stretching ∼ 9.7 µm and O-Si-O bending ∼ 18 µm).The number of sources for S1 and S2 is 1,288 and 3,426, respectively.Splitting the RSG sample into two subgroups is necessary, since fitting all sources to all models without restrictions may lead to stars without observations in longer wavelengths being fitted to models with strong dust emission, even though these sources are supposedly without dust or only with very small amount of dust.If there is no restriction for S1, their DPR may be excessively overestimated, as shown in Figure 4 (see also details below).Moreover, as a precaution, we also visually inspected all the SEDs and tentatively pre-fitted all models for all targets and discovered that, 26 RSGs in the S1 sample showed a significant rise in [3.6], [4.5], [5.8], and [8.0] bands (based on the visual inspection), these sources were then moved to the S2 sample.
Based on the best-fit results, we categorized all RSGs into three distinct types.The first type, named "carbon-rich"(C-rich), is characterized by a significant continuous rise in the MIR region of the SED, without obvious emission peak.The second type, labeled as "oxygen-rich" (O-rich), displays steep infrared emission peaks at 9.7 µm and 18 µm.If a RSG showed no observable dust emission (τ 9.7/11.3= 0.0001, which was the lower limit of optical depth in our model), regardless of whether it was fitted to the C-rich model or the O-rich model, we considered it as "optically thin" (opt-thin).These optically thin RSG stars have SEDs that are similar to the central star, indicating that they have negligible amount of circumstellar dust and mass loss rate.In such cases, there was no difference between the C-rich model and the O-rich model, so we used the O-rich model for all optically thin RSGs.During the fitting process, it was difficult to distinguish the chemical type of dust when the optical depth is small (e.g., τ 9.7/11.3≤ 0.005).Therefore, we classified all sources with an optical depth less than 0.005 as O-rich, since most of the RSGs are expected to be O-rich.Such  weak dust features, regardless of their type, are unlikely to have a significant impact on the calculation of the total DPR.Finally, we visually checked the fitting results of all RSGs and manually modified the fitting results of 141 sources, for which most of these sources had ambiguous dust properties that made it difficult to distinguish their types (we made the judgments based on our experience).

RESULT AND DISCUSSION
Table 4 shows an example of the catalogue with photometric data, chemical composition, optical depth, luminosity, T eff , derived DPR, χ 2 etc. Figure 5 shows several examples of the fittings of optically thin, C-rich, and O-rich RSGs.The model SEDs generally match well with the measurements, with a few ones show poor fittings.The bottom row in Figure 5 shows the SED of two examples of poorly fitted RSGs, which could be due to bad W 3/W 4 observations (left) or binarity (right).We defined a "sf lag," and for those normal RSGs, sf lag = 0.There are 55 sources identified as binaries due to their blue end flux excess, among them, the SEDs of 23 sources exhibit significant anomalies, and we set their sf lag to 2, for the remaining 32 sources whose fitting results were not affected, their sf lag was set to 2.5.Targets similar to the left panel in the bottom row, which have bad W 3 and W 4 observation, their sf lag was set to 3, we identified 245 such sources, which also exhibited anomalies in the luminosity-MLR diagram (see below).The complete table contains initial photometric data and error in all bands, covering a range of 0.3 to 24 µm, and the complete table contains flux after extinction correction in all bands too. c Our defined classification flag, the sources with sf lag = 0 are normal sources, sf lag = 2 represents binaries or contamination, while sf lag = 2.5 indicates only slight flux excess that does not affect the fitting results.sf lag = 3 represents sources with bad W3/W4 band observational data. d The complete table includes more stellar and dust parameters, including total mass of the dust shell, inner radius temperature, mass loss timescale, etc. e The best-fitting model for this RSG corresponds to a specific index number in the model library (referred to as "index" in Table 3).

Different types of the RSG dust shells
The classification statistics of all the sample RSGs are presented in Table 5.As a result, 1,444 RSGs are classified as optically thin, while over half are O-rich with 2,984 RSGs falling under this category.Only 286 RSGs are classified as C-rich.Figure 6 shows the χ 2 distribution.The optically thin group appears to be better than O-rich and C-rich types, as indicated by the smaller χ 2 values.We noted that in our study, we adopted a strict definition of "optically thin", where the optical depth of the model SED was τ 9.7/11.3≤ 0.001, and the infrared dust emission was minimal.If we relax this definition to include an optical depth τ 9.7/11.3≤ 0.001, the proportion of optically thin RSGs would increase to 35.41%.To verify the accuracy of our classification, we crossmatched our catalogue with the spectroscopic data from the Spitzer IRS (InfraRed Spectrograph) Enhanced Products (Houck et al. 2004) using a search radius of 1 ′′ .We identified 50 sources that had matches.Visual inspection of their spectra revealed that 43 of the 50 sources were correctly classified, resulting in an accuracy of 86%.The 7 misclassified sources were all O-rich stars mistakenly identified as C-rich.This is possibly due to that there are both silicate emission peaks and prominent continuous emission in the NIR to MIR bands in their SEDs.We corrected the classification of these sources in our resultant catalogue based on their spectra.Figure 7 shows the examples of IRS spectra with high DPR and their classification in our sample.
The classification of lower mass stars such as AGB stars as C-rich can be attributed to carbon being brought to the surface of the stars through deep convection and dredge-up processes (Herwig 2005;Marigo et al. 2008;Pastorelli et al. 2019).However, the formation mechanism of carbon-rich RSGs, which lacks the mechanism to bring carbon to the surface due to their higher mass, requires further investigation.
The distributions of the resultant values of optical depth, effective temperature, and luminosity of the RSGs with different types are shown in Figure 8.Our results indicate that almost all RSGs (96.80%) have temperatures between 3,300 K and 4,500 K, and luminosities between 3,000 and 400,000 L ⊙ .This is consistent with previous findings (Massey & Johnson 1998;Massey & Olsen 2003;Levesque 2010Levesque , 2017)), which suggests the purity of our sample.No significant differences are visible in the distributions of these parameters among the different types, except for that the C-rich RSGs tend to have smaller optical depths.

The resultant DPR and MLR
We assumed a gas-to-dust ratio of ϕ = 500 (Riebel et al. 2012) for the LMC.The MLR value of each RSG was calculated based on its best-fit DPR value and the assumed gas-to-dust ratio.The  histograms of derived DPR and MLR for the RSG sample are presented in Figure 9. Optically thin RSGs have extremely low MLR.The C-rich RSGs exhibit a higher MLR than the O-rich ones.This could be partly due to our classification method which only recognizes RSGs with obvious carbon emission.
The typical MLR value of our sample RSGs is 10 −8 to 10 −7 M ⊙ yr −1 , while a few can reach to 10 −5 M ⊙ yr −1 .Table 5 provides detailed DPR values for different types of RSGs.The total DPR including is 1.14×10 −6 M ⊙ yr −1 , with O-rich RSGs contributing the most (97.24%),followed by C-rich RSGs and optically thin RSGs.The average DPR of our sample RSGs is about 2.41 × 10 −10 M ⊙ yr −1 , with a median value of about 2.48 × 10 −11 M ⊙ yr −1 .208 RSGs with a DPR higher than 10 −9 M ⊙ yr −1 contributed 76.57% of the total DPR.The fitting SEDs for several RSGs with high DPR are shown in Figure 7.The values of MLR and DPR of RSGs in the LMC were also estimated in previous work.R+12 obtained the global DPR of the LMC from both RSGs and AGB stars, which was around 2.1 × 10 −5 M ⊙ yr −1 .Among them, RSGs accounted for 9.4%, and the total DPR of all RSGs was about 2.0 × 10 −6 M ⊙ yr −1 .B+12 found that the total DPR of 3,908 RSGs in the LMC was about 2.4 × 10 −7 M ⊙ yr −1 with an average DPR of about 6.1 × 10 −11 M ⊙ yr −1 .S+16 estimated that the total DPR of 1,410 SMC RSGs in their sample was 4.6 × 10 −8 M ⊙ yr −1 , with an average DPR of 3.3 × 10 −11 M ⊙ yr −1 .Our results are basically agreement with B+12, R+12 and S+16 within the uncertainties.Table 6 shows our results and the comparison with previous works.Compared to the result of R+12 (all evolved star sammple, including AGB and RSGs), our DPR (only RSGs with sf lag = 0 are included) is slightly lower as shown in Figure 10.But our resultant luminosity values are consistent with theirs.To study the impact of sample size differences on the results, we cross-matched our sample with their RSG sample and identified 4,213 common sources.For these sources, the cumulative R+12 DPR was 1.45 × 10 −6 M ⊙ yr −1 , which is a 27.5% reduction compared to the result of their original result.Sample size undoubtedly directly affects the results.For the 4,213 RSGs, our calculated cumulative DPR is 9.0 × 10 −7 M ⊙ yr −1 , which is still slightly lower than that of R+12.This overall lower trend may be attributed to differences in model settings, fitting strategies, or other factors (see the detailed disscusion in Section 5.1 of Y+23).
Our total DPR is 4.6 times that of B+12, we cautiously regard this difference as within an acceptable range, particularly considering the difference in our research approaches.B+12 employed a hybrid approach in which they initially fitted the DPR for several evolved stars, then established the relationship between DPR and [8.0] excess, finally applied this relationship to the entire sample.The sample they utilized was from Boyer et al. (2011), for which the sample was selected based on 2MASS (Skrutskie et al. 2006) and SAGE (Wright et al. 2010), the spatial coverage of the sample is smaller than that of Ren et al. (2021).We identified 3,855 sources in common.For those sources, our calculated total DPR is 6.74 × 10 −7 M ⊙ yr −1 .
The average MLR of RSGs in M31 and M33 was found to be about 2.0 × 10 −5 M ⊙ yr −1 by W+21.Y+23 found that the typical MLR of RSGs in the SMC is about 10 −6 M ⊙ yr −1 .These values are larger  than our results.The determination of MLR is particularly sensitive to the choice of parameters such as the optical constants, gas-to-dust ratio, and dust shell expansion speed, etc., which can cause the difference in results to be up to one order of magnitude (Srinivasan et al. 2016;Yang et al. 2023).Groenewegen & Sloan (2018, hereafter GS18) fitted the SEDs and IRS spectra of nearly 400 evolved stars in the LMC, SMC, and other dwarf spheroidal galaxies in the Local Group.They estimated the MLR for each target ranging from 10 −4 to 10 −11 M ⊙ yr −1 , which also highly depended on the adopted optical constants.Moreover, we varied some parameters and methods to assess their impact on our DPR and MLR results.The statistical results are shown in Table 7.When we used the optical constant of Dorschner et al. (1995, hereafter Do+95) to generate the model for fitting, the resulted DPR was 22% lower and the chemical classification were also changed.The number of optically thin sources increased to 1,547, while the number of C-rich sources increased to 300, and the number of O-rich sources decreased to 2,867.Their contributions to total DPR were also changed, for which the detailed statistical results are listed in Table 8.
In summary, the impact factors of DPR obtained through the SED fitting method are complex.We can only provide a rough estimation, and the uncertainty of the DPR may reach up to one order of magnitude or even more.

The MLR and stellar parameters relations
The relations between MLR and stellar parameters, such as the luminosity, temperature, radius, and initial mass, etc, are extensively calculated in the literature.In this work, we utilized a third-order polynomial and a piecewise function to fit the data and obtained the expressions for luminosity (both L obs and L mod ) and the MLR, respectively.Figure 11 depicts the relationship between luminosity and MLR.To ensure the accuracy of the fitting, optically thin targets, or targets have a high level of uncertainty (τ 9.7/11.3≤ 0.001 and sf lag = 2, 2.5, 3) are excluded.These sf lag = 3 targets (245 stars) are mostly located in the region where log (L obs /L ⊙ ) < 4.5 and log Ṁ ≥ −7.5, characterized by low luminosity and high MLR, they are considered to be pseudo sources or affected by other nearby sources.The SED of an example star is shown in the bottom row left panel of Figure 5.These stars are excluded from the MLR-luminosity functions6 .As a result, we have third-order polynomial: We compared our results with those from the literature, including de Jager et al. (1988, hereafter D+88), Beasor et al. (2020, hereafter B+20), W+21 and Y+23 as shown in Figure 12.The Figure shows that our resulted MLR values are slightly lower than those of D+88, W+21, and Y+23, particularly in the case of low luminosity.D+88 formulated a correlation from a relatively small sample, while both W+21 and Y+23 used the DUSTY radiation transfer model instead of 2-DUST with various parameter settings and different selections of optical constants, which might significantly affect the results.Our relations appear to be consistent with the results of B+20, especially for regions with log (L obs /L ⊙ ) ≥ 4.0.Our sample also shows similar turning point at log (L obs /L ⊙ ) ≈ 4.4 on the luminosity versus MLR diagram as reported by Y+23, which they call a "knee-like" shape, at log (L obs /L ⊙ ) ≈ 4.6 in the SMC.The similar trend with different turning point of luminosity-MLR relation may be due to the difference in metallicity between the LMC and SMC.Recently, Beasor & Smith (2022) re-evaluated the MLR of Dust-enshrouded RSGs and found that the relations of van Loon et al. (2005) likely overestimated the MLR of RSGs.B+20 also showed that the RSGs with smaller MLRs might not be able to lose all of their dust shells, returned to the blue end in the H-R diagram and became hydrogen-poor supernova.Similarly, our results indicate that the RSGs can not lose several solar masses with static stellar wind.Unlike AGB stars, RSGs may not be the main contributors of dust in the interstellar medium.
We also investigated the relationship between DPR and stellar colors of (J − [8.0]) 0 , (K S − W 3) 0 , ([3.8] − [24]) 0 , as shown in Figure 13 (only RSGs with sf lag = 0 are included).Both the Orich and C-rich RSGs exhibit positive correlations between the DPR and stellar infrared colors.Figure 14 shows our resulted ([3.6] − [8.0]) 0 color versus DPR diagram and the DPR-color relations from previous works.For comparison, we adopted a gas-to-dust ratio of 200 for both our results and the relations from the literature, in order to keep a consistent value of gas-to-dust ratio.The literature relations are mainly based on AGB and carbon stars.The RSGs show a roughly similar DPR-color relation.They are most similar to the M stars from GS18, but not so close to the carbon stars from R+12, Matsuura et al. (2009, M+09), and GS18.

SUMMARY
We establish a pre-computed library of evolved star templates using 2-DUST code (Ueta & Meixner 2003).Through SED fitting, the DPR and dust chemical properties of a complete and uncontaminated RSG sample of 4,714 RSGs in the LMC are determined.Our analysis reveal that the cumulative DPR of RSGs in the LMC is 1.14 × 10 −6 M ⊙ yr −1 .We find that the DPR heavily depends on a few  RSGs, as only 208 RSGs contribute 76.57% of the DPR.Among the sample, 1,444 RSGs are classified as optically thin, meaning they have a minimal dust envelope.Despite accounting for 30.63% of the total number, they only contributed 0.07% of the DPR.Conversely, O-rich RSGs, accounting for 63.30% of the total number, contributed 97.24% of the DPR, making them the primary drivers of mass loss among RSGs.A small proportion of C-rich RSGs (6.07%) contributed 2.69% of total DPR.We verify our classifications with 86% accuracy based on IRS spectrum.We fit the luminosity-MLR relation and confirm the turning point in the relation at log (L obs /L ⊙ ) ≈ 4.4.Our result is consistent with that obtained by Y+23.Given the low mass loss rate derived in this work, whether RSGs can completely lose all their dust envelope is a question that requires further investigation.
Similar to AGB stars, we found a positive correlation between DPR of RSGs and infrared colors.Future work on AGB stars and other galaxies will determine the exact proportion of the contribution of RSGs and AGB to the dust of the whole galaxy and its own evolution.
We are grateful to Hai-Bo Yuan and Jia-Ming Liu for very helpful discussions.This work is supported by the National Key R&D Program of China No. 2019YFA0405500 andNo. 2019YFA0405501, National Natural Science Foundation of China No. 12133002 andU2031209, 12203025, 12173034, 12373048 and11833006, and Yunnan University grant No. C619300A034, Shandong Provincial Natural Science Foundation through project ZR2022QA064.We acknowledge the science research grants from the China Manned Space Project with No. CMS-CSST2021-A09, CMS-CSST-2021-A08 and CMS-CSST-2021-B03.The numerical computations were conducted on the Yunnan University Astronomy Supercomputer.This research made use of the cross-match service provided by CDS, TOP-CAT (Taylor 2005).

Note- a Figure 2 .
Figure 2. Spectra of an example model with T eff = 4, 000 K, log g = 0.00 and M = 5 M ⊙ , [M/H] = −0.50,τ 9.7 = 0.2.The gray spectrum is the MARCS stellar spectrum.The blue and red spectrum are the input and output spectra of 2-DUST, respectively.The black vertical lines above the spectra represent the output wavelength grid.
12 ...... 2.44E... ...... Note-This table is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.aTheSED of the model covers a range of 0.2 to 1000 µm, with 142 grid points, which include dust emission and center star.b The synthetic flux calculated by the convolution of the 2-DUST model flux and the individual filters, we adopt the filter parameters from the Spanish Virtual Observatory (SVO) Filter Profile Service, see http://svo2.cab.inta-csic.es/theory/fps/.c The inner radius of dust shell, measured in units of the radius of the central star.d Temperature at the inner radius of the dust shell.e Density at the inner radius of the dust shell.

Figure 3 .
Figure 3. Histogram of numbers of data points for the RSG sample.

Figure 4 .
Figure 4.The result of direct fitting without restriction at long wavelengths (blue dot-dashed line) and the result of refitting after the restriction (black solid line) for an example RSG in the S1 sample.The observations are marked with different symbols of different colors.The black solid curves are the best fitted models.The black dashed and dot-dashed lines are the spectra contributed by the center star and the circumstellar dust, respectively.The blue dashed line represents the minimum χ 2 model of another chemical component.
table is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.aTheE(B − V ) value comes from the extinction map ofChen et al. (2022).b

Figure 5 .
Figure 5. SED fits of eight example RSGs, the first, second and third row show the good SED fits for optically thin, C-rich and O-rich stars, respectively (sf lag = 0).The bottom row shows anomalous SED fits of two example stars, due to bad W 3/W 4 observation (left, sf lag = 3), and binarity (right, sf lag = 2)).

Figure 6 .
Figure 6.Histograms of the resultant χ 2 .The red, blue and gray bars represent those of the C-rich, O-rich and optically thin stars, respectively.The vertical lines show the 95th percentile χ 2 values of the corresponding samples.

Figure 7 .
Figure 7. Examples of IRS spectra with high DPR and their classification in our sample.The red line represents the IRS spectral data, and other marks are the same as Figure 5.

Figure 8 .
Figure 8. Histograms of the best-fit values of optical depth (left), effective temperature (middle) and luminosity (right).Red, blue, and gray bars represent those of the C-rich, O-rich and optically thin RSGs, respectively.

Figure 9 .
Figure 9. Histograms of the best-fit values of DPR ( Ṁd , left) and MLR ( Ṁ , right).Red, blue, and gray bars represent those of the C-rich, O-rich and optically thin RSGs, respectively.

Figure 10 .
Figure 10.Comparison between our DPR and luminosity and the results of R+12 (left, middle), and the panel on the right is the comparison between the luminosity obtained by integrating the observed SED and the luminosity obtained by integrating the best fit model SED.

Figure 11 .
Figure 11.Luminosity-MLR diagram.The left panel shows the luminosity obtained by integrating the observed SED (L obs ), and the right panel shows the luminosity obtained by integrating the best-fit model L mod .The excluded and retained stars are marked by gray and black, respectively.The red curves represent the best-fitted relations, which are labeled in the diagrams.

Figure 12 .
Figure 12.Comparisons of our luminosity-MLR relations with those from previous works.Left panel show the luminosity obtained by integrating the observation SED (L obs ), and right panel shows the luminosity obtained by integrating the best fit model L mod .The retained stars are marked by black.Different curves are the relations from our work and the literature, which are labelled in the diagrams (see Section 5.3).

Figure 13 .
Figure 13.Stellar colors versus DPR diagrams for the optically thin (gray dots), O-rich (blue dots) and C-rich (red dots) RSG sample.

Table 2 .
Parameters for models

Table 3 .
Template library of evolved stars

Table 4 .
Catalog of the stellar and dust parameters for RSGs in LMC

Table 5 .
Dust-production rate ( Ṁd ) of RSGs in the LMC

Table 6 .
Comparison of total DPRs between different works

Table 7 .
Effect of different parameter settings on DPR