On the Maximum Black Hole Mass at Solar Metallicity

The LIGO/Virgo/KAGRA (LVK) collaboration has delivered a summary of detections of compact object mergers in gravitational waves in three observing runs O1/O2/O3 (The LIGO Scientific Collaboration et al. 2023). The majority of these detections are black hole–black hole (BH–BH) mergers. BHs are found in wide mass range ∼2.6–95 M_⊙. LVK inferred an intrinsic BH mass distribution for more-massive BHs in BH–BH mergers. The distribution decreases steeply with the primary BH mass and it shows two peaks at ∼8 and ∼35 M_⊙. Understanding such a result is a complex task. The number of evolutionary channels may lead to the formation of BH–BH mergers (e.g., Mandel & Broekgaarden 2022), and each channel has its own inherent uncertainties (e.g., Belczynski et al. 2022). In many of the formation channels stars are the building blocks for BH–BH mergers. There is therefore the need for precise predictions from stellar evolution in regard to final state of massive stars at the time of core collapse. Such predictions allow, among other things, to estimate stellar-origin BH masses. The BH mass is predicted to depend sensitively on metallicity (Belczynski et al. 2010a, 2010b). Yet, a good starting point for such calculations is solar metallicity, because it represents a condition at which the effect of wind mass loss is highly noticeable, since many of the stellar wind models prescribe wind-mass-loss rates as a function of chemical abundances and luminosity, which both depend on metallicity levels. Additionally, solar metallicity is one of the most studied and better constrained cases. For instance, the most massive stellar BH found in the Milky Way is Cyg X-1, with M_BH ≃ 21.2 ± 2.2 M_⊙ (Miller-Jones et al. 2021), this value being the lower threshold for theoretical estimations of ${M}_{\mathrm{BH},\max }$.

Recently, two studies (that we are using as a comparison) offered new results in terms of final masses of massive stellar models. The final star mass may be taken as a proxy for BH mass (direct BH formation) or as an upper limit on BH mass (if there is a supernova explosion involved in the BH formation). First, Bavera et al. (2023) performed stellar evolution of massive stars (both as single stars and as parts of a binary system) in the zero-age main-sequence (ZAMS) star mass range M_ZAMS = 15–150 M_⊙ at Z = 0.014 with the MESA evolutionary code. They found that the BH mass from single star evolution increases from M_BH = 3 M_⊙ at M_ZAMS = 15 M_⊙ to a peak value at M_BH = 35 M_⊙ at M_ZAMS = 75 M_⊙, and then BH mass rapidly decreases with initial mass to M_BH = 15 M_⊙ at M_ZAMS = 100 M_⊙. They also found that for more-massive stars (M_ZAMS > 100 M_⊙) the BH mass does not change and remains more or less constant (M_BH = 14 M_⊙). Second, Martinet et al. (2023) provided models performed with the Genec (Geneva evolution code). We only pick the rotating models at Z = 0.014 from this study for our comparison. These include M_ZAMS = 180, 250, 300 M_⊙. Final masses (end of C burning) of these models are M_final = 42, 28, and 37 M_⊙. Additionally these models at final stage are WC/WO spectral subtype stars. For such final masses we may expect that they approximately correspond to BH masses. On one hand, for such high final masses supernova models predict direct BH formation (Fryer et al. 2012); on the other hand, pulsational pair-instability supernovae are predicted to start removing mass for helium cores more massive than ∼35–40 M_⊙ (Woosley 2017).

Both of the above calculations have employed "Dutch"-like wind-mass-loss prescriptions, which have been considered for years the standard setup in 1D stellar evolution. The Dutch wind prescription orbits around the use of three different mass-loss models: one for thin winds (Vink et al. 2001), one for thick winds where the switch from thin to thick winds is placed on the mass fraction of surface hydrogen X_surf > 0.4 (Nugis & Lamers 2000), and one for dust-driven winds at effective temperature T_eff < 10 kK (de Jager et al. 1988). However, recent years offered various detailed studies of massive stars and their mass loss through stellar winds (e.g., Krtička & Kubát 2017; Björklund et al. 2021; Gormaz-Matamala et al. 2023). We take advantage of these new results to update our calculations. Additionally, in Bavera et al. (2023) only nonrotating models were used, and in Martinet et al. (2023) a rather modest overshooting was employed (α_ov ∼ 0.2).

