¹¹institutetext: Département d’Astronomie, Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland ²²institutetext: Gravitational Wave Science Center (GWSC), Université de Genève, CH1211 Geneva, Switzerland ³³institutetext: Institute of Astrophysics, FORTH, N. Plastira 100, Heraklion, 70013, Greece ⁴⁴institutetext: Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong ⁵⁵institutetext: School of Physics and Astronomy, Monash University, VIC 3800, Australia ⁶⁶institutetext: OzGrav: The Australian Research Council Centre of Excellence for Gravitational Wave Discovery, Clayton, VIC 3800, Australia ⁷⁷institutetext: Department of Physics, University of Florida, 2001 Museum Rd, Gainesville, FL 32611, USA ⁸⁸institutetext: Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA ⁹⁹institutetext: Department of Physics & Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA ¹⁰¹⁰institutetext: Electrical and Computer Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA ¹¹¹¹institutetext: Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Magrans, 08193 Barcelona, Spain ¹²¹²institutetext: Institut d’Estudis Espacials de Catalunya (IEEC), Carrer Gran Capità, 08034 Barcelona, Spain

Formation of Wind-Fed Black Hole High-mass X-ray Binaries:
The Role of Roche-lobe-Overflow Post Black-Hole Formation

Zepei Xing ^, e-mail: Zepei.Xing@unige.ch1122 Tassos Fragos 1122 Emmanouil Zapartas 33 Tom M. Kwan 44 Lixin Dai 44 Ilya Mandel 5566 Matthias U. Kruckow 1122 Max Briel 1122 Jeff J. Andrews 77 Simone S. Bavera 1122 Seth Gossage 88 Konstantinos Kovlakas 11111212 Kyle A. Rocha 8899 Meng Sun 88 Philipp M. Srivastava 881010

(Received XXYYZZ; Accepted XXYYZZ)

The three dynamically confirmed wind-fed black hole high-mass X-ray binaries (BH-HMXBs) are suggested to all contain a highly spinning black hole (BH). However, based on the theories of efficient angular momentum transport inside the stars, we expect that the first-born BHs in binary systems should have low spins, which is consistent with gravitational-wave observations. As a result, the origin of the high BH spins measured in wind-fed BH-HMXBs remains a mystery. In this paper, we conduct a binary population synthesis study on wind-fed BH-HMXBs at solar metallicity with the use of the newly developed code POSYDON, considering three scenarios for BH accretion: Eddington-limited, moderately super-Eddington, and fully conservative accretion. Taking into account the conditions for accretion-disk formation, we find that regardless of the accretion model, these systems are more likely to have already experienced a phase of Roche-lobe overflow after the BH formation. To account for the extreme BH spins, highly conservative accretion onto BHs is required, when assuming the accreted material carries the specific angular momentum at the innermost stable orbit. Besides, in our simulations we found that the systems with donor stars within the mass range of $10-20\,M_{\odot}$ are prevalent, posing a challenge in explaining simultaneously all observed properties of the BH-HMXB in our Galaxy, Cygnus X-1, and potentially hinting that the accretion efficiency onto non-degenerate stars, before the formation of the BH, is also more conservative than assumed in our simulations.

Key Words.:

X-rays: binaries — Stars: black holes — Stars: evolution

1 Introduction

Black hole X-ray binaries contain a black hole (BH) accreting material from a companion star, generating X-ray emission. Depending on the donor star mass and accretion processes, black hole X-ray binaries can be classified as low-mass X-ray binaries (BH-LMXBs) and high-mass X-ray binaries (BH-HMXBs). BH-LMXBs consist of a low-mass donor star ( $\sim 2\,M_{\odot}$ ) that transfers mass through Roche lobe overflow onto the BH. BH-HMXBs typically involves a massive ( $\gtrsim 10\,M_{\odot}$ ) OB-type donor star, with the BH accreting a fraction of the companion’s stellar wind. Only three wind-fed BH-HMXBs have been dynamically confirmed: Cygnus X-1, LMC X-1 and M33 X-7 (we list their properties in Table 1). All of them are observed to contain a BH with a high spin and a donor star that is still on its main sequence (MS), almost filling its Roche lobe. Besides, there are several candidates including SS433 (Seifina & Titarchuk 2010), Cygnus X-3 (Zdziarski et al. 2013), HD96670 (Gomez & Grindlay 2021), IC10 X-1 (Prestwich et al. 2007; Laycock et al. 2015), and NGC300 X-1 (Crowther et al. 2010), whose nature remains debated.

To explain why BH-HMXBs are much less numerous than Wolf-Rayet $+$ O star binaries in our Galaxy — the latter being thought to be the progenitor of the former — it has been proposed that the binaries are detectable as BH-HMXBs only when the stellar winds carry sufficient angular momentum to form an accretion disk (e.g. Vanbeveren et al. 2020; Sen et al. 2021). More recently, Hirai & Mandel (2021), through hydro-dynamical simulations, found that an accretion disk is formed when the MS donor’s Roche lobe filling factor, $f_{\rm RL}$ , defined as the ratio between the stellar radius and its Roche lobe radius, is $\gtrsim 0.8-0.9$ , providing an explanation for the high $f_{\rm RL}$ of all three known systems. Sen et al. (2024) analytically derived a disk formation criterion that does not constrain $f_{\rm RL}$ . They predicted a subgroup of X-ray bright wind-fed BH-HMXBs in relatively wide orbits, containing B-type donor stars with low $f_{\rm RL}$ . The disk can form due to the low-velocity winds of these stars. However, no such systems have been observed, and the overabundance of bright BH-HMXBs in their simulations might be related to the disk formation criterion.

Most intriguingly, based on either the disk continuum fitting or the reflection line fitting methods, the BHs in these systems are almost maximally spinning (see, e.g. McClintock et al. 2014; Reynolds 2014, 2021, and references therein). However, these methods are subject to potential uncertainties, which will be discussed in section 4.1. The binary black hole (BBH) population from gravitational-wave observations indicates low aligned spins for BHs (e.g., Abbott et al. 2019, 2023; Miller et al. 2020; Roulet et al. 2021). Arguments have been made that BH-HMXBs and BBHs originate from two distinct binary populations (Fishbach & Kalogera 2022) and due to selection effects, it is natural that these two populations are showing different properties (Liotine et al. 2023). Binary population synthesis (BPS) studies have shown that only a small fraction of BH-HMXBs result in BBH mergers (Neijssel et al. 2021; Gallegos-Garcia et al. 2022; Romero-Shaw et al. 2023).

