Critical Metallicity of Cool Supergiant Formation. I. Effects on Stellar-mass Loss and Feedback

We perform massive-star evolution simulations using the Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019) version No. 10108. We primarily use MESAstar, a 1D stellar evolution code that can model stellar evolution from the pre-main-sequence stage to the zero-age main-sequence (ZAMS) stage, and from the ZAMS to the iron core-collapse stage. It hydrodynamically evolves stars using microphysics such as nuclear reaction networks, equations of state, and opacities, which are the key physics of stellar evolution and are available in MESA.

MESA provides several stellar wind functions that can be used to calculate the mass-loss rates using stellar parameters at each time step during the stellar evolution. When a mass-loss routine is turned on, the code adjusts the stellar structure by rescaling the mass coordinates of the mass cells, according to the mass-loss rate given by the wind functions (Paxton et al.2011).

Our simulations include two steps. First, we evolve a pre-main-sequence star to reach the ZAMS stage, and then, we use the ZAMS model as the initial condition in the stellar evolution and evolve the model from the main sequence to the end of the star's life. Different nuclear reaction networks are employed in these two steps. For the pre-main-sequence evolution, we simply adopt the "basic" network with eight isotopes. For the evolution from the ZAMS stage onward, we employ the more complicated nuclear reaction network "approx21 (Cr60)," which includes 21 isotopes ⁴

In the stellar evolution simulations, we mainly use the default parameters for massive stars in the "inlist_massive_defaults" file in the MESAstar module. Additionally, the Ledoux criterion and Henyey's mixing-length theory (Henyey et al. 1965) are applied, and the mixing-length scaling parameter α (i.e., the ratio of mixing length to pressure scale height) is set at 1.5. The Type 2 OPAL opacity tables (Iglesias & Rogers 1996) are applied during and after the helium burning. In addition to the default parameters, we limit the change of helium abundance per time step to ≤10% to prevent abrupt changes.

Although the mass-loss rates can be calculated using theoretical approaches (e.g., Vink et al. 2000, 2001), such calculations are time-consuming and only available for hot winds. Instead, simple formulae of the mass-loss rates derived from observations or theoretical calculations are usually adopted in stellar evolution. Most of these widely used mass-loss prescriptions have been built in as modules in MESA. In our stellar evolution simulations, we apply the hot-wind prescription of Vink et al. (2000, 2001, hereinafter V01), the cool-wind prescription of de Jager et al. (1988, hereinafter dJ88), and the WR wind prescription of Nugis & Lamers (2000, hereinafter NL00). These prescriptions are described below.

For hot OB stars, metal ions in the stellar atmospheres absorb photons in the ultraviolet (UV) resonance lines and are radiatively accelerated, yielding steady-state winds (Lucy & Solomon 1970; Castor et al. 1975; Puls et al. 2008). Using multiline radiative transfers calculated with the Monte Carlo methods (Abbott & Lucy 1985), Vink et al. (2000, 2001) presented a grid of wind models, fitted the models with a linear regression, and established a mass-loss prescription for hot stars with effective temperatures T_eff > 10⁴ K.

When a massive star evolves into an RSG, its cool wind is driven by the radiation pressure on dust grains, similar to the winds from the asymptotic-giant-branch stars. To date, although empirical laws have been established, no well-established theoretical predictions for cool winds exist, because the origin and structure of RSG winds are very complex and the cool-wind mechanism is poorly understood (Smith 2014). Among the mass-loss rates suggested by the various empirical laws, up to a factor of 10 scatter has been seen (dJ88; Reimers 1975; Nieuwenhuijzen & de Jager 1990; van Loon et al. 2005; Mauron & Josselin 2011; Goldman et al. 2017; Beasor & Davies 2018; Beasor et al. 2020; Humphreys et al. 2020). In our simulations, we adopt the most commonly used cool-wind prescription by dJ88, which is an empirical function based on the fitting of the observed mass-loss rates of a set of OM stars. Although $\dot{M}\sim {Z}^{m}$ has been included in some recent cool-wind functions (e.g., Mauron & Josselin 2011; Groh et al. 2019), the exponent m is uncertain. Recent observations have shown that the mass-loss rates of RSGs are nearly independent of metallicity (Goldman et al. 2017; Beasor et al. 2020). Therefore, we adopt the dJ88 prescription without metallicity dependence available in MESA, similar to other works such as Choi et al. (2016).

WR stars have T_eff > 10⁴ K, but their mass-loss rates are a factor of 10 higher than those of hot O stars exhibiting the same luminosity (Nugis & Lamers 2002); thus, the hot-wind prescription is not suitable for them. We adopt the WR wind prescription by NL00, which is an empirical law of mass-loss rate as a function of luminosity, helium abundance, and metallicity. It was derived from the observations of 64 Galactic WR stars with the correction of clumping in the wind. This prescription is supported by an optically thick radiation-driven wind theory (Nugis & Lamers 2002).

The hot-, cool-, and WR wind schemes operate under different stellar evolution stages. The switches between different wind schemes depend on T_eff and the surface hydrogen abundance (X_H ) of the simulated stars. The wind prescriptions are explained in detail in the Appendix.

In the "Dutch" hot-/WR wind scheme in MESA, a "low-T" regime is present where the dJ88 prescription is adopted rather than the V01 prescription. For simplicity, we use the term "hot wind" to denote only the mass loss derived from the V01 prescription and use the term "cool wind" to denote all the mass loss derived from the dJ88 prescription.