Here we propose an updated collection of wind-mass-loss prescriptions that can be applied to all massive stars that are expected to form BHs. We also advocate for the use of rotating models to better represent the reality of stellar populations and the adoption of high overshooting values (α_ov ∼ 0.5) for massive stars as showed by Scott et al. (2021). We perform most of our models with our updated version of MESA, and we cross-check our MESA results with a couple of calculations performed with our Genec models. Both codes are made to run possibly the same physics (same winds, rotation, convection/overshooting, angular momentum, and chemical transport).

In our study we use two different detailed evolutionary codes, version 23.05.1 of Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023), and Genec (Eggenberger et al. 2008) for stars at solar metallicity Z = 0.014 (Asplund et al. 2009) within an initial mass range between 10 and 300 M_⊙. We unified the key input physics in both codes. We stop all our simulations at the depletion of carbon in the stellar core, since the timescale between the end of core C burning and the core collapse is negligible in terms of star and core mass evolution (see, e.g., Woosley et al. 2002), and it does not therefore change the final estimates for BH masses.

2.1. Mass-loss Prescription

The theoretical calculation of wind-mass-loss rates ($\dot{M}$) depends on the physical properties of a star, and stellar winds change with the spectral type and evolutionary stage changes. In Figure 1 we show our adopted wind-mass-loss prescription. Our choices are motivated by the recent work on various phases of wind mass loss for massive stars. Yet we use older prescriptions in parts of the parameter space where no updates are available or when the change from an older prescription to a more recent one would not be sufficiently motivated.

Zoom In Zoom Out Reset image size

For line-driven winds, which are the dominant winds for the mass range we are evaluating in this work, we make the division between optically thin and optically thick winds.

Optically thin winds models imply that the stellar radius is well defined (such as OBA-type stars), whereas optically thick wind implies that there is not a clear separation between the photosphere and the expanded atmosphere, such as in Wolf–Rayet (W-R) stars. For stellar evolution models, this distinction is made when the star is hot and X_surf is below certain value, normally X_surf < 0.3 or <0.4 (Yusof et al. 2013; Köhler et al. 2015). However, very massive stars such as WNh develop thick winds in early stages of their evolution (Tehrani et al. 2019; Martins & Palacios 2022) because of their proximity to the Eddington limit due to their extreme luminosity (Vink et al. 2011; Bestenlehner et al. 2014; Gräfener 2021). For this reason, for stars with large hydrogen mass fraction at surface we set the transition from thin to thick winds when the Eddington factor Γ_e reaches some Γ_e,trans. Studies analyzing the wind efficiency have found a smooth transition between thin and thick winds, without abrupt jumps in $\dot{M}$ (Vink & Gräfener 2012; Sabhahit et al. 2023), but using enhanced values for the mass loss in the low Γ_e regime. For this work, we set Γ_e,trans = 0.5, based on Γ_e,trans ≃ 0.47 found by Bestenlehner et al. (2020).

For optically thin winds we adopt the mass-loss rate from Gormaz-Matamala et al. (2022b), based on their self-consistent m-CAK prescription (Gormaz-Matamala et al. 2019, 2022a).

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{\dot{M}}_{\mathrm{GM}23} & = & -40.314+15.438\,w+45.838\,x-8.284\,w\,x\\ & & +\,1.0564\,y-w\,y/2.36-1.1967\,x\,y\\ & & +\,z\times \left[0.4+\displaystyle \frac{15.75}{{M}_{* }/{M}_{\odot }}\right],\end{array}\end{eqnarray} \tag{ 1 }$$