However, the question of where the extreme spins come from has not yet been answered. The prescriptions of efficient angular momentum (AM) transport within stellar interiors, such as Tayler-Spruit (TS) dynamo (Spruit 2002), the revised prescription by Fuller et al. (2019), and the calibrated TS dynamo presented in Eggenberger et al. (2022), can explain the rotation profile of the Sun (Eggenberger et al. 2005) and the core-rotation rate of low-mass red giants and subgiants from asteroseismic measurements (Beck et al. 2012; Mosser et al. 2012; Deheuvels et al. 2014; Gehan et al. 2018). These theories predict low natal BH spins in general (Fuller & Ma 2019), but especially so for the first-born BHs in binary systems (Fragos & McClintock 2015). The second-born BHs can obtain spins through tidal interactions between their progenitors and the first-born BHs, which usually occur in tight orbits following common-envelope events. In the framework of efficient AM transport, the effective inspiral spin distribution of BBHs can be reproduced in the isolated binary evolution channel (Qin et al. 2018; Bavera et al. 2020, 2021; Belczynski et al. 2020; Olejak & Belczynski 2021).

Looking at BH-HMXBs, Qin et al. (2019) investigated two possible scenarios that can produce spinning BHs, the case A mass transfer (MT) channel and the chemical homogeneous evolution channel, and concluded that only if the internal AM transport is inefficient during the late evolution of the progenitor stars can high-spin BHs be produced, a scenario not supported yet by any other observational evidence or theoretical arguments. Another way for BHs to acquire significant spins is through super-Eddington accretion (Moreno Méndez 2011). The observed, isotropic-equivalent, X-ray luminosity of ultraluminous X-ray sources (e.g. Winter et al. 2006; Kovlakas et al. 2020) can exceed the Eddington limit by a few orders of magnitude. However, the possible geometrical beaming of the X-ray emission means that the actual accretion rate onto the BH does not have to be significantly super-Eddington (King et al. 2001; Poutanen et al. 2007; Kovlakas et al. 2022). Qin et al. (2022) explored the possibility of highly super-Eddington accretion and its impact on BH spin in the context of BH-HMXBs. This study, however, was limited in that they only considered one specific binary configuration, did not model the evolution of the binary prior to the BH formation and neglected the fact that observed BH-HMXBs, like Cygnus X-1, are currently detached.

In this work, we use the newly-developed BPS code POSYDON (Fragos et al. 2023), which incorporates extensive grids of detailed stellar structure and binary evolution models, to investigate the population of wind-fed BH-HMXBs. Among other advantages, detailed binary evolution simulations can accurately model the case A mass transfer phase. The case A mass transfer phase after the BH formation is of great importance for the formation of wind-fed BH-HMXBs, but often treated poorly in rapid BPS codes due to the lack of information about the internal structure of stars (Dorozsmai & Toonen 2024).

The paper is organized as follows. In section 2 we briefly introduce the new binary-star model grids and the details of our BPS study. In section 3 we present the results of our binary models and population properties. We discuss the implications of this work on wind-fed BH-HMXBs in section 4 and summarize the conclusions of this paper in section 5.

Table 1: Properties of Wind-fed BH-HMXBs

BH-HMXB	$M_{\rm BH}\,[M_{\odot}]$	$M_{\rm donor}\,[M_{\odot}]$	$P_{\rm{orb}}\,[\rm d]$	$\chi_{\rm BH}$	$f_{\rm RL}$	$L_{\rm x}\,[\rm erg\,s^{-1}]$	Reference
Cygnus X-1	$21.2\pm{2.2}$	$40.6^{+7.6}_{-7.1}$	$5.60$	$0.95^{+0.04}_{-0.084}$	0.997	$2.1\times 10^{37}$	1, 2, 3
LMC X-1	$10.91\pm{1.41}$	$31.79\pm{3.48}$	$3.91$	$0.897^{+0.077}_{-0.176}$	0.971	$2.3\times 10^{38}$	3, 4
M33 X-7	$11.4^{+3.3}_{-1.7}$	$38^{+22}_{-10}$	$3.45$	$0.84\pm{0.05}$ ^aa^aa $a$ Ramachandran et al. (2022) updated a BH spin of $\approx 0.6$	0.899	$(0.5-2)\times 10^{38}$	5, 6, 7

¹¹1\tablebib

(1)Orosz et al. (2011); (2)Miller-Jones et al. (2021); (3)Draghis et al. (2023); (4)Orosz et al. (2009); (5)Pietsch et al. (2006); (6)Liu et al. (2008, 2010); (7)Ramachandran et al. (2022)

2 Method

For all the computations presented in this paper, we use the publicly available BPS code POSYDON²²2We utilized the version of the POSYDON code identified by the commit hash 83a458d9d7c14cea2e371a8857bd061ac8c555c0, available at https://github.com/POSYDON-code/POSYDON/, along with the POSYDON v1.0 dataset published at https://zenodo.org/records/6655751., which incorporates single- and binary-star model grids at solar metallicity simulated with the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). The specific details of the stellar and binary physics employed in the production of precalculated grids are outlined in Fragos et al. (2023). Here we highlight the wind prescription adopted in our models. We utilized the Dutch scheme in MESA for stars with initial mass exceeding $8\,M{\odot}$ . Specifically, we applied the wind-loss prescription of de Jager et al. (1988) for cool stars with effective temperatures below $10,000\,\rm{K}$ , and Vink et al. (2001) scheme for hot stars with effective temperatures above $11,000\,\rm{K}$ . For temperatures between these two thresholds, we implemented a linear interpolation of the wind-loss rates. Additionally, for hot stars with a surface hydrogen fraction below $0.4$ , we adopted the Wolf-Rayet wind prescription proposed by Nugis & Lamers (2000).

By default, in POSYDON, accretion onto BHs is assumed to be Eddington-limited. As part of this work, beyond the publicly available grids, we compute two additional binary-star model grids of binaries consisting of a compact object and a hydrogen rich, MS star (CO-HMS grids), where we vary the accretion efficiency onto BHs. The two additional grids share the same initial parameters and input physics with the original grid. In the first new grid, we consider the limiting case of fully conservative accretion onto the BH, while in the second one, the accretion efficiency $\eta$ is established based on general relativistic radiation magnetohydrodynamic (GRRMHD) simulations (Kwan et al. in prep.). In that study, the authors carried out a series of simulations of super-Eddington disks in the magnetically arrested disk (MAD) state following the set up in Dai et al. (2018) and Thomsen et al. (2022) but around stellar-mass BHs. They varied the BH mass, spin and accretion rate $\dot{M}_{\rm acc}$ achieved in equilibrium states across different simulations to explore how these parameters affect the energy output and outflow properties. Based on the simulation results, they obtained an averaged fit for the relation between $\eta$ and the MT rate $\dot{M}_{\mathrm{tr}}$ (which is assumed to equal the sum of the mass accretion rate and the outflow rate obtained in the GRRMHD simulation):

\eta\left(\dot{m}\right)\equiv\frac{\dot{M}_{\mathrm{acc}}}{\dot{M}_{\mathrm{tr}}}=0.33\left(\dot{m}\right)^{-0.08},

(1)

