The Formation of a 70 M_⊙ Black Hole at High Metallicity

LB-1 is reported as a detached binary system containing a B star with a mass of 8 ${M}_{\odot }$ (−1.2/ + 0.9 ${M}_{\odot }$) and a black hole (BH) with a mass of 68 ${M}_{\odot }$ (−13/ + 11 ${M}_{\odot }$). The binary system orbit is almost circular, with e = 0.03 (−0.01/ + 0.01 ${M}_{\odot }$), and has an orbital period of P_orb = 78.9 days (−0.3/ + 0.3 days). This corresponds to a physical semimajor axis of a = 300–350 ${R}_{\odot }$ and a Roche lobe radius of the BH R_BH,lobe ≲ 200 ${R}_{\odot }$. This system is one of the widest known binary system hosting a stellar-origin BH; see https://stellarcollapse.org. Two other binaries, proposed to host BH candidates, were also discovered by the radial velocity method by Thompson et al. (2019) and Giesers et al. (2018; although this is in a globular cluster and has a very large period, P = 167 days, an eccentric orbit, e = 0.6, and it must have formed by capture).

LB-1 was discovered by the 4 m class telescope LAMOST and the spectroscopic orbit was confirmed by the 10 m class Gran Telescopio Canarias and Keck telescopes. Chandra non-detection places the X-ray emission at the very low level, <2 × 10³¹ erg s⁻¹. An H_α emission line was observed, however, and since it follows a BH (small accretion disk around the BH from the B star wind) the double spectroscopic orbital solution was obtained. The system is on the outskirts of the Galactic disk, in the anti-Galactic center direction, about 4 kpc away from Sun. There is no globular cluster nearby (<4 kpc). The chemical composition of the B star indicates a metal abundance Z = 0.02 that is slightly over solar, assuming Z_⊙ = 0.017. The full information on the system parameters and the discovery is reported in Liu et al. (2019).

Since the publication of the discovery paper, there are a number of studies that attempt either to reject specific formation scenarios of LB-1 (the massive BH is the BH–BH merger product or a very close BH–BH binary; see Shen et al. 2019) or to explain it with some specific scenarios: stellar evolution of a massive magnetic star (Groh et al. 2019), merger of two unevolved stars (Tanikawa et al. 2019), merger of a BH and an unevolved star (Banerjee 2019; Olejak et al. 2019). It was also pointed out that the existence of LB-1 (and its future evolution) may be in tension with the non-detection of ultraluminous X-ray sources or BH neutron star systems in the Galaxy (Safarzadeh et al. 2019). Alternatively, the nature of LB-1 is questioned with a reanalysis of observational data and results that support the idea that either the BH or both components are of lower mass than originally claimed (Abdul-Masih et al. 2019; El-Badry & Quataert 2019; Eldridge et al. 2019; Irrgang et al. 2019; Simón-Díaz et al. 2019). This would allow the classical scenario of an isolated binary evolution at high metallicity to explain the formation of LB-1.

In fact, the existence of a 70 ${M}_{\odot }$ BH in a high-metallicity environment seems challenging. The current consensus is based on mass-loss rate estimates and their dependence on metallicity for H-rich stars (Vink et al. 2001) and He-rich stars (Vink & de Koter 2005; Sander et al. 2020) that seem to limit BH mass to about 20 ${M}_{\odot }$ at solar metallicity (Belczynski et al. 2010). Existing electromagnetic observations seem to support this paradigm (Casares & Jonker 2014). Note the masses of the two most massive stellar-origin BHs that are known to have formed at relatively high metallicity are the well known Cyg X-1 (M_BH = 14.8 ± 1.0 Z ≈ 0.02; Orosz et al. 2011) and M33 X-7 (M_BH = 15.7 ± 1.5, Z ≈ 0.1 Z_⊙; Valsecchi et al. 2010).

The mass of the BH in LB-1 seems to contradict pair-instability pulsation supernovae (PPSNe) and pair-instability supernova (PSN) theory, which limit BH mass to about M_BH < 40–50 ${M}_{\odot }$ (Bond et al. 1984; Heger & Woosley 2002; Woosley 2017; Farmer et al. 2019; Leung et al. 2019). This limit was recently proposed to be as high as ∼55 ${M}_{\odot }$ for non-zero metallicity stars (Population I/II) by Belczynski et al. (2017). Note that for the LIGO/Virgo most massive BH–BH merger in O1/O2 (GW170729), the primary BH mass was reported to be 51.2 ${M}_{\odot }$. This high mass (not the merger itself) is likely a statistical fluctuation (Fishbach et al. 2019). However, even that mass can be explained as long as the BH was formed at low metallicity. The PPSN/PSN instability can be avoided (at best) for a 70 ${M}_{\odot }$ star that produces as 69 ${M}_{\odot }$ BH if only small neutrino mass loss takes place at the BH formation. This was envisioned to be plausible for an ultra-low metallicity and Population III stars, as they can keep massive H-rich envelopes (Heger & Woosley 2002; Woosley 2017).