We establish a set of models to investigate the influence of the initial mass and metallicity on the stellar-mass loss. In our grid, the initial mass values are sampled successively in the range of 1160 M_⊙ with 1 M_⊙ steps. Moreover, 38 different values of initial absolute metallicities Z are chosen: Z = (1, 2,...,9) × 10⁻⁵, (1, 2,...,9) × 10⁻⁴, (1, 2,...,9) × 10⁻³, and (1.0, 1.1, 1.2,...,2.0) × 10⁻². For reference, the solar metallicity is Z_⊙ = 1.34 × 10⁻² (Asplund et al. 2009). Overall, the total number of models in our grid is 1900.

Our stellar models are typically evolved until the iron core-collapse stage. However, some models with initial masses of 11–17 M_⊙ stall at the carbon- or oxygen-burning stages, usually because of abrupt changes in abundance that cannot be easily treated numerically. As these stars have already passed the RSG phase, we adopt the masses at this point as the final masses.

We perform our simulations using the TIARA-TC cluster ⁵ maintained by the Institute of Astronomy and Astrophysics, Academia Sinica (ASIAA). This cluster has 50 nodes, 600 CPU cores, and 2.3 TB memory. Our grid simulations were completed in about 100,000 core-hours.

We perform simulations of stellar evolution with different initial masses and metallicities. Their evolutionary tracks in the theoretical Hertzsprung–Russell (HR) diagram, i.e., luminosity (L) versus T_eff diagram, are shown in the left panels of Figure 1. During the main-sequence phase of 20–60 M_⊙ stars, T_eff decreases by 0.1–0.25 dex and L increases by 0.3–0.15 dex over the main-sequence lifetime of ∼8–4 Myr. A star reaches the terminal-age main sequence (TAMS) when the hydrogen in its stellar core is exhausted. At this time, a hydrogen-burning shell forms, the star starts expanding, and T_eff quickly decreases. Then, the core-helium burning starts, and the outer envelope continues to expand and becomes much cooler. In the HR diagram, the star quickly crosses the "Hertzsprung gap" and moves to the right side. After these post-main-sequence evolution stages, the star becomes a supergiant. Supergiants with $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})\lt 3.6$ (or T_eff ≲ 4000 K) are often called RSGs, $3.6\lt \mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})\lt 3.8$ (or 4000 K ≲ T_eff ≲ 6300 K) are called yellow supergiants (YSGs), and $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})\gt 3.8$ (or T_eff ≳ 6300 K) are called BSGs (e.g., Georgy 2012).

Zoom In Zoom Out Reset image size

A star's lowest T_eff in its evolutionary track depends on its initial mass and metallicity. As shown in Figure 1, for stars with Z = 0.02, those with initial masses ≳50 M_⊙ can only become YSGs instead of RSGs. Furthermore, a comparison of the models of Z = 0.02 and Z = 0.0005 shows that the evolutionary tracks of lower-Z stars usually end at higher T_eff. This feature is critical in the relation between mass loss and metallicity and will be elaborated on in later sections. For convenience, we use the term "cool supergiant" as long as the dJ88 wind prescription is fully active, i.e., T_eff < 10,000 K, or $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})\lt 4.0$.

During the advanced burning stages, stars with initial masses ≳50 M_⊙ can reach luminosities ≳10⁶ L_⊙. These stars can be regarded as luminous blue variables (LBVs) (Humphreys & Davidson 1994), but our simulations do not include eruptive mass loss, which is a common feature of LBVs. Furthermore, some LBVs lose their envelopes and evolve back to a hotter region in the HR diagram to become WR stars. We use the criterion of surface hydrogen abundance (X_S ) below 0.4 to identify WR stars in our work.

As a star evolves through different burning phases, its T_eff and radius vary, resulting in different types of stellar winds and mass-loss rates. As described in Section 2.3, the mass-loss rates are calculated using a combination of mass-loss prescriptions chosen according to T_eff and X_S . In the HR diagrams in Figure 1, most high-Z and some low-Z stars shift from the hot-wind region to the cool-wind region in the post-main-sequence stage. The massive stars of 50 and 60 M_⊙ with Z = 0.02 further evolve back to the hot region and enter the WR phase. The criterion X_S < 0.4 is used to activate the WR wind scheme.

The right panels of Figure 1 present the evolution of the mass-loss rates in our stellar evolution simulations. During the main-sequence phase, $\mathrm{log}{T}_{\mathrm{eff}}$ is in the range of ∼4.45–4.75, and thus, hot winds dominate. For models of the same Z, higher-mass stars have higher luminosities and higher mass-loss rates; furthermore, a main-sequence star's mass-loss rate is usually low and stable because the T_eff and luminosity changes are not large. For models of different Z but the same initial mass, a higher Z leads to a higher mass-loss rate, but the $\dot{M}$ dependence on Z is not a simple Z^m because $\dot{M}$ also depends on L and T_eff, which are also dependent on Z. For Z = 0.0005, $\dot{M}$ for 20–60 M_⊙ main-sequence stars is in the range ∼10⁻⁸–10⁻⁶ M_⊙ yr⁻¹, and for Z = 0.02, ∼10^−7.5–10^−5.5 M_⊙ yr⁻¹. When a star evolves into the post-main-sequence stage, T_eff suddenly decreases. The change in T_eff causes the enhancement of the mass-loss rate by about 13 orders of magnitude because the T_eff passes through two jumps in the mass-loss rate included in the wind prescriptions. These two jumps correspond to two physical transitions in the supergiant winds:

(1) The "bi-stability jump" in early-type supergiant winds. When the T_eff is less than ∼25,000 K, the mass-loss rate suddenly increases and the wind velocity decreases, which is called the "bi-stability jump" (Pauldrach & Puls 1990). This is caused by the recombination of Fe iv. When the temperature decreases below ∼25,000 K, the ionization fraction of Fe iii below the sonic point significantly increases, and Fe iii lines are more efficient drivers of winds than Fe iv lines (Vink et al. 1999). The mass-loss prescription of V01 illustrates the bi-stability jump by setting different mass-loss functions in the two regimes with T_eff values higher and lower than the jump temperature, which ranges between 22,500 and 27,500 K.

(2) The transition from hot winds to cool supergiant winds. When the T_eff decreases below 10,000–12,000 K, the star evolves into a cool supergiant. The cool, dust-driven wind becomes the dominating means of mass loss, and the mass-loss rate significantly increases.

These two jumps in the mass-loss rate occur at the post-main-sequence stage, marked as the bi-stability jump and the hot-to-cool-wind transition region in Figure 1. The stellar evolution timescale between these two jumps is usually so short, <10⁴ yr, that these two jumps are unresolvable in the evolution of mass-loss rate plots (the right panels of Figure 1).

Nevertheless, for some low-Z stars with Z = 0.0005, their T_eff only decreases to the range between 10,000 and 25,000 K at the onset of the post-main-sequence stage, and they remain as early-type supergiants. Some of the other low-Z stars do cool down to T_eff < 10,000 K, but this happens in the very late evolutionary stage after core-carbon-burning commences; thus, the separation between the two jumps is about several hundred thousand years. For example, for the star with 30 M_⊙ and Z = 0.0005, the second jump in the mass-loss rate occurs at the very end of its lifetime.

Figure 2 displays representative relations between a star's initial mass (M_i ) and the total mass loss (ΔM), which is the integrated mass loss over the stellar lifetime. The total mass loss generally increases with the initial stellar mass, but the detailed relation depends on the metallicity. At first glance, three different trends of mass loss are observed at different metallicities:

(1) For Z ≳ 5 × 10⁻³, the high-Z cases, the total mass loss of a star generally increases with its initial mass. For example, for Z = 0.02 (>1 Z_⊙ = 0.0134), a 15 M_⊙ star loses ∼2 M_⊙ and a 60 M_⊙ star loses >30 M_⊙ during its lifetime.

(2) For Z ≲ 3 × 10⁻⁴, the low-Z cases, the total mass loss of a star also increases with its initial mass, but it is considerably lower than that of the Z ≳ 10⁻² cases. For example, for Z = 10⁻⁴, a 15 M_⊙ star loses ∼0.01 M_⊙, and a 60 M_⊙ star loses ∼0.5 M_⊙. Some low-Z high-mass stars have high mass losses, but their total mass loss is still <20% of the initial mass.

(3) For stars with Z ∼ 10⁻³, the total mass loss is not a smooth function of the initial mass; instead, it oscillates between the high values for the high-Z cases and the low values for the low-Z cases.

Zoom In Zoom Out Reset image size

We have shown that the total mass loss during a star's lifetime largely depends on its metallicity. To further illustrate the metallicity effects on mass loss, we plot the total mass loss as a function of metallicity in Figure 3.

Zoom In Zoom Out Reset image size

Figure 3 displays a significant jump in the total mass loss at Z ∼ 10⁻³. When Z ≲ 3 × 10⁻⁴, the mass loss is generally low. When Z approaches ∼10⁻³, the mass loss jumps to a significantly higher level.

Herein, we call the Z ∼ 10⁻³ the critical metallicity (Z_c) for the mass loss of massive stars. Generally, significant mass loss occurs only when Z > Z_c. Nonetheless, Z_c is not a sharp boundary as shown in the stellar initial mass versus metallicity (M_i − Z) diagram of mass-loss fraction (ΔM/M_i) in Figure 4, where all 1900 models are included. Two distinct regions are present in this diagram: the high-Z region with a mass-loss fraction ≳40% and the low-Z region with a mass-loss fraction ≲10%. The separation of the high-Z and low-Z regions is clear, and the boundary at Z ∼ 10⁻³ is what we call Z_c. For initial masses ≲30 M_⊙, the boundary between these two regions is sharp; for initial masses ≳30 M_⊙, the boundary is fuzzy.

Zoom In Zoom Out Reset image size

To understand the origin of Z_c, we compare the individual contributions of mass loss by hot, cool, and WR winds in order to identify a specific stellar evolutionary stage that is most accountable for Z_c. Figure 5 displays the contributions of different winds for initial stellar masses of 25 and 50 M_⊙. In both cases, a distinct jump in mass loss occurs only in the cool wind near Z_c, indicating that the cool wind is responsible for Z_c.

Zoom In Zoom Out Reset image size

The M_i−Z diagrams of the mass-loss fraction for hot, cool, and WR winds are presented in Figure 6. For cool winds, going from high Z to low Z, the high mass loss makes a sharp transition into a low mass loss near a threshold Z of ∼1 × 10⁻³. The threshold Z is slightly lower for stars with M_i = 25–30 M_⊙ and higher for the more massive stars. For hot winds, the transition of high mass loss to low mass also occurs near a threshold Z of 1 × 10⁻³, but the transition is gradual. For WR winds, the threshold Z is higher, near 3 × 10⁻³. As the total contribution of hot and WR winds to the mass loss is small, ΔM/M_i ≲ 20%, the M_i−Z diagram of cool wind is very similar to that of the total mass loss (all winds combined) in Figure 4; furthermore, the Z_c of the total mass loss is similar to the threshold Z of the cool winds. It should be noted that at low Z, the hot winds may lose more mass than the cool winds up to the threshold Z of the cool winds; however, the mass lost via hot winds is so much smaller than the mass lost via cool winds that the cool winds' threshold Z dominates.