for $\dot{m}\geq 10$ , where $\dot{m}$ is the MT rate in units of the Eddington accretion rate $\dot{M}_{\mathrm{Edd}}$ . In cases where $0<\dot{m}<10$ , we adopted an accretion rate of $\mathrm{min}(\dot{M}_{\mathrm{tr}},10\,\dot{M}_{\mathrm{Edd}}\,\eta(10))$ . Generally, $\eta$ ranges within $10-30\%$ and is inversely correlated with $\dot{M}_{\mathrm{tr}}$ .

For our synthetic HMXB population models, we simulate $500,000$ binaries for each accretion model, starting from two stars at zero-age MS. The primary star masses range within $10-120\,M_{\odot}$ , with its distribution determined by the initial mass function from Kroupa (2001). The secondary masses fall within the interval of $7-120\,M_{\odot}$ , following a flat distribution from the minimum to the primary star mass. The initial orbital period distribution is uniform in logarithmic space from $0.35\,\rm{d}$ to $6000\,\rm{d}$ , which is an extended prescription from Sana et al. (2013). We identify wind-fed BH-HMXBs by searching for detached binaries containing a BH and a companion star with its Roche lobe filling factor, $f_{\rm RL}$ , falling between $0.9\lesssim f_{\rm RL}<1$ , which is valid across most of the parameter variations explored by Hirai & Mandel (2021). We set a minimal mass of $10\,M_{\odot}$ for the companion star during the wind-fed phase. In all presented population models, we used the nearest-neighbor scheme for the interpolation between models in a grid (see Section 7 in Fragos et al. 2023).

To determine the explodability of the pre-supernova core and the BH mass from the core properties, we adopted the Patton & Sukhbold (2020) prescription and supernova engine from Sukhbold et al. (2016). We assume the supernova kick velocities of the BHs from core-collapse supernova follow a Maxwellian distribution with a velocity dispersion of $\sigma_{\mathrm{CCSN}}=265\,\mathrm{km\,s^{-1}}$ (Hobbs et al. 2005), rescaled by a factor of $1.4\,M_{\odot}/M_{\mathrm{BH}}$ . The natal spin of a BH is calculated based on the internal rotational profile of the BH progenitor star, right before core collapse (see Section 8.3.4 in Fragos et al. 2023). The further evolution of the BH spin parameter $\chi_{\rm BH}$ is calculated by assuming the accreted material carries the specific angular momentum at the innermost stable circular orbit (ISCO). After the accretion of an infinitesimal mass, the dimensionless BH spin parameter $\chi_{\rm BH,f}$ can be calculated as (Bardeen 1970):

\chi_{\rm BH,f}=\frac{r_{\mathrm{ISCO,i}}^{1/2}}{3}\frac{M_{\mathrm{i}}}{M_{\mathrm{f}}}\biggl{\{}4-\left[3\,r_{\mathrm{ISCO,i}}\left(\frac{M_{\mathrm{i}}}{M_{\mathrm{f}}}\right)^{2}-2\right]^{1/2}\biggl{\}},

(2)

for where $M_{\mathrm{i}}$ and $M_{\mathrm{f}}$ are the BH masses before and after accretion, respectively (with ( $M_{\mathrm{f}}-M_{\mathrm{i}})\ll\ M_{\mathrm{f}}$ ), $r_{\mathrm{ISCO,i}}$ is the initial normalized ISCO radius that can be calculated as $r_{\mathrm{ISCO,i}}=3+Z_{2}-[(3-Z_{1})(3+Z_{1}+2Z_{2})]^{1/2}$ , with $Z_{1}=1+(1-\chi_{\rm{BH,i}}^{2})^{1/3}[(1+\chi_{\rm{BH,i}})^{1/3}+(1-\chi_{\rm{BH,i}})^{1/3}]$ and $Z_{2}=(3\chi_{\rm{BH,i}}^{2}+Z_{1}^{2})^{1/2}$ . The BH spin is updated at every timestep for the binary-star models in the CO-HMS grid.

To estimate the X-ray luminosity of wind-fed BH-HMXBs, we calculate the wind mass accretion rate $\dot{M}_{\rm{w,acc}}$ onto the BHs based on the results of Hirai & Mandel (2021), who integrate the mass flux through the accretion radius in their simulations. The mass accretion fraction, defined as the ratio of the mass accretion rate to the total mass-loss rate from the donor star, exhibits a similar trend as the Bondi–Hoyle accretion (Bondi & Hoyle 1944), albeit subjects to large model uncertainties related to wind acceleration parameters and Eddington factors (Hirai & Mandel 2021). As a rough estimate, we adopt the values of the mass accretion fraction at $f_{\rm RL}\to 1$ from Hirai & Mandel (2021) based on their reference parameter configuration, and subsequently interpolate or extrapolate across various mass ratios. Then, the X-ray luminosity is given as

L_{\mathrm{x}}=\epsilon\dot{M}_{\rm{w,acc}}c^{2},

(3)

where $\epsilon=0.1$ is the radiative efficiency of accretion.

3 Results

3.1 Individual Wind-fed BH-HMXBs

We first illustrate a reference binary-star evolutionary model, which at some point in its evolution appears as a bright wind-fed BH-HMXB system. The reference binary initially consists of a primary star of $34.9\,M_{\odot}$ , a secondary star of $28.7\,M_{\odot}$ , with an orbital period of $10.4\,\rm{d}$ . The binary went through a case A MT phase and subsequently a case B MT phase. The primary star formed a BH of $12.2\,M_{\odot}$ with an initial spin of $0.047$ . The secondary lost mass due to stellar winds but gained $\approx 0.8\,M_{\odot}$ from accretion. In Figure 1, we show the evolution of the binary from $\approx 0.3\,\rm{Myr}$ before it enters the wind-fed HMXB phase until the stage when the donor star undergoes significant contraction during its helium-burning phase. We can see that as the donor star evolves, it keeps expanding until it fills its Roche lobe. It stays at the wind-fed HMXB phase for a short period of $\approx 0.1\,\rm{Myr}$ before the Roche-lobe overflow (RLO) phase. After a fast RLO phase, the donor star loses nearly half of the envelope mass and its surface helium mass fraction increases by approximately a factor of two. This is consistent with the observations of Cygnus X-1, where the donor star’s surface helium abundance is enhanced by more than a factor of two (Shimanskii et al. 2012). Its center helium mass fraction has increased to $0.9$ , indicating that the star is approaching the end of its MS. At this stage, the change in the hydrogen and helium abundance profiles of the partially stripped star leads to changes in the mass-radius relation during the MS (see, e.g., Quast et al. 2019; Farrell et al. 2022). The donor star contracts due to the mass loss and the binary becomes detached again. Then, the donor star stays very close to the Roche-lobe radius and the binary appears as a wind-fed HMXB for a longer period of $\approx 0.3-0.5\,\rm{Myr}$ . Under the assumption of Eddington-limited accretion, the BH barely accretes any mass during the fast MT phase. Alternatively, the BH accretes $\simeq 2\,M_{\odot}$ and spins up to $\chi_{\rm BH}\approx 0.45$ assuming the moderately super-Eddington accretion scenario, while it reaches nearly maximum spin with fully conservative accretion. Because the transferred material that is not accreted by the BH leaves the system carrying away the BH’s specific orbital AM, the Eddington-limited accretion leads to the largest binary orbital shrinkage due to MT. Moreover, Eddington-limited accretion results in the most significant stripping of the envelope, which subsequently leads to the quickest introduction of the Wolf-Rayet wind. Consequently, the duration of the second wind-fed phase is the shortest for the Eddington-limited accretion scenario. In the scenario of fully conservative accretion, the donor star experiences the least stripping and continues expanding after the fast RLO. As a result, following a long wind-fed phase, the binary enters a second slow RLO phase during the late MS stage of the donor star. Regardless of the differences, for all the three accretion scenarios, the binary spends more time as an observable wind-fed HMXB after the RLO phase, rather than before.