where ${\dot{M}}_{\mathrm{GM}23}$ is in M_⊙ yr⁻¹; and where w, x, y, and z are defined as

$$\begin{eqnarray*}&&w=\mathrm{log}\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{\mathrm{kK}}\right),x=\displaystyle \frac{1}{\mathrm{log}g},y=\displaystyle \frac{{R}_{* }}{{R}_{\odot }},z=\mathrm{log}\left(\displaystyle \frac{{Z}_{* }}{{Z}_{\odot }}\right).\end{eqnarray*}$$

The range of validity of m-CAK prescription is $\mathrm{log}g$ ≥3.2 and T_eff ≥ 30 kK (Gormaz-Matamala et al. 2022a). However, if we extend ${\dot{M}}_{\mathrm{GM}23}$ for compact and cooler temperatures we can easily match with other mass-loss recipes such as Krtička & Kubát (2017, 2018) at cooler stars, the reason why we can revise the limit of validity as $\mathrm{log}g\geqslant 3.0$. Hence, ${\dot{M}}_{\mathrm{GM}23}$ covers the most part of the main sequence (MS) and a fraction of post-MS expansion. For the cases of $\mathrm{log}g\lt 3.0$, we keep the classical formula of Vink et al. (2000, 2001, hereafter Vink's formula), thus keeping the bistability jump in the mass loss around 25 kK (Vink et al. 1999) only for evolved stars. The existence of a jump in the regime of blue supergiants has been supported by recent theoretical Krtička et al. (2021, 2024) and empirical (Bernini-Peron et al. 2023) mass-loss studies.

For optically thick winds, i.e., when X_surf < 0.3 or Γ_e ≥ 0.5, we use the $\dot{M}\propto {{\rm{\Gamma }}}_{{\rm{e}}}$ relationship calculated by Bestenlehner (2020) for WNh stars.

$$\begin{eqnarray}&&\mathrm{log}{\dot{M}}_{{\rm{R}}136}=-5.19+2.69\mathrm{log}({{\rm{\Gamma }}}_{{\rm{e}}})-3.19\mathrm{log}(1-{{\rm{\Gamma }}}_{{\rm{e}}}).\end{eqnarray} \tag{ 2 }$$

The constants stem from the calibration done over R136 by Brands et al. (2022), which is located in the LMC (Z = 0.006) and thus is valid for that metallicity only. Therefore, to make Equation (2) applicable for different metallicities, we include an extra term to readapt Bestenlehner's formula for Z = 0.014.

$$\begin{eqnarray}&&\mathrm{log}{\dot{M}}_{{\rm{B}}20}=\mathrm{log}{\dot{M}}_{{\rm{R}}136}+\mathrm{log}\left(\displaystyle \frac{{Z}_{* }}{{Z}_{\mathrm{LMC}}}\right)\times \left[0.4+\displaystyle \frac{15.75}{{M}_{* }/{M}_{\odot }}\right].\end{eqnarray} \tag{ 3 }$$

The extra metallicity dependent is the same as the term introduced in Equation (1), normalized for Z_LMC. This formula is valid only for T_eff ≥ 30 kK, the reason why for cooler stars we kept the validity of Vink's formula.

For H-poor regime (X_surf ≤ 10⁻⁷), we adopt the formula from Sander & Vink (2020), based on the hydrodynamically consistent wind calculations for classical W-R stars Sander et al. (2020).

$$\begin{eqnarray}&&\mathrm{log}\dot{M}=a\mathrm{log}[-\mathrm{log}(1-{{\rm{\Gamma }}}_{{\rm{e}}})]-\mathrm{log}2{\left(\displaystyle \frac{{{\rm{\Gamma }}}_{{\rm{e}},{\rm{d}}}}{{{\rm{\Gamma }}}_{{\rm{e}}}}\right)}^{{c}_{d,b}}+\mathrm{log}{\dot{M}}_{\mathrm{off}},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray*}&&a=2.932,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{{\rm{\Gamma }}}_{{\rm{e}},{\rm{b}}}=-0.324* \mathrm{log}\left(\displaystyle \frac{{Z}_{* }}{{Z}_{\odot }}\right)+0.244,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{c}_{d,b}=-0.44* \mathrm{log}\left(\displaystyle \frac{{Z}_{* }}{{Z}_{\odot }}\right)+9.15,\end{eqnarray*}$$