Here, we propose that a similar mechanism may also work at high metallicity. The modification that we need to introduce to stellar evolution is to lower wind mass-loss rates for (at least some) massive stars. The empirical diagnostics of the winds of massive stars are complex, particularly because of the wind clumping (Fullerton et al. 2006; Oskinova et al. 2007), and the agreement between theory and observations is not always conclusive (Keszthelyi et al. 2017).

In lower-metallicity environments, such as in the LMC and the SMC galaxies, some studies have (Massa et al. 2017) indicated that wind mass-loss rates may actually be higher than the values typically adopted in evolutionary predictions (Vink et al. 2001; Belczynski et al. 2010). Other studies, however, seem to agree with standard calculations (Ramírez-Agudelo et al. 2017), and yet others point to much lower mass-loss rates than expected (Bouret et al. 2003; Ramachandran et al. 2019; Sundqvist et al. 2019). In the upper stellar mass regime, Hainich et al. (2013, 2019) determined mass-loss rates that are in broad agreement with the theoretical expectations.

In this work we consider the mass regime 70–100 ${M}_{\odot }$ at solar metallicity. Vink & Gräfener (2012) have shown that for stars in the transitional regime (from optically thin to thick winds) the standard mass-loss rates should apply. The empirical studies that include hydrogen-rich Wolf–Rayet stars (Hamann et al. 2019) find mass-loss rates that are lower than the values theoretically predicted by Nugis & Lamers (2000) for the most luminous objects. However, what the mass-loss rates of such massive stars are when they are very young is not well known. Gruner et al. (2019) found that the mass-loss rates of the earliest O-type star in the Galaxy (HD 93129A, the primary mass is ∼100 ${M}_{\odot }$) compares well with the theoretical expectations, but this result depends on an assumed clumping parameters. Furthermore, about 7% of OB stars are known to have (mostly) dipolar magnetic fields (Fossati et al. 2015; Wade et al. 2016; Grunhut et al. 2017). Some of these known magnetic stars are massive, but do not quite reach the mass regime considered here unless errors on mass estimates are considered: 61 ± 33 ${M}_{\odot }$ for CPD-28 2561 or ∼60 ${M}_{\odot }$ for HD 148937 (David-Uraz et al. 2019). These magnetic fields may capture wind particles and reduce wind mass-loss rates independent of star mass and metallicity (Owocki et al. 2016; Georgy et al. 2017; Petit et al. 2017; Shenar et al. 2017). Here we show two things: (1) the decrease of wind mass-loss rates (independent of the reduction origin) allows some models to avoid pair-instability associated mass loss and allows for the formation of high-mass BHs (∼50–70 ${M}_{\odot }$) at high metallicity; and (2) we are not able to make such a massive BH progenitor star fit within the binary orbit of LB-1, if in fact LB-1 hosts a 70 ${M}_{\odot }$ BH.

2.1. Simple StarTrack Simulation

We used the population synthesis code StarTrack (Belczynski et al. 2002, 2008) to quickly test the possibility of the formation of a 70 ${M}_{\odot }$ BH with decreased wind mass loss. We employed the rapid core-collapse supernova (SN) engine NS/BH mass calculation (Fryer et al. 2012), with strong PPSN/PSN mass loss (Belczynski et al. 2016). Standard winds for massive stars are used as the base model: O/B star (Vink et al. 2001) winds and LBV winds (specific prescriptions for these winds are listed in Section 2.2 of Belczynski et al. 2010). In wind mass-loss prescriptions we introduce a multiplication factor of f_wind = 1.0 as our standard calculation. Note that this approach produces a maximum of ∼15 ${M}_{\odot }$ for BHs at high metallicity (Z = 0.02 assuming Z_⊙ = 0.017), as demonstrated in Figure 1. We also calculate the evolution of single stars for decreased winds for two extra models with f_wind = 0.5, 0.2. It is clear from Figure 1 that winds need to be reduced by a factor of ∼5 to produce a ∼70 ${M}_{\odot }$ BH at high metallicity.

Zoom In Zoom Out Reset image size

