On the Maximum Black Hole Mass at Solar Metallicity

In high-metallicity environments the mass that black holes (BHs) can reach just after core collapse widely depends on how much mass their progenitor stars lose via winds. On one hand, new theoretical and observational insights suggest that early-stage winds should be weaker than what many canonical models prescribe. On the other hand, the proximity to the Eddington limit should affect the formation of optically thick envelopes already during the earliest stages of stars with initial masses M ZAMS ≳ 100 M ⊙, hence resulting in higher mass-loss rates during the main sequence. We use the evolutionary codes MESA and Genec to calculate a suite of tracks for massive stars at solar metallicity Z ⊙ = 0.014, which incorporate these changes in our wind-mass-loss prescription. In our calculations we employ moderate rotation, high overshooting, and magnetic angular momentum transport. We find a maximum BH mass MBH,max=28.3 M ⊙ at Z ⊙. The most massive BHs are predicted to form from stars with M ZAMS ≳ 250 M ⊙, with the BH mass directly proportional to its progenitor’s M ZAMS. We also find in our models that at Z ⊙ almost any BH progenitor naturally evolves into a Wolf–Rayet star due to the combined effect of internal mixing and wind mass loss. These results are considerably different from most recent studies regarding the final mass of stars before their collapse into BHs. While we acknowledge the inherent uncertainties in stellar evolution modeling, our study underscores the importance of employing the most up-to-date physics in BH mass predictions.


INTRODUCTION
The LIGO/Virgo/KAGRA (LVK) collaboration has delivered summary of detections of compact object mergers in gravitational-waves in 3 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 M ⊙ and ∼35 M ⊙ .Understanding such 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. 2010b,a).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 chemi- † Deceased on 13th January 2024.amedeoromagnolo@gmail.comcal 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), being this value the lower threshold for theoretical estimations of M 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 non-rotating 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).