Zoom In Zoom Out Reset image size

Therefore, Z_c of mass loss mainly stems from cool winds. This result is reasonable because the mass-loss rates of cool winds are usually considerably higher than those of hot winds, as included in the wind prescriptions (Section 3.2).

We have shown that the mass-loss fraction (ΔM/M_i) of a massive star is high for high Z and low for low Z, making a transition at Z_c that is caused mainly by the cool winds. To gain more insights into the dichotomy of mass-loss fraction and its scatter around Z_c, we use 25 M_⊙ and 50 M_⊙ stars with different metallicities that are near Z_c as examples and look into their post-main-sequence evolution. For each stellar mass, five models with Z = 0.0005, 0.0008, 0.001, 0.002, and 0.005 are examined. These models are marked by white dots in Figure 4. It can be seen that the 25 M_⊙ stars show a monotonic increase of mass-loss fractions with increasing Z, while the 50 M_⊙ stars show nonmonotonic variations of the mass-loss fraction across Z_c. We deliberately choose these two masses with contrasting mass-loss variations across Z_c and expect the comparison between them may shed light on the physical origin of the transition and the scatter of the mass-loss fraction near Z_c.

To investigate the crucial differences among these models, we plot their stellar evolutionary tracks and mass-loss rates in Figure 7. For the 25 M_⊙ stars, the three models yielding high ΔM (Z = 0.001, 0.002, and 0.005) are those that evolve to radii R > 1000 R_⊙ and $\mathrm{log}({T}_{\mathrm{eff}}/{\rm{K}})\lt 3.6$ (i.e., T_eff < 4000 K). In other words, these stars end their lives as RSGs. Due to stellar expansion, their mass-loss rates increase by two to three orders of magnitude in the supergiant stages. In contrast, the two models yielding low ΔM (Z = 0.0005 and 0.0008) only expand to ∼100 R_⊙. After moderate expansions, they bounce back to the left (i.e., higher T_eff) in the HR diagram, and their T_eff remains higher than ∼10,000 K. Such a high T_eff prevents the cool wind from operating, and thus, their mass-loss rates only slightly increase during the post-main-sequence stages. In summary, the two different behaviors of stellar expansion result in the dichotomy of total mass loss.

Zoom In Zoom Out Reset image size

For the 50 M_⊙ stars, the evolutionary tracks are more complicated. At the age of 4 Myr (marked by "+" in the HR diagram in Figure 7), the two models yielding high ΔM (Z = 0.0005 and 0.005) have already evolved into the cool-wind domain and eventually expand to R > 1000 R_⊙. On the other hand, the three models yielding low ΔM (Z = 0.0008, 0.001, and 0.002) are still in the hot-wind domain at 4 Myr. These three models eventually enter the cool-wind domain at the very end of their lifetimes when the core-carbon burning commences, and approach a high mass-loss rate of ∼10⁻⁴ M_⊙ yr⁻¹. However, their cool winds only last a few ×10⁴ yr. Thus, the integrated ΔM is still small, yielding ΔM/M_i < 10%.

The above examples demonstrate that the crucial process that determines the total mass loss is the expansion toward the formation of a supergiant. If a star expands to a cool supergiant with R > 1000 R_⊙ at the beginning of the post-main-sequence stage, significant mass loss will occur. Otherwise, if a star never expands to the extent of a cool supergiant or only expands to that extent at a very late stage (i.e., after the core-carbon burning begins), only a minimal amount of mass loss will occur.

The crucial role played by the supergiant expansion in the dichotomy of mass loss is further illustrated in Figure 8, which presents maps of maximum radii in the M_i−Z grid. We have considered maximum radii over two time spans: the entire lifetime (left panel of Figure 8) and the time span from the ZAMS to the beginning of core-carbon burning (right panel of Figure 8). The left panel shows a complex map with a visible dichotomy only for stellar masses of 14–27 M_⊙. The complexity stems from the maximum radii arising from different evolutionary stages and paths for stars of different masses. For high metallicities (Z > Z_c), stars reach maximum radii during the core-helium burning, while for low metallicities (Z < Z_c), 14–27 M_⊙ stars reach the maximum radii during the core-helium burning, but the higher-mass stars reach maximum radii after the core-carbon burning. The core-carbon-burning stage is pretty late in a massive star's lifetime, and the duration for the high mass-loss rate is short; thus, massive stars that reach maximum radii during the core-carbon-burning stage have small mass-loss fractions. The complications at the core-carbon-burning stage make the maximum radius over the entire lifetime an ineffective indicator of total mass loss.

Zoom In Zoom Out Reset image size

As the core-helium-burning stage lasts the longest after the main sequence, and the bulk of mass loss occurs during this stage, we next consider the maximum stellar radii during the time span from the ZAMS to the beginning of core-carbon burning. In the right panel of Figure 8, a cleaner dichotomy is seen and the transition occurs near Z_c ∼ 0.001. In fact, this map of maximum stellar radii resembles closely the map of the mass-loss fraction (Figure 4): Both show dichotomy across the Z_c ∼ 0.001 and a larger scatter around Z_c at high stellar masses. Examined closely, even the scatter exhibits one-to-one correspondences between these two maps. These similarities indicate that the maximum stellar radius before core-carbon burning is an effective indicator of total mass loss.