3.2 Population of Wind-fed BH-HMXBs

Refer to caption — Figure 1: Evolution of an example binary through the wind-fed HMXB phase under three BH accretion scenarios. The top panels show the evolution of the donor star radius (black solid line), Roche lobe radius (orange dashed line), 0.9 times the Roche lobe radius (orange dotted line), and orbital period (grey solid line). The green and blue areas indicate the duration of the wind-fed BH-HMXB phase before and after MT via RLO, respectively. The plus symbol denotes the end of the MS for the donor star. The bottom panels show the evolution of the BH mass (black solid line), donor star mass (orange solid line), donor star convective core mass (orange dotted line), and BH spin (grey solid line).

The individual binary models described in Section 3.1 provide us with some insights into the formation of wind-fed BH-HMXBs. However, it is a population study that can reveal whether the features identified in the evolution of a few individual binary models are representative of the whole population and allow us to evaluate their significance. Binaries that remain in a bright wind-fed BH-HMXB phase over a longer period of time are more likely to be detected. Consequently, we present distributions of the population properties for wind-fed BH-HMXBs in our simulations, weighted by the duration of the wind-fed HMXB phase for each binary.

While in Section 3.1 we described the evolutionary history of one reference binary, several variations of evolutionary pathways are present in our BH-HMXB population models. For the moderately super-Eddington accretion model, for example, we found that, prior to BH formation, $48.3\%$ of binaries went through both case A and case B MT, $21.6\%$ only experienced case A MT, $11.9\%$ avoided any MT, $11.6\%$ underwent case B MT, and $6.5\%$ entered a contact phase. The first two channels dominate the population, characterized by similar initial conditions. The binaries that avoided MT contain a relatively massive primary star above $\approx 60\,M_{\odot}$ (Bavera et al. 2023). However, due to strong stellar winds, the BHs originating from these massive stars do not exhibit a significantly different mass distribution. Among these, the case B MT channel is the only one showing a longer pre-RLO duration. These binaries tend to have wide orbits before RLO and more evolved secondary stars, resulting in a shorter post-RLO wind-fed phase.

Figure 2 shows the expected distributions of orbital periods, component masses, BH spins and X-ray luminosity for an observable wind-fed BH-HMXB population under the moderately super-Eddington accretion model. The figure distinguishes wind-fed BH-HMXBs that are observable before and after a RLO phase, shown with orange and blue color respectively. Notably, there is a higher probability of observing a wind-fed BH-HMXB after the occurrence of MT via RLO, constituting $63\%$ of the total duration. Romero-Shaw et al. (2023) also mention that, in their BPS models, some binaries would engage in a second wind-fed BH-HMXB phase after MT via RLO, when the mass ratio is close to unity. Our detailed binary models show that, when weighted by duration, the binaries that are in the post-RLO wind-fed BH-HMXB phase are predominant, characterized by a broader range of mass ratios. This can be seen in the subplot of $M_{\rm{d}}$ versus $M_{\rm{BH}}$ of Figure 2, where the mass ratio can reach $\sim 2$ . Compared to rapid BPS codes, which typically move across single stellar tracks during MT phases, POSYDON models the readjustment of mass-losing stars. Consequently, the partially stripped MS donor stars can shrink more substantially than rapid BPS codes predict in general, leading to a higher contribution from the post-RLO wind-fed BH-HMXB phase. The wind-fed BH-HMXB phase following RLO onto the BH constitutes an even higher fraction (almost 80% by duration) of bright BH-HMXBs, with X-ray luminosities above $10^{37}\,\rm erg\,s^{-1}$ (see appendix C).

As the MS donor stars should be in close orbits, nearly filling the Roche lobes, the majority of the binaries exhibits orbital periods of $\sim 1-20\,\rm{d}$ , which is anticipated. The binaries that are in wide orbits extending to thousands of days in period contain post-MS donor stars. Those supergiant donor stars evolve much faster than MS stars, resulting in small contributions to the whole population. Besides, a smaller fraction of stripped stars, fitting in close orbits with orbital periods less than $1\,\rm{d}$ , meet the requirement after experiencing common envelope evolution. Overall, MS donor stars account for $\approx 90\%$ of the whole population. However, it is important to note that Hirai & Mandel (2021) focused on line-driven winds and their conclusions about disk formation are only applicable to massive MS stars. Nevertheless, our population results indicate that post-MS donor stars are unlikely to make significant contributions due to their short evolutionary timescales, even if they have less strict requirements on the Roche-lobe filling factor for disk formation.

Regarding the donor star masses, there is a clear difference between the observations and the simulation samples. The three observed binaries contain donor stars above $\approx 30\,M_{\odot}$ while majority of donor stars are below $\approx 20\,M_{\odot}$ in our simulations. The three observed BH-HMXBs have different metallicities. As a result, direct comparison of donor star masses is not robust, as the simulations were conducted at solar metallicity. The implications of the distribution of the donor star masses are discussed in Section 4. The donor stars are mostly close to the terminal-age MS. We have included the distributions of the central and surface helium abundance of the donor stars in appendix A.

The BH masses are mostly within the range of $\sim 8-17\,M_{\odot}$ , where moderately super-Eddington accretion has increased the current BH masses by a few solar masses compared to their birth values. For the BH spins, a bimodal distribution is evident, with the two peaks corresponding to the population before and after RLO. Finally, we show a rough estimate of the X-ray luminosity, acknowledging the model uncertainties associated with the accretion fraction calculated in Hirai & Mandel (2021). We find that on average the post-RLO sub-population exhibits higher X-ray luminosities compared with the pre-RLO one. This is a result of the fact that on average the post-RLO sub-population has stronger stellar winds and higher mass ratios ( $q=M_{\rm{BH}}/M_{\rm{d}}$ ) compared to the pre-RLO one.