$$\begin{eqnarray*}&&\mathrm{log}{\dot{M}}_{\mathrm{off}}=0.23* \mathrm{log}\left(\displaystyle \frac{{Z}_{* }}{{Z}_{\odot }}\right)-2.61.\end{eqnarray*}$$

Hydrodynamically consistent solutions from Sander & Vink (2020) have been recently upgraded by adding a temperature-dependent term (Sander et al. 2023). That term, however, corresponds to the temperature calculated for a stellar radius, which is not the same radius for W-R stars adopted by Genec (Meynet & Maeder 2005). Therefore, we adopt Equation (4) without the temperature dependence, because its incorporation for stellar evolution models is beyond the scope of this work. Additionally, for the sake of consistency we also did not adopt the temperature dependence in our MESA setup. We nevertheless acknowledge that the addition of this temperature-dependent term represents an improvement for the modeling of H-poor W-R stars, and future publications with the presented MESA models will incorporate it in our wind-driven mass-loss prescription.

2.2. Evolutionary Codes

The implementation of Γ_e as a complementary condition for the initiation of thick winds translates into increased mass-loss rates for stars above 60 M_⊙ at ZAMS. As noticeable in Figure 2, thick winds tend to be initiated comparatively earlier in the evolution of a star the higher its M_ZAMS is. The most massive stars almost spend no time at all in their thin winds phase, which leads them to promptly become WNh stars without any prior giant phase. This could potentially explain the lack of red supergiants (RSG) beyond the Humphrey–Davidson limit (Humphreys & Davidson 1994). This comes in agreement with the conclusions of Mennekens & Vanbeveren (2014), who reached a similar result through the alteration of the mass-loss rates during the luminous blue variable (LBV) phase. We nevertheless highlight that the difference between our results and the ones we cite arises from the fact that our results attribute the WNh winds mass loss as the prohibitor of massive RSGs to being observed beyond the Humphrey–Davidson limit, while Mennekens & Vanbeveren (2014) suggested it was due to the combined effect of wind mass loss and LBV outbursts.

Zoom In Zoom Out Reset image size

The earlier transition from thin to thick winds at Γ_e = 0.5 (instead of the previous X_surf = 0.3 or X_surf = 0.4) makes our new evolutionary tracks to find stars with optically thick envelopes and still a large fraction of hydrogen abundance, as the WNh stars observed in the Arches cluster (Martins et al. 2008; Martins 2023). We also show that the proposed mass-loss prescription leads every star at M_ZAMS ≳ 20 M_⊙ (roughly the minimum initial mass of BH progenitors at Z = 0.014 from single star evolution) to evolve into a W-R star during the latest stages of its evolution. Below that initial mass, the mass-loss rates and the internal mixing are not enough for stars to enter the thick winds regime. We highlight that this scenario may change depending on the adopted prescription for dust-driven winds (e.g., van Loon et al. 2005; González-Torà et al. 2024, 2023b). For instance Beasor et al. (2021, 2023) show considerably lower mass-loss rates than de Jager et al. (1988) and more specifically Beasor et al. (2023) claims that winds under this regime "are inefficient at removing the H envelope."

We report the final masses of stars as a function of their initial mass at solar metallicity calculated with MESA and Genec in Figure 3. We have unified the input key physics in both codes in terms of rotation, convection, overshooting, and wind-mass-loss rates. We did not calibrate the codes any further for the results to match. Yet, for the two test models (at M_ZAMS = 60 and 200 M_⊙) both codes produce almost identical final star masses that show robustness of our predictions. We also show as a comparison the estimates from Martinet et al. (2023) and an approximation of the BH mass distribution from Bavera et al. (2023). Most of the stars in our models, according to Fryer et al. (2012), collapse directly into BHs without any supernova event.