We further find that models with high mass loss are almost exactly those that have expanded to cool supergiants with R > 1000 R_⊙ at the core-helium-burning stage, and these models are usually of Z ≳ 10⁻³. Therefore, we conclude that the Z_c of mass loss is essentially the critical metallicity of cool supergiant formation.

Generally, cool supergiant formation and the corresponding mass loss only occur when Z is higher than the critical metallicity (Z_c ∼ 10⁻³). For stars with Z > Z_c, they usually expand to R > 1000 R_⊙ and become cool supergiants at the beginning of the post-main-sequence stage. These stars can lose large amounts of mass due to cool winds. In contrast, stars with Z < Z_c either do not reach R > 1000 R_⊙ throughout their lifetimes or only expand to this extent at the very final stage of stellar evolution. For these low-Z stars, the cool-wind scheme is not activated in most of their stellar lifetime; thus, their total mass loss is very low. For Z ∼ Z_c, a star either succeeds or fails to become a cool supergiant, and thus shows erratic variations of the total mass loss as a function of metallicity (Figures 2 and 4). The physical mechanism that determines whether a star can become a cool supergiant will be investigated in our Paper II.

The flux of kinetic energy carried by a stellar wind is called the mechanical luminosity of the wind, defined as ${L}_{{\rm{w}}}\equiv \tfrac{1}{2}\dot{M}{v}_{\infty }^{2}$, where v_∞ is the wind terminal velocity. The determination of $\dot{M}$ has already been described in Section 2. Our stellar evolution models do not provide v_∞ directly but provide escape velocity (v_esc) that can be scaled to determine v_∞. For hot winds, we follow V01 and directly use the v_∞/v_esc scaling relations given in Equations (A1) and (A2) in the Appendix. For cool winds, we follow Leitherer et al. (1992) and set the RSG wind terminal velocity at a constant 30 km s⁻¹. More details of v_∞/v_esc can be found in the Appendix. A star's $\dot{M}$ and v_∞ are used to compute L_w, then L_w is integrated over the star's lifetime to determine the total kinetic energy injected by its winds into the ISM.

Figure 9 illustrates v_∞ and L_w of 25 M_⊙ stars with different Z values around Z_c. Main-sequence O stars produce fast winds of ∼1500–3000 km s⁻¹. When a star evolves into an RSG, its v_∞ suddenly drops by a factor of 100. At the same time, its $\dot{M}$ goes up by two to three orders of magnitude (Figure 7). As a result, L_w decreases by one to two orders of magnitude when the star evolves from a main-sequence O star to an RSG. Although cool RSG winds release large amounts of mass, they only carry relatively small amounts of kinetic energy because of their slow wind velocities. For stars with higher initial masses and metallicities, the domination of hot winds in the energy output is even more significant.

Zoom In Zoom Out Reset image size

Figure 10 presents the total kinetic energy output over the stellar lifetime of each model. This figure also shows the contributions of the hot and cool winds separately. It is evident that the hot winds dominate the energy output, while the contributions from cool winds are almost negligible. The dominance of hot winds is attributed to their velocities being much faster than the cool winds' velocities, and their duration (including the main-sequence stage) spanning ∼90% of the stellar lifetime. It can also be seen from the top panel of Figure 10 that the most massive stars with near-solar metallicities each inject 10⁵⁰–10⁵¹ erg of kinetic energy into the surrounding ISM, which is comparable to the explosion energy of an SN. In contrast, the cool winds of such a massive star contribute only up to ∼10⁴⁷ erg to the total kinetic energy.

Zoom In Zoom Out Reset image size

It is interesting to note that the Z_c for the dichotomy of total mass loss shows some effects only on the cool winds' total kinetic energy. This is understandable because Z_c is closely associated with the onset of cool winds, as discussed in Section 4.3. The stellar winds' total kinetic energy injected into the ISM is dominated by contributions from hot winds. As cool winds play a negligible role in the total wind energy, its closely associated Z_c does not affect the total wind energy or total hot-wind energy. In other words, whether a star becomes an RSG has no bearing on the total wind kinetic energy injected into the ISM.

The mass loss and kinetic energy output from stellar winds of individual stars can be used to evaluate the feedback from a star cluster. We consider a star cluster from a single burst of star formation with an initial mass function (IMF) ξ(M) = ξ₀ M^−α , where M is the initial stellar mass, ξ₀ is a scaling factor, and α is a constant. Here, ξ₀ is scaled by the number of stars in the system N so that ${dN}={\xi }_{0}{M}^{-\alpha }{dM}$. In our calculations, we apply α = 2.35 (the Salpeter IMF; Salpeter1955). For comparison, the widely used IMF in Kroupa (2001) gives α = 2.3 for stars with masses >1 M_⊙, which is similar to the Salpeter IMF.

We have calculated the feedback for a cluster with a total mass of 10⁵ M_⊙. Although this mass is associated with a super star cluster, the feedback scales linearly with the cluster mass and can be used to determine feedback of clusters of any other masses. For the cluster members, we adopt the canonical initial stellar mass range of 1–150 M_⊙. We note that the upper mass limit is not important because the population of very massive stars is much smaller than that of lower-mass stars. With the above assumptions, we integrate the $\dot{M}$ and L_w of individual stars to calculate the wind feedback from the entire cluster. The feedback we present only includes the contributions from massive stars of M_i = 11–60 M_⊙, which is the mass range of our simulations; thus, our results show the lower limit of wind feedback.