The HMXB population properties resulting from the other two accretion models can be found in appendix B. The most notable distinctions caused by accretion are evident in BH masses and spins. In Figure 3, the top panels show the time-weighted probability distributions of the BH mass and spin for the three accretion scenarios, distinguishing between the wind-fed phases before and after MT via RLO. In the pre-RLO HMXB sub-populations, the distributions are the same, as expected. Most of the BH masses are below $\approx 15\,M_{\odot}$ , and the spins are below $\approx 0.2$ . In the case of Eddington-limited accretion, the BH mass and spin ranges remain nearly unchanged after RLO, but the mean of the BH mass shifts towards higher values. The higher normalization of the distribution denotes that on average the duration of the post-RLO wind-fed HMXB phase is longer, and that the probability of observing a wind-fed BH-HMXB after it has experienced a RLO phase is higher. Moderately super-Eddington accretion slightly increases the BH masses and spins up the BHs to a range of $\chi_{\rm{BH}}\sim 0.2-0.6$ . For fully conservative accretion, a bimodal distribution of BH masses emerges, as the BHs could accrete a substantial amount of mass through RLO to be in a range of $\sim 15-25\,M_{\odot}$ . After accretion, BH spins become extremely high, exceeding $\approx 0.9$ for the majority of them. In the bottom panels of Figure 3, we present the probability distributions of the donor star mass and the orbital period, in which the three scenarios exhibit no significant differences. The fully conservative scenario leads to more stable MT, resulting in the disappearance of stripped stars coming from common envelope evolution. Regardless of the accretion scenarios, the total duration of the wind-fed BH-HMXB phase following MT via RLO is always longer than that before MT. The fraction of post-RLO duration relative to the total duration is greater than $60\%$ for all the three accretion scenarios.

4 Discussion

4.1 Black Hole Spins

Our study has shown that, assuming efficient AM transport inside stars, the BHs in wind-fed BH-HMXBs are born with spins $\lesssim 0.2$ . Eddington-limited accretion does not lead to significant changes in the spins. The moderately super-Eddington accretion model can produce BHs with spins of $\chi_{\rm{BH}}\sim 0.3-0.6$ whereas the fully conservative MT can spin up most of the BHs to extreme spins if the accreted material carries the specific AM at ISCO. It is evident that explaining the observed BH spins in BH-HMXBs poses a challenge in the absence of substantial accretion onto the BHs. Even when considering the impact of supernova processes on the natal BH spins, such as those arising from slow ejecta (Schrøder et al. 2018) or shock instability (Moreno Méndez & Cantiello 2016), the attainment of extreme spins remains elusive. Batta et al. (2017) found that BHs with high spins of $\sim 0.8$ can be produced, but only if the fallback material reaches the position of the companion star and extracts AM from the orbit.

The moderately super-Eddington accretion model is also insufficient to reproduce the extreme BH spins in BH-HMXBs. In the GRRMHD simulations, which are used to estimate the BH accretion efficiency in our moderately super-Eddington model, ordered magnetic fields thread the super-Eddington disks. These disks eventually enter MAD state, where the AM carried by the transferred material near the BH is less than the AM at the ISCO (e.g. Lowell et al. 2024). Therefore, the increase in BH spins in our simulations should be considered as an upper limit. Additionally, BHs fed by MADs launch powerful jets that can effectively spin down the BHs when the spin is very large (Tchekhovskoy et al. 2012; Lowell et al. 2024). If the disks around the BHs in HMXBs easily reach the MAD state, it naturally excludes the possibility that the AM of the BH is obtained through gas accretion. However, if the disk does not contain enough magnetic flux to enter a MAD state, the accretion process can spin up the BH to a very high equilibrium spin value (Gammie et al. 2004). Such a disk does not have significant magnetic pressure to power magnetic outflow, and it is expected that BHs accrete a higher fraction of the disk gas from non-MAD state super-Eddington disks, which lead to higher accretion efficiency and, thus, higher BH spins. In the fully-conservative accretion model, we find that a large fraction of BHs accreted more material than their initial masses, indicating that the accretion is more than adequate to produce extreme spins.

It is essential to note that the two primary methods for measuring BH spins in HMXBs may be susceptible to systematic uncertainties (Taylor & Reynolds 2018; Salvesen & Miller 2021; Falanga et al. 2021) and model dependency (Zdziarski et al. 2024a). For example, Zdziarski et al. (2024b) argue that when considering a warm, optically-thick, thermal Comptonizing layer on top of the disk, the spin of Cygnus X-1 is measured to be low, $\chi_{\rm{BH}}\lesssim 0.3$ . Yet, if the extreme measured spins are real, to explain them in the context of super-Eddington accretion, mass and AM accretion efficiencies higher than those predicted by GRRMHD simulations of MAD-state super-Eddington disks should be considered.

If both the mass accretion and AM transport to BHs are efficient enough to produce highly spinning BHs, it would likely also influence the properties of coalescing BBH populations. During the evolution of merging BBHs, RLO accretion onto the first-born BH is anticipated. The BHs can undergo spin-up due to efficient accretion. In this case, we would not expect the first-born BHs to exhibit nearly-zero spins in BBH mergers. Shao & Li (2022) employed mildly super-Eddington accretion of case A MT onto the BHs, highlighting its capability to account for the observed highly spinning BHs in BBH mergers. Bavera et al. (2021) found that highly conservative accretion onto BHs effectively reduces the contribution from the stable MT channel for BBH mergers, as the binaries do not shrink sufficiently during MT to coalesce due to gravitational waves within the Hubble time. However, most rapid BPS codes, such as COSMIC, which was utilized by Bavera et al. (2021) to model the MT phase, are recognized for their inadequate treatment of case A mass transfer (Dorozsmai & Toonen 2024). In contrast, BPS codes that utilize detailed binary models and assume fully conservative accretion onto black holes (e.g. Briel et al. 2023) indicate that the stable MT channel remains viable. However, a direct comparison between the models presented in these studies should be done with caution, as they also differ in several other key physics assumptions. The specific properties of BBHs, distribution of effective inspiral spins, and their formation channels under efficient BH accretion scenarios should be investigated in future studies.

4.2 Donor Star Masses

The three observed wind-fed BH-HMXBs exhibit donor stars with masses $\gtrsim 30\,M_{\odot}$ . In our population studies, we found that donor stars within $\sim 10-20\,M_{\odot}$ are expected to be prevalent in wind-fed BH-HMXBs at solar metallicity. One caveat in comparing our population models to the observed properties of the three known wind-fed BH-HMXBs is that our population models assume solar metallicity, while the observed systems have varied metallicities. Specifically, M33 X-7 has been suggested to have a metallicity between $0.1\,Z_{\odot}$ (Orosz et al. 2007) and $0.5\,Z_{\odot}$ (Ramachandran et al. 2022). Similarly, LMC X-1 is expected to have sub-solar metallicity, consistent with the LMC metallicity. For these two systems, comparisons with our simulations are not straightforward, as lower metallicities can lead to the formation and survival of more massive stars. For Cygnus X-1, however, a solar or even a supersolar metallicity (Shimanskii et al. 2012) is suggested, making its $40\,M{\odot}$ donor star hard to reconcile with our simulations (but see discussion in Neijssel et al. (2021)).