Zoom In Zoom Out Reset image size

In contrast with other models like Bavera et al. (2023) we do not report a quasi-flat distribution of BH masses for the most massive stars in our models, nor a BH mass peak in the low-mass range. On the contrary, the BH mass increases monotonically as a function of ZAMS mass to M_BH = 27.7 M_⊙ for our most massive model at M_ZAMS = 300 M_⊙. On the other hand, our models predict lower final star masses than the ones from Martinet et al. (2023) for the highest ZAMS masses. This is due to the fact that in canonical Genec simulations the switch from thin to thick winds is based only on H surface abundance, plus the upgrade of the mass-loss prescription for the late evolutionary stages (Higgins et al. 2021). This makes even the most massive stars spend a significant part of the lifetime in a thin (low) wind-mass-loss regime. In our models very massive stars enter the thick (high) wind-mass-loss regime very early in the evolution (note that we also use the Eddington factor in addition to surface abundance for the switch). As pointed out in Section 2.1 the upgraded switch using both H surface abundance and the Eddington factor seems in better agreement with observations.

Our study shows that in solar metallicity environments, assuming the current defaults for dust-driven mass loss, most of the BH progenitors in their latest stages of their evolution naturally become W-R stars even without invoking tidally induced rotational mixing or a mass transfer event with a stellar companion. This is due to the combined effect of internal mixing and strong stellar winds at high metallicity levels. Our results substantially differ from recent results in the context of BH masses at solar metallicity. We find a different initial–final star mass relation for our updated stellar models. First, our models predict a maximum BH mass of ∼28 M_⊙ that is smaller by at least ∼7 M_⊙ if compared to Bavera et al. (2023) and Martinet et al. (2023). Second, our most massive BHs are formed from significantly different initial star masses (M_ZAMS > 250 M_⊙) than the ones found by Bavera et al. (2023) from single star evolution (M_ZAMS ∼ 75 M_⊙). These differences are understood in terms of the updated wind-mass-loss prescriptions and of the internal mixing prescription (see Appendix B for more details on the role of overshooting and rotational velocity). Note that we have not only adopted different wind-mass-loss prescriptions but we have also proposed transition criteria among prescriptions that are more complex than traditionally adopted in the modeling of massive stars. Following these results we claim that the formation of a BH with a mass above 27 M_⊙ can likely only come from (a) interaction with other objects, be it accretion from interstellar medium or a mass transfer event with a star or merger with other stellar/compact objects, (b) the Big Bang (i.e., primordial BHs), or (c) the isolated evolution of a star with M_ZAMS ≳ 250 M_⊙.

The shape of the initial–final star mass relation at solar metallicity is an important calibration point in the predictions of the overall BH mass distribution in the Universe. LVK observatories are sensitive to detecting BHs from a wide range of metallicities (i.e., already reaching galaxies to redshift of z ∼ 1). Such studies as presented here need to be extended to a full range of metallicities before any reliable initial mass conclusions can be made about the BH mass distribution observed in gravitational waves. For lower metallicities wind-mass-loss rates are smaller and final star masses are larger than predicted for Z = 0.014. This will introduce another key factor/uncertainty in BH mass calculation: pair-instability pulsations mass loss and stellar disruptions (e.g., Farmer et al. 2020).

We are not claiming that our calculations represent a definitive answer to the question of BH mass distribution at solar metallicity. Yet, we argue that given the employed recent advancements in massive stellar evolution, this is an incremental step forward in understanding stellar-origin BHs.

Authors acknowledge support from the Polish National Science Center grant Maestro (2018/30/A/ST9/00050). We would also like to thank number of people for their helpful and/or critical comments on our modeling: Raphael Hirschi, Sylvia Ekström, Georges Meynet, Chris Fryer, Jakub Klencki, Simone Bavera, Tassos Fragos, Joachim Bestenlehner, and Gemma Gonzalez-Tora. The authors also thank the MESA community at large for their helpful feedback on the models creation. Computations for this article have been performed using the computer cluster at CAMK PAN. We dedicate this paper to Krzysztof Belczyński, who contributed to this research before his untimely passing on 2024 January 13.