In Figure 11, we choose a high-Z (Z = 0.02) case and a low-Z (Z = 0.0002) case and plot their integrated mass-loss rates and energy-injection rates of a 10⁵ M_⊙ cluster before the age of 10 Myr. At about 3–4 Myr, the most massive stars evolve into the post-main-sequence stage. In the high-Z cluster, these massive stars successively become cool supergiants and significantly enhance the overall mass-loss rate. However, the lifetimes of cool supergiants are short, so each star can only contribute to the enhancement of the mass-loss rate for a short time. When stars of lower masses become cool supergiants, the more massive stars have already ended their lifetimes. The overall trend of the mass-loss rate is a gradual decrease after the sudden enhancement at about 3–4 Myr. The jagged appearance of the curve in the mass-loss rate is caused by the limit of mass resolution (1 M_⊙) in our simulations. For the low-Z cluster, the stars do not evolve into cool supergiants at the helium-core-burning stage, so the mass-loss rate after 4 Myr can only increase to 10⁻⁵ M_⊙ yr⁻¹ due to the bi-stability jumps. Some sudden spikes of the mass-loss rate occur at about 4 Myr because some low-Z stars reach cool supergiants in the very late stage of their lifetimes, and thus they enhance the mass-loss rate for very short times.

Zoom In Zoom Out Reset image size

While the mass-loss rate of the high-Z cluster is enhanced by cool winds at the ages of 3–4 Myr, cool winds do not contribute much to the energy-injection rate. Thus, for both high- and low-Z clusters, their energy-injection rates do not rise drastically but continue to decline after 3–4 Myr, when the most massive stars end their lifetimes.

We compare our mass-loss rates and energy-injection rates with those adopted by the FIRE-3 cosmological simulations (Hopkins et al. 2023), which are shown by the blue dashed lines in Figure 11. The mass-loss rate in FIRE-3 is expressed by a function of time and Z with four segments separated by 1.7, 4, and 20 Myr. For the high-Z model, the post-main-sequence mass-loss rate from our simulations is in reasonable agreement with the FIRE-3 curve, while for the low-Z model, our mass-loss rate after ∼4 Myr is about one order of magnitude lower than the FIRE-3 rate. The essential difference is because the low-Z stars in our models do not evolve into cool supergiants when they enter the post-main-sequence stage, so the mass-loss rate is low. For the times before 4 Myr, the comparison of mass-loss rates is out of our scope as we do not consider the very massive stars with >60 M_⊙. We also compare our energy-injection rate to the FIRE-3 function. For both high- and low-Z cases, our energy-injection rate is lower than the FIRE-3 rate by about an order of magnitude after ∼4 Myr. FIRE-3 uses the entire mass-loss rate and the average velocity, which is ∼1000 km s⁻¹ at ∼4 Myr, to estimate the wind feedback. In contrast, we separately consider hot and cool winds, with hot winds having velocities of >1000 km s⁻¹ but low mass-loss rates and cool winds having high mass-loss rates but a velocity of 30 km s⁻¹. Thus, our calculations lead to a lower energy-injection rate than the FIRE-3 function.

We have calculated the feedback from stellar winds in a cluster using the Salpeter IMF. The critical metallicity of cool supergiant formation significantly affects the total mass loss within a cluster, and we suggest that this effect be considered in order to refine the treatment of stellar feedback in cosmological simulations. In contrast, the kinetic energy feedback is dominated by hot winds, which are irrelevant to the critical metallicity of cool supergiants.

We further integrate the mass and energy injected by a 10⁵ M_⊙ cluster over time, and the results are shown in Figure 12. The contribution of wind feedback in our calculations is from stars with M_i = 11–60 M_⊙. For comparison, we also plot the feedback that only includes the contribution of 11–30 M_⊙ stars. From these plots, hot winds dominate the kinetic energy and momentum feedback for all of the metallicities, while cool winds surpass hot winds in the contribution to mass loss when Z > Z_c. The effect of critical metallicity occurs in the output of mass, but not in that of kinetic energy. The WR wind is never the dominating component, although its contribution to kinetic energy is larger than that of the cool wind in an environment of solar metallicity. We note that these integrated quantities are the lower limit, as we do not consider the very massive stars with >60 M_⊙.

Zoom In Zoom Out Reset image size

Our study finds that the critical metallicity of cool supergiant formation causes a significant increase in mass loss. This effect arises from stellar evolution rather than the treatment of mass-loss prescription. We discern two different effects of metallicity on RSG evolution and review the previous theoretical works that have shown a similar result of critical metallicity.

From the evolutionary tracks shown in Figure 7, there are two apparent effects of metallicity on the evolution of cool supergiants. First, at a lower metallicity, the overall evolutionary tracks, from the ZAMS to the supergiant stage, shift to higher T_eff or bluer locations. Second, the evolutionary tracks of stars with Z < Z_c end at hotter T_eff, whereas stars with Z > Z_c can successfully evolve into cooler regimes. This is a sharp transition rather than a continuous change. In the following paragraphs, we discuss these two different effects.