The prevalence of $\sim 10-20\,M_{\odot}$ stars in our simulations can be explained as follows. Firstly, the accretion onto non-degenerate stars in POSYDON is regulated by stellar rotation, typically leading to low accretion efficiency (Fragos et al. 2023). Consequently, the secondary stars generally do not gain substantial mass before the BH formation. Such inefficient accretion is suggested to prevent the formation of merging NSBH with first-born NSs (Xing et al. 2024) and potentially contributes to discrepancies between the predicted and observed Be X-ray binary donor mass distributions (Vinciguerra et al. 2020; Akira Rocha et al. 2024). In the case of wind-fed BH-HMXBs, if the mass accretion efficiency is higher during earlier evolutionary stages, the donor star could have accreted more mass before BH formation, potentially alleviating the tension between our simulations and the observations.

Secondly, at solar metallicity, when the primary star becomes a BH, the secondary star should have lost a significant fraction of its mass through stellar winds if it is initially a massive star $\gtrsim 40\,M_{\odot}$ . Moreover, many BH-HMXBs should be in the second wind-fed phase after the donor stars have transferred mass to the BHs through RLO. Romero-Shaw et al. (2023) produced wind-fed BH-HMXBs with donor stars that can reach $\sim 70\,M_{\odot}$ at solar metallicity. In our simulations, when the BH forms, the massive donor star has lost a substantial amount of mass due to stellar winds and mass transfer, reducing its mass to below $\approx 50\,M_{\odot}$ . The major difference arises from the modelling of stellar winds of massive stars at the MS stage. We found that massive stars $\gtrsim 80\,M_{\odot}$ at solar metallicity have strong winds that can deplete the hydrogen envelope and induce Wolf–Rayet winds during their MS stage (Bavera et al. 2023). As a result, the maximum donor star mass in our simulations is below $50\,M_{\odot}$ .

Moreover, according to the IMF, the number of low-mass stars surpasses the number of high-mass stars. With a flat distribution of the initial mass ratio and without any selection effect on the donor stars, it is natural to expect a prevalence of $\sim 10-20\,M_{\odot}$ donor stars as they are more common and have longer lifetimes compared with a $\sim 40\,M_{\odot}$ star. The absence of massive donor stars was also indicated by Sen et al. (2024), who found X-ray bright binaries do not include contributions from initially $\gtrsim 40\,M_{\odot}$ donor stars.

The differences between modelled and observed donor star masses might indicate that certain physical assumptions or treatments in our stellar and binary evolution are not robust. Wind prescriptions, mass accretion efficiency for non-degenerate stars, and AM loss from non-conservative MT, as suggested in Romero-Shaw et al. (2023), could all impact the population properties. Therefore, it is imperative to investigate them further, seeking a consistent model that can be applied to diverse binary systems.

Another avenue worth exploring involves investigating the selection effects that favor massive donor stars. Our estimation of X-ray luminosity is too rough to be quantitatively analyzed as we only adopt one set of parameters for estimating the accretion efficiency for all binaries from Hirai & Mandel (2021) and interpolate or extrapolate among different mass ratios. Moreover, the mass accretion rate obtained in Hirai & Mandel (2021) is an upper limit because not all mass passing through the accretion radius is guaranteed to be accreted. In our simulations, most donor stars with masses within $\sim 10-20\,M_{\odot}$ can produce X-ray luminosity over $10^{37}\,\rm{erg\,s^{-1}}$ . However, if we apply the low-end efficiency in Hirai & Mandel (2021) to all binaries, those with X-ray luminosity above $10^{37}\,\rm{erg\,s^{-1}}$ account for $\approx 21\%$ of the total population and $\approx 56\%$ of them contain a donor star greater than $20\,M_{\odot}$ . A corner plot showing the property distributions for these X-ray bright binaries can be found in appendix C.

5 Conclusions

The three observed wind-fed BH-HMXBs demonstrate high BH spins, the origin of which remains an open question. We used the new-generation BPS code POSYDON to investigate the population properties of wind-fed BH-HMXBs, adopting three CO-HMS grids with different assumed accretion efficiencies onto BHs. We summarize the key conclusions here:

1.

Wind-fed BH-HMXBs at solar metallicity are more likely to have already been though a fast RLO phase, after the BH formation. Regardless of the mass-accretion models, in our simulations, the total duration of wind-fed BH-HMXBs in a post-RLO phase is higher than that in a pre-RLO phase. After mass transfer via RLO, the partially stripped donor star contracts due to the changes in its hydrogen and helium abundance profile, remaining close to its Roche lobe. At this stage, the binary stays detached for a long period, appearing as a wind-fed BH-HMXB. We highlight the importance of employing detailed modelling in BPS studies, which includes a self-consistent treatment of MT and subsequent readjustment of stars experiencing mass loss or gain.
2.

An accretion efficiency prescription inspired by GRRMHD simulations can produce BHs with moderate spins in the range of $\sim 0.3-0.6$ , while fully conservative MT can spin up BHs to near-maximal spins, when assuming that the accreted material carries the specific AM at the ISCO. To explain the high spins observed in the three wind-fed BH-HMXBs by super-Eddington accretion, a mass and AM accretion efficiency onto BHs higher that that predicted by GRRMHD simulation of super-Eddington disks in the MAD state is required.
3.

In our simulations, MS donor stars with masses in the range of $\sim 10-20\,M_{\odot}$ dominate the wind-fed BH-HMXB population at solar metallicity. This is in tension with the fact that the only confirmed wind-fed BH-HMXB in our Galaxy, Cygnus X-1, contains a donor star of $\approx 40\,M_{\odot}$ . Assuming the mass measurements are reliable, we need to explore further observational selection effects that favor massive donor stars, as well as consider an accretion efficiency model onto non-degenerate stars that leads to more conservative MT.

Acknowledgements.

The POSYDON project is supported primarily by two sources: the Swiss National Science Foundation (PI Fragos, project numbers PP00P2_211006 and CRSII5_213497) and the Gordon and Betty Moore Foundation (PI Kalogera, grant award GBMF8477). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie RISE action, grant agreement No 873089 (ASTROSTAT- II). Z.X. acknowledges support from the China Scholarship Council (CSC). E.Z. acknowledges funding support from the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the ”3rd Call for H.F.R.I. Research Projects to support Post-Doctoral Researchers” (Project No: 7933). S.S.B., T.F., and Z.X. were supported by the project number PP00P2_211006. T.K. and L.D. acknowledge the support from the National Natural Science Foundation of China and the Hong Kong Research Grants Council (12122309, N_HKU782/23, 17305523, 17314822). S.S.B. was also supported by the project number CRSII5_213497. M.B. ackowledges support from the Boninchi Foundation. K.A.R. and M.S. are supported by the Gordon and Betty Moore Foundation (PI Kalogera, grant award GBMF8477) K.A.R. is also supported by the Riedel Family Fellowship. K.A.R. also thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant No. 1829740, the Brinson Foundation, and the Moore Foundation; their participation in the program has benefited this work. K.K. acknowledges support from the Spanish State Research Agency, through the María de Maeztu Program for Centers and Units of Excellence in R&D, No. CEX2020-001058-M. J.J.A. acknowledges support for Program number (JWST-AR-04369.001-A) provided through a grant from the STScI under NASA contract NAS5-03127. I.M. acknowledges support from the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE230100016.