The most remarkable difference between old and new wind prescriptions is the value of the mass-loss rate. These lower values of $\dot{M}$ show a better agreement with the current understanding of stellar winds, which are known to be clumped (Bouret et al. 2005; Šurlan et al. 2012, 2013). More detailed theoretical studies on the wind structure of O-type stars, self-consistently coupling radiative acceleration and hydrodynamics in the comoving frame and non–local thermodynamic equilibrium conditions (Sander et al. 2017; Gormaz-Matamala et al. 2021), have confirmed the same diagnostics.

For thin winds, we observe in Figure 4 that the self-consistent $\dot{M}$ from the m-CAK prescription from Gormaz-Matamala et al. (2023) agrees very well with the formula from Krtička & Kubát (2017, 2018) especially for MS stars as outlined by Gormaz-Matamala et al. (2022a). Wind solution from Björklund et al. (2021) also is lower than Vink's formula, although their $\dot{M}$ is extremely low at low luminosities, together with not evidencing the increase of $\dot{M}$ around ∼25 kK due to the bistability jump. However, for thick winds, we observe in Figure 4 that both Bestenlehner (2020) and Sander & Vink (2020) provide mass-loss rates below the former Nugis & Lamers (2000). Mass loss for H-rich stars is lower, because it describes the winds of stars in an intermediate evolutionary stage between O-type and W-R stars.

Zoom In Zoom Out Reset image size

However, notice that both panels of Figure 4 only show the relationship $\dot{M}\sim {L}_{* }$, and straightforward differences in mass-loss rates at the same luminosity may be misleading. For example, lowering mass-loss rates makes stars reach higher luminosities during their evolution (Gormaz-Matamala et al. 2022b; Björklund et al. 2023; Gormaz-Matamala et al. 2023), and therefore the mass-loss rates at subsequent evolutionary stages may be larger rather than smaller after adopting updated "lower" wind-mass-loss prescriptions. As another example, $\dot{M}$ for H-rich W-R stars (Bestenlehner et al. 2020) is lower than $\dot{M}$ for H-poor stars (Sander & Vink 2020) at the same luminosity, because they describe stars with a different L/M ratio (and thus different mass regimes).

We show in Figure 5 our parameter study for different initial rotation velocities (while keeping α_ov = 0.5) and convective overshooting parameter (with the initial Ω/Ω_crit = 0.4). Within the considered parameter space a more efficient mixing causes stars to enter their thick winds and/or dust-driven winds phases earlier in their evolution, which nullifies any noticeable difference in terms of helium core mass at the end of the MS. This explanation gains also a stronger physical motivation for very massive stars M_ZAMS ≳ 100 M_⊙, which are almost fully convective (i.e., they do not need overshooting nor rotation to mix elements throughout their whole volume) and enter the thick winds regime almost as early as ZAMS. As a reference, we also show some extreme cases on the right-hand plot of Figure 5 with very inefficient mixing within stellar interiors (no rotation and α_ov ≤ 0.2), if compared with our reference model described in Section 2. We report that for nonrotating stars there is no smooth transition as a function of α_ov from the M_final behavior with a local maximum at M_ZAMS ∼ 60 M_⊙ shown for α_ov ≤ 0.15 to the monotonic behavior that is noticeable in all other models. We explain that this drastic transition is due to the fact that for these models with relatively low mixing efficiencies stars do not reach the W-R phase. To support this statement we show Table 1, where we report the cumulative mass lost from the adopted wind-mass-loss prescriptions for a 60 M_⊙ star (roughly the M_ZAMS where we see the localized M_final maximum for the lower mixing efficiency cases) for different initial conditions of α_ov and Ω/Ω_crit.

Zoom In Zoom Out Reset image size