(1) The overall shift of the evolutionary tracks with metallicity. This effect has been discussed with regard to the two ends of stellar evolution: the shift of ZAMS and the color of supergiants. First, the displacement of ZAMS toward the bluer locations with decreasing Z has been reported by stellar models (e.g., Schaller et al. 1992; Mowlavi et al. 1998; Groh et al. 2019). Among several physical parameters that directly change with metallicity (opacity, energy generation rate, molecular weight, etc.), the ZAMS shift has been attributed to opacity κ (see, e.g., Maeder 2009). In general, decreasing Z leads to lower κ, transferring more radiation outward and leading to higher luminosities, and the stars also need to contract further to maintain the energetic equilibrium. With a larger L and a smaller R, a higher T_eff(∝ L/R²) is required for lower-Z stars.

At the end of stellar evolution, supergiants with lower Z also shift toward bluer locations in the HR diagram. This feature has been observed in stars in the Local Group (Humphreys 1979a,1979b; Elias et al. 1985) and has been reproduced in stellar models (e.g., Schaller et al. 1992; Meynet et al. 1994; Groh et al. 2019). Surveys of massive stars have confirmed that the RSGs in galaxies with lower metallicities have earlier spectral types on average: M2 I in the Milky Way, M1-M1.5 I in the LMC, and K3K7 I in the Small Magellanic Cloud (SMC) (Massey & Olsen 2003; Levesque et al. 2006). This can be explained as the shift of the Hayashi line, which represents the coolest extent of giant stars, to higher T_eff at lower Z (Elias et al. 1985; Levesque et al. 2006). These stars shift to bluer positions at lower Z due to opacity (Hayashi & Hoshi 1961; Kippenhahn et al. 2013), which is a similar physical effect to the ZAMS shift.

(2) The critical metallicity of cool supergiant formation. Instead of causing continuous shifts in the HR diagram, the critical metallicity acts as a boundary between two different types of evolution tracks in the supergiant stage. The stars with Z > Z_c successfully become RSG with R ≳ 1000 R_⊙, and those with Z < Z_c remain as smaller BSGs. This effect causes a significant increase in mass loss.

Some previous models have shown that low-Z stars can never evolve into RSGs (Arnett 1991; Baraffe & El Eid 1991; Brocato & Castellani 1993; Hirschi 2007; El Eid et al. 2009; Limongi 2017; Groh et al. 2019). The metallicity value has been identified as Z ∼ 10⁻³ (Baraffe & El Eid 1991;El Eid et al. 2009), which is consistent with our results. We start with the investigation of mass loss and identify the same phenomenon in cool supergiant formation.

However, the physical origin of the critical metallicity of cool supergiant formation is still poorly understood. A profound question remains: Why do massive stars of Z > Z_c successfully expand to RSGs, but those of Z < Z_c fail to become RSGs? In the companion Paper II, we will present some experiments that study the physical origin of critical metallicity. The explanation of its mechanism will help us understand the fundamental physics of RSG formation.

In this subsection, we discuss the uncertainties in our models from both stellar models and mass-loss prescriptions. We then evaluate the impacts of these factors on our main results.

6.2.1. Uncertainties from Stellar Models

6.2.2. Uncertainties from Mass-loss Prescriptions

Our study reveals a critical metallicity of cool supergiant formation that leads to a significant difference in mass loss between high- and low-Z environments. The low-Z stars do not evolve into cool supergiants, and thus only have minimal mass loss through steady-state winds over their lifetimes. In the following paragraphs, we discuss the consequences of critical metallicity.

It is widely known that pre-SN evolution and mass loss set up the physical condition for SN explosions. The types of progenitor stars often determine the types of SNe. For example, RSGs are the progenitors of SNe II-P, as identified in pre-SN images (see Smartt 2009, 2015 for review). It has also been shown that the types of SNe can be significantly affected by the interaction with CSM produced by RSGs (e.g., van Loon 2010; Ekström et al. 2012; Smith 2014; Smartt 2015; Beasor et al. 2020, 2021; Moriya 2021). For example, the early-phase light curves of SNe II-P can be better explained if dense CSM is considered (Moriya et al. 2011, 2017, 2018). Moreover, it is suggested that CSM not only affects SN types but also makes SN explosions more likely to happen (Morozova et al. 2018).

Based on these understandings and the critical metallicity we found, we can expect the low-Z stars to behave differently from high-Z stars in SN explosions. First, if the low-Z stars indeed end up as BSGs instead of RSGs, whether they will still explode as SNe is uncertain. Even if they successfully explode, the resulting SN type may differ from SNe II-P. Low-Z stars never go through the RSG phase, or only become RSGs in the very final stage, so they only lose a negligible amount of mass and may not end as typical SNe II-P. Their circumstellar material can be too dilute for interacting SNe, although we have not considered eruptive mass loss that can generate dense CSM.

On a larger scale, the critical metallicity also impacts the stellar feedback in a star cluster. From our results, the total mass returned to the ISM by stellar winds shows a significant increase at Z_c. In the high-Z (Z > Z_c) domain, these results are consistent with the current treatment of feedback in cosmological simulations such as FIRE-3. However, if the low-Z (Z < Z_c) stars do not evolve into cool supergiants, as our simulations predict, the current treatment of feedback may overestimate the mass injected into the ISM in these low-Z environments. This discrepancy exists only in mass injection but not in energy output, as the kinetic energy feedback is dominated by hot winds, which are irrelevant to the critical metallicity.

As the critical metallicity Z_c ∼ 0.001 is lower than 0.1 Z_⊙, it is still difficult to test it using current observational instruments. Nevertheless, it will be promising to test the stellar evolution below this metallicity with the next-generation telescopes. For example, observations of massive stars in the low-Z galaxies are listed as one of the scientific goals of the LUVOIR observatory (Garcia et al. 2021).