References

Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, ApJ, 882, L24
Abbott et al. (2023) Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Physical Review X, 13, 011048
Akira Rocha et al. (2024) Akira Rocha, K., Kalogera, V., Doctor, Z., et al. 2024, arXiv e-prints, arXiv:2403.07172
Bardeen (1970) Bardeen, J. M. 1970, Nature, 226, 64
Batta et al. (2017) Batta, A., Ramirez-Ruiz, E., & Fryer, C. 2017, ApJ, 846, L15
Bavera et al. (2020) Bavera, S. S., Fragos, T., Qin, Y., et al. 2020, A&A, 635, A97
Bavera et al. (2023) Bavera, S. S., Fragos, T., Zapartas, E., et al. 2023, Nature Astronomy, 7, 1090
Bavera et al. (2021) Bavera, S. S., Fragos, T., Zevin, M., et al. 2021, A&A, 647, A153
Beck et al. (2012) Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55
Belczynski et al. (2020) Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, A&A, 636, A104
Bondi & Hoyle (1944) Bondi, H. & Hoyle, F. 1944, MNRAS, 104, 273
Briel et al. (2023) Briel, M. M., Stevance, H. F., & Eldridge, J. J. 2023, MNRAS, 520, 5724
Crowther et al. (2010) Crowther, P. A., Barnard, R., Carpano, S., et al. 2010, MNRAS, 403, L41
Dai et al. (2018) Dai, L., McKinney, J. C., Roth, N., Ramirez-Ruiz, E., & Miller, M. C. 2018, ApJ, 859, L20
de Jager et al. (1988) de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259
Deheuvels et al. (2014) Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27
Dorozsmai & Toonen (2024) Dorozsmai, A. & Toonen, S. 2024, MNRAS[arXiv:2207.08837]
Draghis et al. (2023) Draghis, P. A., Miller, J. M., Costantini, E., et al. 2023, arXiv e-prints, arXiv:2311.16225
Eggenberger et al. (2005) Eggenberger, P., Maeder, A., & Meynet, G. 2005, A&A, 440, L9
Eggenberger et al. (2022) Eggenberger, P., Moyano, F. D., & den Hartogh, J. W. 2022, A&A, 664, L16
Falanga et al. (2021) Falanga, M., Bakala, P., La Placa, R., et al. 2021, MNRAS, 504, 3424
Farrell et al. (2022) Farrell, E., Groh, J. H., Meynet, G., & Eldridge, J. J. 2022, MNRAS, 512, 4116
Fishbach & Kalogera (2022) Fishbach, M. & Kalogera, V. 2022, ApJ, 929, L26
Fragos et al. (2023) Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2023, ApJS, 264, 45
Fragos & McClintock (2015) Fragos, T. & McClintock, J. E. 2015, ApJ, 800, 17
Fuller & Ma (2019) Fuller, J. & Ma, L. 2019, ApJ, 881, L1
Fuller et al. (2019) Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661
Gallegos-Garcia et al. (2022) Gallegos-Garcia, M., Fishbach, M., Kalogera, V., L Berry, C. P., & Doctor, Z. 2022, ApJ, 938, L19
Gammie et al. (2004) Gammie, C. F., Shapiro, S. L., & McKinney, J. C. 2004, ApJ, 602, 312
Gehan et al. (2018) Gehan, C., Mosser, B., Michel, E., Samadi, R., & Kallinger, T. 2018, A&A, 616, A24
Gomez & Grindlay (2021) Gomez, S. & Grindlay, J. E. 2021, ApJ, 913, 48
Hirai & Mandel (2021) Hirai, R. & Mandel, I. 2021, Publications of the Astronomical Society of Australia, 38, e056
Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., & Kramer, M. 2005, MNRAS, 360, 974
Jermyn et al. (2023) Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15
King et al. (2001) King, A. R., Davies, M. B., Ward, M. J., Fabbiano, G., & Elvis, M. 2001, ApJ, 552, L109
Kovlakas et al. (2022) Kovlakas, K., Fragos, T., Schaerer, D., & Mesinger, A. 2022, A&A, 665, A28
Kovlakas et al. (2020) Kovlakas, K., Zezas, A., Andrews, J. J., et al. 2020, MNRAS, 498, 4790
Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
Kwan et al. (in prep.) Kwan, T. M., Dai, L., Kan, Z. C. K., et al. in prep.
Laycock et al. (2015) Laycock, S. G. T., Maccarone, T. J., & Christodoulou, D. M. 2015, MNRAS, 452, L31
Liotine et al. (2023) Liotine, C., Zevin, M., Berry, C. P. L., Doctor, Z., & Kalogera, V. 2023, ApJ, 946, 4
Liu et al. (2008) Liu, J., McClintock, J. E., Narayan, R., Davis, S. W., & Orosz, J. A. 2008, ApJ, 679, L37
Liu et al. (2010) Liu, J., McClintock, J. E., Narayan, R., Davis, S. W., & Orosz, J. A. 2010, ApJ, 719, L109
Lowell et al. (2024) Lowell, B., Jacquemin-Ide, J., Tchekhovskoy, A., & Duncan, A. 2024, ApJ, 960, 82
McClintock et al. (2014) McClintock, J. E., Narayan, R., & Steiner, J. F. 2014, Space Sci. Rev., 183, 295
Miller et al. (2020) Miller, S., Callister, T. A., & Farr, W. M. 2020, ApJ, 895, 128
Miller-Jones et al. (2021) Miller-Jones, J. C. A., Bahramian, A., Orosz, J. A., et al. 2021, Science, 371, 1046
Moreno Méndez (2011) Moreno Méndez, E. 2011, MNRAS, 413, 183
Moreno Méndez & Cantiello (2016) Moreno Méndez, E. & Cantiello, M. 2016, New A, 44, 58
Mosser et al. (2012) Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10
Neijssel et al. (2021) Neijssel, C. J., Vinciguerra, S., Vigna-Gómez, A., et al. 2021, ApJ, 908, 118
Nugis & Lamers (2000) Nugis, T. & Lamers, H. J. G. L. M. 2000, A&A, 360, 227
Olejak & Belczynski (2021) Olejak, A. & Belczynski, K. 2021, ApJ, 921, L2
Orosz et al. (2011) Orosz, J. A., McClintock, J. E., Aufdenberg, J. P., et al. 2011, ApJ, 742, 84
Orosz et al. (2007) Orosz, J. A., McClintock, J. E., Narayan, R., et al. 2007, Nature, 449, 872
Orosz et al. (2009) Orosz, J. A., Steeghs, D., McClintock, J. E., et al. 2009, ApJ, 697, 573
Patton & Sukhbold (2020) Patton, R. A. & Sukhbold, T. 2020, MNRAS, 499, 2803
Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
Pietsch et al. (2006) Pietsch, W., Haberl, F., Sasaki, M., et al. 2006, ApJ, 646, 420
Poutanen et al. (2007) Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G., & Abolmasov, P. 2007, MNRAS, 377, 1187
Prestwich et al. (2007) Prestwich, A. H., Kilgard, R., Crowther, P. A., et al. 2007, ApJ, 669, L21
Qin et al. (2018) Qin, Y., Fragos, T., Meynet, G., et al. 2018, A&A, 616, A28
Qin et al. (2019) Qin, Y., Marchant, P., Fragos, T., Meynet, G., & Kalogera, V. 2019, ApJ, 870, L18
Qin et al. (2022) Qin, Y., Shu, X., Yi, S., & Wang, Y.-Z. 2022, Research in Astronomy and Astrophysics, 22, 035023
Quast et al. (2019) Quast, M., Langer, N., & Tauris, T. M. 2019, A&A, 628, A19
Ramachandran et al. (2022) Ramachandran, V., Oskinova, L. M., Hamann, W. R., et al. 2022, A&A, 667, A77
Reynolds (2014) Reynolds, C. S. 2014, Space Sci. Rev., 183, 277
Reynolds (2021) Reynolds, C. S. 2021, ARA&A, 59, 117
Romero-Shaw et al. (2023) Romero-Shaw, I., Hirai, R., Bahramian, A., Willcox, R., & Mandel, I. 2023, MNRAS, 524, 245
Roulet et al. (2021) Roulet, J., Chia, H. S., Olsen, S., et al. 2021, Phys. Rev. D, 104, 083010
Salvesen & Miller (2021) Salvesen, G. & Miller, J. M. 2021, MNRAS, 500, 3640
Sana et al. (2013) Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107
Schrøder et al. (2018) Schrøder, S. L., Batta, A., & Ramirez-Ruiz, E. 2018, ApJ, 862, L3
Seifina & Titarchuk (2010) Seifina, E. & Titarchuk, L. 2010, ApJ, 722, 586
Sen et al. (2024) Sen, K., El Mellah, I., Langer, N., et al. 2024, arXiv e-prints, arXiv:2406.08596
Sen et al. (2021) Sen, K., Xu, X. T., Langer, N., et al. 2021, A&A, 652, A138
Shao & Li (2022) Shao, Y. & Li, X.-D. 2022, ApJ, 930, 26
Shimanskii et al. (2012) Shimanskii, V. V., Karitskaya, E. A., Bochkarev, N. G., et al. 2012, Astronomy Reports, 56, 741
Spruit (2002) Spruit, H. C. 2002, A&A, 381, 923
Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38
Taylor & Reynolds (2018) Taylor, C. & Reynolds, C. S. 2018, ApJ, 855, 120
Tchekhovskoy et al. (2012) Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2012, in Journal of Physics Conference Series, Vol. 372, Journal of Physics Conference Series, 012040
Thomsen et al. (2022) Thomsen, L. L., Kwan, T. M., Dai, L., et al. 2022, ApJ, 937, L28
Vanbeveren et al. (2020) Vanbeveren, D., Mennekens, N., van den Heuvel, E. P. J., & Van Bever, J. 2020, A&A, 636, A99
Vinciguerra et al. (2020) Vinciguerra, S., Neijssel, C. J., Vigna-Gómez, A., et al. 2020, MNRAS, 498, 4705
Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
Winter et al. (2006) Winter, L. M., Mushotzky, R. F., & Reynolds, C. S. 2006, ApJ, 649, 730
Xing et al. (2024) Xing, Z., Bavera, S. S., Fragos, T., et al. 2024, A&A, 683, A144
Zdziarski et al. (2024a) Zdziarski, A. A., Banerjee, S., Chand, S., et al. 2024a, ApJ, 962, 101
Zdziarski et al. (2024b) Zdziarski, A. A., Chand, S., Banerjee, S., et al. 2024b, ApJ, 967, L9
Zdziarski et al. (2013) Zdziarski, A. A., Mikolajewska, J., & Belczynski, K. 2013, MNRAS, 429, L104