Our specific example is a star with M_zams = 104 ${M}_{\odot }$ at Z = 0.02 and the star is evolved with Hurley et al. (2000) analytic formulae (used in many population synthesis and globular cluster evolutionary codes). H-rich wind mass-loss rates are decreased with f_wind = 0.2. The star keeps its H-rich envelope throughout the entire evolution. After 3.8 Myr of evolution, the star has a mass of M_tot = 69.8 ${M}_{\odot }$ with an H-rich envelope mass of M_env = 24.8 ${M}_{\odot }$, He core mass of M_He = 44.99 ${M}_{\odot }$, and CO core mass of M_CO = 34.8 ${M}_{\odot }$. According to the simplistic population synthesis prescription (no PPSN/PSN for stars with M_He < 45.0 ${M}_{\odot }$; Woosley 2017) this star is not yet subject to PPSN/PSN. The star undergoes core collapse and with 1% neutrino mass loss it forms a BH through direct collapse: M_BH = 69.1 ${M}_{\odot }$.

2.2. Single-star Evolutionary Models

To explore the possibility of the LB-1 BH being the descendant of a single star and to test simple estimates from Section 2.1, we ran a series of stellar evolution models using the MESA code revision 11701 (Paxton et al. 2011; 2013; 2015; 2018; 2019). We used a solar initial composition of Z = 0.014 for all models with an Asplund et al. (2009) metal mixture (initial_zfracs = 6), and the corresponding opacity tables (kappa_file_prefix=''a09'') including low-temperature tables (kappa_lowT_prefix = ''lowT_fa05_a09p'') and C/O-enhanced (type 2) opacity tables (kappa_CO_prefix = ''a09_co''). For convection, we used the Schwarzschild boundary location condition and included convective boundary mixing with a value of the exponentially decaying diffusion coefficient parameter f and f₀ everywhere equal to 0.004. For the reaction network, we used the basic.net and auto_extend_net = .true., with which MESA adapts the network along the evolution to the smallest network needed to trace energy generation. The main "stabilizing" setting/approximation was the use of extra pressure at the surface of the star by setting Pextra_factor=2. Another one was the use of MLT++ (see Paxton et al. 2013, Section 7). These settings might underestimate the radius of the star in our models. The models were evolved at least until the end of core He burning and generally stopped due to convergence issues near the end of core carbon burning.

We used the "Dutch" scheme for mass loss with a default Dutch_scaling_factor = 0.85. The two main mass-loss prescriptions experienced by our hydrogen-rich models are from Vink et al. (2001) for hot stars and from de Jager et al. (1988), which we used for the cool "Dutch" wind. In order to reduce the mass-loss rates, we lowered the Dutch_scaling_factor by introducing a multiplication factor in front of wind mass-loss rates and changing it over a wide range, f_wind = 1.0–0.0. We calculated nonrotating and rotating models (see Table 1). The standard rotation settings were used (Heger et al. 2000). Rotation is set on the zero-age main-sequence and the initial rotation rate, in terms of Ω/Ω_crit, is given in Table 1. We include the following rotation-induced instabilities; Eddington–Sweet circulation, secular shear instability, and Taylor–Spruit dynamo (Spruit 2002). Table 1 gives the key properties of representative stellar models. Using the physical ingredients described above and considering that the main uncertainty in the models is mass loss, we reduced the mass loss with a multiplication factor given in the table in an attempt to produce a final total mass equal to that of LB-1, i.e., around 70 ${M}_{\odot }$.

Table 1.  Initial Mass, Rotation, and Mass-loss Rescaling Factor (Columns 1–3) and Final Total, He and CO Core Masses, and Maximum Radius (Columns 4–7) of the Stellar Models

Mzams	Ω/Ωcrit	fwind	Mtot	MHe	MCO	Rmax/ ${R}_{\odot }$
Nonrotating models
100	0.0	0.576	70.8	41.5	36.9	711.1
85	0.0	0.333	70.9	31.6	27.6	653.9
70	0.0	0.0	70.0	30.8	27.0	637.5
Rotating models
100	0.6	0.576	61.6	49.5	43.9	260.8
85	0.6	0.576	58.2	40.3	35.4	363.9
85	0.6	0.333	62.9	46.8	41.3	235.0
75	0.6	0.576	53.9	34.5	30.1	376.5
70	0.6	0.576	50.2	32.1	27.8	324.1
70	0.4	0.282	58.5	32.5	28.3	611.8
Rotating models losing entire H-layer
100	0.8	1.0	40.5	40.5	36.8	170.9
100	0.8	0.882	43.4	43.4	37.5	165.5