In future works, we plan to apply these 1D mass-loss models to 2D and 3D radiation hydro simulations to further understand how mass loss shapes the CSM and ISM around massive stars. With the advancement of numerical models and observational data, we may reveal the mystery of the mass loss of massive stars and their feedback.

The settings of mass-loss rate in our simulations follow the wind schemes established in MESA. We apply the "Dutch" scheme for hot and WR winds and the dJ88 prescription for cool winds. The criteria to activate these winds are as follows:

• 1.  
  If T_eff > 12,000 K, the "Dutch" scheme (hot/WR wind) is fully operating.
• 2.  
  if T_eff ≤ 8000 K, the cool wind is fully operating, and the prescription of dJ88 is used.
• 3.  
  Between 80,000 K and 12,000 K, a linear interpolation of the dJ88 (cool) and "Dutch" (mainly hot/ WR) winds is applied.

The "Dutch" scheme is not a single function, but a combination of V01's hot-wind prescription, NL00's WR wind prescription, and the extension of dJ88's cool-wind prescription. In detail, the wind configuration of the "Dutch" scheme is as follows:

• 1.  
  If T_eff ≥ 11,000 K and X_H ≥ 0.4, V01's hot-wind prescription is adopted.
• 2.  
  If T_eff ≥ 11,000 K and X_H < 0.4, NL00's WR wind prescription is adopted.
• 3.  
  If T_eff ≤ 10,000 K (low-T "Dutch" scheme), the dJ88 prescription, which is also used to express the cool wind, is adopted.
• 4.  
  If 10,000 ≤ T_eff ≤ 11,000 K, a linear interpolation of the dJ88 and V01/NL00 prescriptions is adopted.

There are criteria to turn off the mass loss at the very end of the stellar lifetime. If the central temperature is higher than 2 × 10⁹ K, mass loss will be fully turned off. If the central temperature is between 1 × 10⁹ and 2 × 10⁹ K, the mass-loss rate weakens linearly with the central temperature. These settings simply follow the default in MESA.

As our MESA simulations do not include wind models, the wind terminal velocity (v_∞) is estimated using empirical functions of v_esc and other stellar parameters.

• 1.  
  Hot wind. If V01's prescription of mass loss is turned on, we follow the wind velocity formulae adopted by Vink et al. (2001). These expressions originate from Leitherer et al. (1992) and Lamers et al. (1995).

  • (a)  
    If T_eff > 27,500 K, which is higher than the bi-stability jump, then

    $\begin{eqnarray}&&{v}_{\infty }/{v}_{\mathrm{esc}}=2.6{(Z/0.019)}^{0.13}.\end{eqnarray} \tag{ A1 }$
  • (b)  
    If T_eff < 22,500 K, which is lower than the bi-stability jump, then

    $\begin{eqnarray}&&{v}_{\infty }/{v}_{\mathrm{esc}}=1.3{(Z/0.019)}^{0.13}.\end{eqnarray} \tag{ A2 }$
  • (c)  
    If 22,500 < T_eff < 27,500, the mass-loss rate is the linear combination of the high-T mass-loss rate (${\dot{M}}_{h}$) and the low-T mass-loss rate (${\dot{M}}_{l}$):

    $\begin{eqnarray}&&\dot{M}=\alpha {\dot{M}}_{h}+(1-\alpha ){\dot{M}}_{l},\end{eqnarray} \tag{ A3 }$

    where α ≡ (T_eff − 22,500)/500. The wind velocity we adopt is the weighted average of the high-T wind velocity (v_∞,h ) and the low-T wind velocity (v_∞,l ), and weighting coefficients are based on the ratio of $\dot{M}$ from high- and low-T components. The weighted average velocity ${\bar{v}}_{\infty }$ can be expressed as

    $\begin{eqnarray}&&{\bar{v}}_{\infty }=\displaystyle \frac{\alpha {\dot{M}}_{h}{v}_{\infty ,h}+(1-\alpha ){\dot{M}}_{l}{v}_{\infty ,l}}{\alpha {\dot{M}}_{h}+(1-\alpha ){\dot{M}}_{l}}.\end{eqnarray} \tag{ A4 }$
• 2.  
  Cool wind. We follow Leitherer et al. (1992) to set the cool-wind velocity at a constant 30 km s⁻¹. As long as the dJ88 prescription is used, even if it is called from the "Dutch" low-T scheme, we regard this part of mass loss as the cool wind.
• 3.  
  WR wind. We use the wind velocities given by Nugis & Lamers (2000) for WR stars.

  • (a)  
    WN stars:

    $\begin{eqnarray}&&\mathrm{log}({v}_{\infty }/{v}_{\mathrm{esc}})=0.61-0.13\mathrm{log}L+0.3\mathrm{log}Y.\end{eqnarray} \tag{ A5 }$
  • (b)  
    WC stars:

    $\begin{eqnarray}&&\mathrm{log}({v}_{\infty }/{v}_{\mathrm{esc}})=-2.37-0.43\mathrm{log}L-0.07\mathrm{log}Z.\end{eqnarray} \tag{ A6 }$

In some transition temperatures, the mass-loss rate is composed of more than one wind component. In such cases, we use the ratios of the mass-loss rates of these components as the weighting coefficients to calculate the average wind velocity. The expression is similar to Equation (A4), but with a linear combination of the hot-/WR and cool-wind velocities.