Appendix A Surface and Center Helium Abundance

In this section, we present the population property distributions, in the form of a corner plot, that includes center helium abundance $Y_{\rm{center}}$ and surface helium abundance $Y_{\rm{surf}}$ of the donor star under the moderately super-Eddington accretion scenario. In Figure 4, we show the time-weighted distribution of orbital period, donor star mass, black hole spin, donor star’s center helium abundance, and surface helium abundance. For both pre-RLO and post-RLO systems, the majority of the donor stars have $Y_{\rm{center}}$ above $\approx 0.8$ , which means that they are approaching the terminal-age MS. As a consequence of RLO, the outer layer of the envelope of the donor stars is stripped, leading to an evident increase in $Y_{\rm{surf}}$ .

Appendix B Population Property Distributions for Eddington-limited and Fully Conservative BH Accretion

Here we display population property distributions, in the form of corner plots, for the population properties of wind-fed BH-HMXBs under the model assumptions of Eddington-limited and fully conservative BH accretion. Figure 5 shows the same distributions as Figure 2 but for the Eddington-limited accretion model. The BH spins are $\lesssim 0.3$ , with no significant BH spin-up observed to distinguish the populations before and after MT via RLO. Figure 6 shows the case of fully conservative BH accretion. Bimodal distributions can be seen for BH masses and spins due to accretion. Most of the BHs can be spun up to have spins above $0.8$ after MT via RLO and the binaries remain at this stage longer. The distributions of orbital periods and donor star masses appear similar across different BH accretion models.

Appendix C Population Property Distributions for X-ray bright sources

When adopting the low-end wind accretion efficiency from Hirai & Mandel (2021), a large fraction of the binaries with donor stars below $20\,M_{\odot}$ exhibit X-ray luminosities below $10^{37}\,\rm{erg\,s^{-1}}$ . To illustrate the characteristics of X-ray bright systems, we display a corner plot of the binaries with X-ray luminosity exceeding $10^{37}\,\rm{erg\,s^{-1}}$ in Figure 7. About $79\%$ of them are post-RLO systems. The donor stars exhibit a relatively flat distribution from $10$ to $30\,M_{\odot}$ , extending up to $\approx 40\,M_{\odot}$ . The X-ray bright systems are characterized by a higher proportion of massive donor stars because of their strong winds.

Formation of Wind-Fed Black Hole High-mass X-ray Binaries: The Role of Roche-lobe-Overflow Post Black-Hole Formation