Download table as:  ASCIITypeset image

Considering nonrotating models, a model without mass loss (M_zams = 70 ${M}_{\odot }$, f_wind = 0.0) is also included for reference as the most extreme (and unrealistic) case. With the rescaled wind, f_wind = 0.576, a model with an initial mass of 100 ${M}_{\odot }$ ends up having a total mass 70.8 ${M}_{\odot }$. This model has final core masses that will experience pair-instability pulsation mass loss, however, and thus lose more mass before it produces a BH. Furthermore, its radius is too large to fit in the orbit of LB-1. The most interesting model is M_zams = 85 ${M}_{\odot }$, with f_wind = 0.333. The final total mass is 70.9 ${M}_{\odot }$ and very importantly the final CO core mass is below the limit for PPSN mass loss. Indeed, the CO core mass of this model is M_CO = 27.6 ${M}_{\odot }$ (see Figure 2), which is below the CO core mass threshold for PPSN according to Table 1 in Woosley (2017; no pulsations for models with CO core masses below 28 ${M}_{\odot }$). It is thus possible for this model to produce a 70 ${M}_{\odot }$ BH. Unfortunately, the maximum radius of this model (R_max ≈ 650 ${R}_{\odot }$; see Figure 3) is too large to fit in the Roche lobe of the LB-1's BH (<200 ${R}_{\odot }$) and this model thus cannot provide a full solution for the origin of LB-1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Considering rotating models, rotation-induced mixing leads to more massive cores and more mass loss (e.g., Hirschi et al. 2004). Thus, the rotating models with similar initial parameters end up with smaller total masses and larger core masses, which makes them less suitable candidates to explain LB-1. The only advantage of rotating models over nonrotating ones is that they end with smaller radii that could possibly fit in the LB-1. So we may then ask the question: what is the most massive final single-star model that would always fit in LB-1? Considering models that lose entire H-rich layers (e.g., 100 ${M}_{\odot }$ models at the bottom of Table 1) or pure He-star models of Woosley (2017) or Farmer et al. (2019), BH masses up to 45–50 ${M}_{\odot }$ can be produced. Since He-stars are very compact, these would fit within the LB-1 BH Roche lobe but the BH mass would be below the current lower-mass limit of 55 ${M}_{\odot }$ for LB-1. We also consider rotating models with f_wind = 0.576 that do not lose H-rich layers: M_zams = 70, 75, 85, 100 ${M}_{\odot }$. The 70 ${M}_{\odot }$ model has a final CO core mass below the pair-instability pulsation mass range and thus is likely to collapse to a BH with little mass loss. The final radius, however, is not so small (R_max = 324 ${R}_{\odot }$) and that model would not fit in the Roche lobe of LB-1's primary; the total mass is smaller than the lower-mass estimate of BH mass in LB-1. The 75 and 85 ${M}_{\odot }$ have larger final masses but also larger CO core masses and radii and thus will likely lose some mass by pair-instability supernova pulsations and thus not fit in LB-1. The 100 ${M}_{\odot }$ model produces a relatively small maximum radius (R_max = 261 ${R}_{\odot }$), but it would still not fit in the Roche lobe of LB-1's primary. Although the final model mass is large (above the lower limit on LB-1 BH mass), this model has a massive CO core and is subject to strong PPSN mass loss. A similar case is found for the 85 ${M}_{\odot }$ rotating model with f_wind = 0.333. It thus seems very unlikely that a single star or a noninteracting star in a binary system could produce the BH in LB-1 with the currently derived properties.

It is generally believed that Population I/II stars cannot form BHs in the mass range ∼55–135 ${M}_{\odot }$, the so-called second mass gap, due to mass loss in PPSNe and due to total star disruption by pair-instability supernovae. It has been noted, however, that in one specific case the lower bound of the second mass gap can be shifted to ∼70 ${M}_{\odot }$. Such a case was proposed for metal-poor (Population III) stars, for which wind mass loss is negligible even for high-mass stars and then such stars can retain massive H-rich envelopes throughout their evolution. The retention of a massive H-rich envelope allows a star to ignite an H-burning shell, which supports the outer stellar layers and helps the density/temperature in the stellar interior avoid the pair-instability regime (where the adiabatic index becomes small γ < 4/3). In principle, one can imagine a stable (against PPSN/PSN) stellar configuration with a 70 ${M}_{\odot }$ star at the core collapse, with an He core mass of ≲40 ${M}_{\odot }$ and H-rich envelope of ≳30 ${M}_{\odot }$ for a metal-poor star (for which mass loss is expected to be low, at least lower than that at high metallicity).