MODELING
In our study we use two different detailed evolutionary codes, the version 23.05.1 of Modules for Experiments in Stellar Astrophysics (MESA Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 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.

Mass loss prescription
The theoretical calculation of wind mass loss rates ( Ṁ ) 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.
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 (WR) 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 Ṁ (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).
where ṀGM23 is in M ⊙ yr −1 ; and where w, x, y, and z are defined as The range of validity of m-CAK prescription is log g ≥ 3.2 and T eff ≥ 30 kK (Gormaz-Matamala et al. 2022a).However, if we extend ṀGM23 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, reason why we can revise the limit of validity as log g ≥ 3.0.Hence, ṀGM23 covers the most part of the Main Sequence (MS) and a fraction of post-MS expansion.For the cases of log g < 3.0, we keep the classical formula of Vink et al. (2000Vink et al. ( , 2001, hereafter Vink's formula), thus keeping the bi-stability 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 BSGs has been supported by recent theoretical Krtička et al. (2021Krtička et al. ( , 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 Ṁ ∝ Γ e relationship calculated by   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 Eq. 2 applicable for different metallicities, we include an extra term to readapt Bestenlehner's formula for Z = 0.014.
(3) The extra metallicity dependent is the same as the term introduced in Eq. 1, normalized for Z LMC .This formula is valid only for T eff ≥ 30 kK, reason why for cooler stars we kept the validity of Vink's formula.
For H-poor regime (X surf ≤ 10 −7 ), we adopt the formula from Sander & Vink (2020), based on the hydrodynamically consistent wind calculations for classical WR stars Sander et al. (2020).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 WR stars adopted by Genec (Meynet & Maeder 2005).Therefore, we adopt Eq. 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 didn't 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 WR stars, and future publications with the presented MESA models will incorporate it in our wind-driven mass loss prescription.

Evolutionary Codes
For MESA, we adopt the Ledoux criterion for convective boundaries (Ledoux 1947) and mixing length l = 1.82 (Choi et al. 2016).We use step-overshooting, with a value of α ov above every convective region that represents a modified version of the POSYDON prescription described in Fragos et al. (2023).For M ZAMS < 4 M ⊙ we use α ov = 0.16, which was adopted in the MIST project from calibrations of the Sun and open clusters (Choi et al. 2016).For 8 M ⊙ < M ZAMS < 20 M ⊙ we use α ov = 0.415 from the Brott et al. (2011) work.Between 4 M ⊙ and 8 M ⊙ at ZAMS we interpolate α ov between 0.16 and 0.415.Finally, for M ZAMS ≥ 20 M ⊙ we adopt α ov = 0.5 (Scott et al. 2021).Below each convective region we also adopt a value of undershooting α under that is α under = α ov /5.In order to reduce superadibaticity in regions near the Eddington limit we use the use superad reduction method, which is stated to be a more constrained and calibrated prescription than MLT++ in MESA (Jermyn et al. 2023).For each star we adopt an equatorial velocity over critical velocity at ZAMS equal to Ω/Ω crit = 0.41 .To model the rotation and mixing of material we use the calibrations from Heger et al. (2000), with the addition of Tayler-Spruit magnetic field-driven diffusion (Tayler 1973;Spruit 2002;Heger et al. 2005).
For Genec we also employ the Ledoux criterion as Georgy et al. (2014) and Sibony et al. (2023) for convective boundaries.The treatment of rotation and its respective impact on mass loss is given by Maeder & Meynet (2000), with initial rotation of our models: V rot /V crit = 0.4.The transport of angular momentum inside the star follows the prescription of Zahn (1992) complemented by Maeder & Zahn (1998), whereas here we use Taylor-Spruit dynamo following Maeder & Meynet (2004).Same as for MESA we also use α ov = 0.5, whereas the abundances correspond to the values given by Asplund et al. (2009), whereas opacities come from Iglesias & Rogers (1996).

RESULTS
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 winds mass loss and LBV outbursts.
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 Wolf-Rayet 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. 2023b,a).For instance Beasor et al. (2021Beasor et al. ( , 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 M ⊙ 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 distribu-tion from Bavera et al. (2023).Most of the stars in our models, according to Fryer et al. (2012), collapse directly into BHs without any SN event.
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 to spend significant part of the lifetime in 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.

CONCLUSION
Our study shows that in solar metallicity environments, assuming the current defaults for dust-driven mass loss, most of the black hole progenitors in their latest stages of their evolution naturally become Wolf-Rayet 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 black hole masses at solar metallicity.We find a different initial-final star mass relation for our updated stellar models.Firstly, our models predict a maximum black hole mass of ∼ 28 M ⊙ that is smaller by at least ∼ 7 M ⊙ if compared to Bavera et al. (2023) and Martinet et al. (2023).Secondly, our most massive black holes 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 black hole with a mass above 27 M ⊙ can likely only come from a) interaction with other objects, be it accretion from interstellar medium, mass transfer event with a star or merger with other stellar/compact objects b) the Big Bang (i.e.primordial black holes) 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 black hole mass distribution in Universe.LIGO/Virgo/Kagra observatories are sensitive to detecting black holes 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 black hole 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 black hole 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 black hole 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 black holes.
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 modelling: 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 the 13 th of January 2024.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 r.h.s.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 non-rotating 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 WR 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 .

Fig. 2 .
Fig. 2.-Left: H-R diagram for selected MESA models.Colors correspond to various wind mass loss prescriptions used in our study.Note that stars with initial masses M ZAMS ≲ 60 M ⊙ are subject to significant radial expansion, while more massive stars do not expand during evolution.Right: Fraction of ZAMS mass lost through given wind mass loss prescription (color coded as shown in the figure).Note that majority of mass loss for most massive stars (M ZAMS ≳ 60 M ⊙ ) is expected by H-rich thick winds and by H-poor thick winds (for which we use recent updated formulae).

Fig. 3 .
Fig. 3.-ZAMS mass -final star mass (end of C-burning) relation at solar metallicity Z = 0.014.Our updated MESA models are shown with blue points (final mass of a star equals mass of a BH: no supernova) and blue stars (upper limit on mass of a BH, where blue squares show our estimates of NS/BH mass after supernova).Our two updated Genec models are shown with red squares.Predictions from both codes match.We predict a monotonic growth of the BH mass as a function of its ZAMS mass.The predicted maximum BH mass in our models is M BH,max = 28.3M ⊙ and comes from a ZAMS mass M ZAMS ∼ 300 M ⊙ .As a comparison we also show the models from Bavera et al. (2023) for single star evolution (yellow line) and from Martinet et al. (2023) (green triangles): note that these models predict considerably different BH masses than our models.log Ṁoff = 0.23 * log Z * Z ⊙ − 2.61 .

Fig. 5 .
Fig. 5.-ZAMS mass -final star mass (end of core C-burning) relation at Z ⊙ .On the l.h.s plot we compare the results of our simulations between different initial rotation velocities at ZAMS.On the r.h.s.plot, instead, we compare the results from the simulations made with different values for the convective overshooting parameter αov.The dotted lines represent non-rotating stars with αov ≤0.2.