We found that a similar configuration can be achieved for high-metallicity stars if wind mass-loss rates are decreased in stellar evolution models. For one model, a nonrotating M_zams = 85 ${M}_{\odot }$ and Z = 0.014 star, we can form a 70 ${M}_{\odot }$ BH as a single-star or a binary component in a very wide noninteracting binary if standard wind mass-loss rates are reduced by factor of ∼5. This is a rather surprising and unexpected result on its own. Note that this result is totally independent of LB-1 and its true nature, whether it hosts a massive BH or not. This model, however, is not useful in the context of LB-1, as the stellar radius of this star (≳650 ${R}_{\odot }$) is too large to fit within LB-1's orbit.

The main uncertainty in the massive star models is mass loss. We reduced the mass-loss rates in order to produce higher final masses. Note that reduced wind mass loss does not have to be in effect for all stars, but it may be possible that wind is quenched only for some fraction of very massive stars (e.g., via magnetic capture of wind particles: see Section 1). Other studies (e.g., Limongi & Chieffi 2018; Chieffi & Limongi 2019), however, show that a higher mass loss is needed in the red supergiant (RSG) phase to reproduce the absence of observed SNe II above a certain luminosity (Smartt 2009). Evolved massive stars are also expected to lose mass via eruptive events, e.g., LBV-type mass loss, beyond the Humphreys–Davidson limit (Humphreys & Davidson 1979; Langer 2012; Smith 2014). This extra mass loss was suggested to explain the apparent lack of cool luminous massive stars in the Milky Way (Mennekens & Vanbeveren 2014). Note that our model of a nonrotating 85 ${M}_{\odot }$ star that can produce a 70 ${M}_{\odot }$ BH enters the cool (log₁₀(T_eff) ≈ 3.9) and luminous (log₁₀(L/L_⊙) ≈ 6.3) regions of the H-R diagram (see Figure 3). Even at low metallicity, Small Magellanic Cloud stars are not found at such low temperatures and high luminosities (see Figure 13 of Ramachandran et al. 2019).

Therefore, the existence of LB-1, if it really hosts a massive 70 ${M}_{\odot }$ BH, points to some other possibilities. (i) Pair-instability does not operate in stars as expected. This would allow a rapidly rotating massive star to evolve homogeneously, retaining a small radius and forming a 70 ${M}_{\odot }$ helium-rich object that would directly collapse to a BH. (ii) Or the BH is a descendant of a BH–BH or BH–star merger in the inner binary and LB-1 was originally a triple system. Note that this would also require homogeneous evolution of two ∼30–50 ${M}_{\odot }$ stars to not affect a nearby B star, but this would not require violating pair-instability theory. However, a gravitational-wave kick during a BH–BH merger or any natal kick at BH formation might be incompatible with the very low eccentricity of LB-1. (iii) Perhaps some stars expand less due to an exotic composition and modifications of opacities or to an unknown additional mixing process. Alternatively, LB-1 may have lower-mass components than was claimed in the discovery paper and then standard stellar/binary evolution can account for the formation of such system.

Note that if BHs as massive as 70 ${M}_{\odot }$ exist in young and metal-rich environments, e.g., Galactic disk, they would most likely have low spins, since our models employ effective angular momentum transport by a magnetic dynamo (a ≲ 0.15; see Belczynski et al. 2017). If such a massive BH could catch a companion, e.g., in an open cluster, or have formed in a wide binary with another BH that then evolves into close/merging system, e.g., by a "lucky" natal kick injection into a short period and eccentric orbit, then LIGO/Virgo will sooner or later discover such massive BHs. LIGO/Virgo detection of objects of such mass will be burdened with large errors, ∼20–30 ${M}_{\odot }$ up and down, so in principle even a detection of a 100 ${M}_{\odot }$ BH could be possibly explained by one of our models.

The authors would like to thank Lida Oskinova and the anonymous reviewer for useful comments. K.B. acknowledges support from the Polish National Science Center (NCN) grant Maestro (2018/30/A/ST9/00050). K.B. and R.H. acknowledge support from the World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. R.H. and E.K. acknowledge support from the ChETEC COST Action (CA16117), supported by COST (European Cooperation in Science and Technology). A.H. acknowledges support from the National Science Foundation under grant No. PHY-1430152 (JINA Center for the Evolution of the Elements) and from the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. J.C. acknowledges support by the Spanish Ministry of Economy, Industry and Competitiveness (MINECO) under grant AYA2017-83216-P.