VFTS 243 as predicted by the BPASS fiducial models

ABSTRACT

The recent discovery of an unambiguous quiescent black hole (BH) and main-sequence O star companion in VFTS 243 opens the door to new constraints on theoretical stellar evolution and population models looking to reproduce the progenitors of BH–BH binaries. Here, we show that the binary population and spectral synthesis fiducial models natively predict VFTS 243-like systems: We find that VFTS 243 likely originated from a binary system in an ∼15 d orbit with primary mass ranging from 40 to 50 |$\mathrm{\, M}_\odot$| and secondary star with initial mass of 24–25 |$\mathrm{\, M}_\odot$|⁠. Additionally, we find that the death of the primary star must have resulted in a low-energy explosion E < 10⁵⁰ erg. With a uniform prior, we find that the kick velocity of the newborn BH was ≤10 |$\, \text{km s}^{-1}$|⁠. The very low eccentricity reported for VFTS 243 and the subsequent conclusion by the authors that the supernova kick must have been very small are in line with the peak in the posterior distribution between 0 and 5 |$\, \text{km s}^{-1}$| found from our numerical simulations performed with a uniform prior. Finally, the reduced Hobbs kick distribution commonly used in BH population synthesis is strongly disfavoured.

stars: evolution, stars: binaries: general, stars: black holes, stars: individual: VFTS 243

1 INTRODUCTION

The majority of massive stars (born with a mass greater than 8 |$\mathrm{\, M}_\odot$|⁠) are present in binary systems (Sana et al. 2012, 2014; Moe & Di Stefano 2017), and on exhausting their nuclear supply, they collapse under their gravitational pressure to form a compact remnant. If the collapsing star is sufficiently massive, this can create a binary system composed of a bright massive star in orbit around a black hole (BH). Such binary systems could be strong emitters of X-rays fuelled by the accretion of matter from the star on to the BH companion (Shakura & Sunyaev 1973). Such a distinctive observational feature makes the identification of the binary system hosting the BH possible. Additionally, attempts have been made to identify star + BH binaries in systems where the BH is quiescent (not X-ray loud) through spectroscopic studies. In recent years, a number of such systems have been reported (LB-1, Liu et al. 2019; HR 6819, Rivinius et al. 2020; NGC 1850-BH1, Saracino et al. 2021), only to be later suspected to harbour a stripped star as a BH impostor (Bodensteiner et al. 2020; El-Badry & Quataert 2020; Shenar et al. 2020; El-Badry & Burdge 2021; Stevance, Parsons & Eldridge 2022; see Bodensteiner et al. 2022, for a review).

The VFTS 243 system, however, seems to be the first undisputed discovery of a quiescent BH in a stellar binary. Reported by Shenar et al. (2022a), it is found to be an O7 + BH in the 30 Doradus region. The very precise measurements of the orbital components (P = 10.4031 ± 0.0004 d, e = 0.017 ± 0.012) and extensive observational data surrounding this object have the potential to provide improved observational constraints or clues on a number of bottlenecks in stellar evolution theory such as natal supernova (SN) kicks and the genealogies of the progenitor systems of BH + BH binaries. Detectable gravitational wave sources are, for the most part, BH–BH merger systems (The LVK Collaboration 2021), but our collective understanding of their progenitor routes remains hazy (Mandel & Broekgaarden 2021). This is in part due to the wide array of possible channels and the degeneracy between these routes when the only observations available are at the beginning (binary star systems) and the end (merger of compact objects) of the life of the systems. Being able to study binaries comprised of a star and a compact remnant is therefore essential to anchor our theoretical models to observations at a crucial intermediate stage in the life of gravitational wave progenitors.

In the supplementary material of the VFTS 243 BH discovery paper, Shenar et al. (2022a) present an exemplar system calculated with mesa (Paxton et al. 2011, 2013, 2015, 2018, 2019) to show a possible progenitor route to the formation of VFTS 243-like objects. The flexibility of mesa allows observational parameters to be reproduced rather nicely (generally speaking and in the case of VFTS 243), but in order to truly understand the most likely genealogies of such systems population synthesis is required. The binary population and spectral synthesis (BPASS) project combines detailed binary stellar evolution (using a custom version of the Cambridge stars code; Eggleton 1971; Eldridge, Izzard & Tout 2008) with population synthesis based on observationally inferred binary distributions. The fiducial models (BPASS v.2.2.1) were released in 2018 (Stanway & Eldridge 2018) and have since been shown to reproduce a large number of observable phenomena, particularly in the realm of massive stars (e.g. Massey et al. 2021). The key feature is that BPASS results are not tailored to fit specific observations – they are pre-calculated models from a grid of initial parameters. This allows us to take a different approach when comparing observational systems and theory: instead of creating a model to match the observations as best as possible, we can search through our existing predictions and find systems that are similar to those observed (e.g. Xiao et al. 2019; Stanway, Eldridge & Chrimes 2020; Stevance, Eldridge & Stanway 2020; Tang et al. 2020; Byrne, Stanway & Eldridge 2021; Stevance & Eldridge 2021; Briel, Stevance & Eldridge 2022; Chrimes et al. 2022; Ghodla et al. 2022; Stevance et al. 2022).

In this work, we search the fiducial BPASS models for predictions of VFTS 243-like systems and compare our genealogies to those of the exemplar mesa model presented in Shenar et al. (2022a). In addition, we investigate the natal kick velocities that best reproduce the orbital configuration of VFTS 243 and compare these to commonly used natal kick prescriptions for BH remnants. In Section 2, we describe our models and numerical methods; in Section 3, we show our posterior kick velocities and best genealogies to VFTS 243 as predicted by BPASS. In Section 4, we discuss the differences between our channel and the mesa model in Shenar et al. (2022a), as well as the implications of our kick posteriors. Finally, we summarize and conclude in Section 5.

2 MODELS AND NUMERICAL METHODS

We use the fiducial BPASS v2.2.1 models (described in Eldridge et al. 2017; Stanway & Eldridge 2018): each population is born from an instantaneous starburst of 10⁶|$\mathrm{\, M}_\odot$| at a single metallicity; in the present case, we use Z = 0.006 as it best corresponds to the average metallicity of the Large Magellanic Cloud (LMC; Massey et al. 2021). Although a wide scatter of metallicities is present in the LMC (Narloch et al. 2022), searches performed with Z = 0.010 BPASS models revealed analogous pathways and so the specific choice of metallicity does not impact our results. The initial mass function is a Kroupa prescription (Kroupa 2001) with a maximum mass of 300 |$\mathrm{\, M}_\odot$|⁠. The binary fraction and initial parameters of the binary models (mass ratio and period) are taken from the observational survey of Moe & Di Stefano (2017). SN explosions cannot be modelled in detail, but BPASS data products include ejecta masses and remnant masses for three likely scenarios: a weak SN explosion (10⁵⁰ erg), a typical SN (10⁵¹ erg), and a very energetic explosion (10⁵² erg). For the purposes of this investigation, we performed numerical simulations and model searches for all three scenarios as BHs can originate from very powerful SNe (e.g. in the collapsar scenario; Woosley 1993; Woosley & Heger 2012) or result in failed SNe (Smartt 2009; Reynolds, Fraser & Gilmore 2015). These alternative scenarios will change the remnant and ejecta masses of the post-SN model that will only indirectly affect the final orbit when using reduced Hobbs or ultra-stripped SN kick prescriptions (see below), but directly affect the kick value if using a Bray & Eldridge (2016) prescription where the ejecta mass and remnant mass are parameters included in the definition.

In order to take into account the effect of SN kicks, we use an auxiliary code called tui (previously used and described, although not mentioned by name at the time, in Ghodla et al. 2022). In this work, we implement new kick prescriptions that focus on the specific case of BH remnants. A common kick prescription is a reduced kick (e.g. Team COMPAS et al. 2021) where the BH remnants are assumed to receive the same momentum as the neutron star remnants. In this scenario, the kick value is scaled by 1.4/M_BH, where 1.4 is taken as the typical mass of a neutron star. One of the most commonly used neutron star kick distribution is a Maxwellian with σ = 265|$\, \text{km s}^{-1}$|⁠, also called a Hobbs distribution (Hobbs et al. 2005); additionally, for low ejecta masses some compact object population synthesis models include the effects of ultra-stripped SNe (Tauris et al. 2013), often with a default value of 30 |$\, \text{km s}^{-1}$|⁠. We combine these considerations into a split kick prescription, where for an SN with an ejecta mass >0.35 |$\mathrm{\, M}_\odot$| (Yao et al. 2020), we use the BH mass-scaled Hobbs distribution, while below 0.35 |$\mathrm{\, M}_\odot$| we impart a kick of 30 |$\, \text{km s}^{-1}$| that we also scale by the mass of the remnant. In the final matches we find to VFTS 243, no progenitor matched the ultra-stripped SN criteria. To further investigate the kick distribution of VFTS 243-like systems without any prior assumptions, we also perform our analysis with kicks drawn from a uniform distribution between 0 and 100 |$\, \text{km s}^{-1}$| – note that no mass scaling is applied in this case.

For each model undergoing an SN (i.e. obeys the condition that its CO core mass >1.38 |$\mathrm{\, M}_\odot$|⁠), we draw 2000 kicks and record the models that remain bound. In addition, we also consider a rotational velocity criterion as the O star in VFTS 243 was inferred to have a spin velocity of at least 180 |$\, \text{km s}^{-1}$| (Shenar et al. 2022a). The original stars models in BPASS do not take into account rotation in the evolution of the stellar interior (although it is considered in the orbital evolution). Using a disc-to-star angular momentum accretion efficiency (ν) of 0.5 (i.e. 50 per cent), any star accreting mass |$\Delta M \gt M r_{\mathrm{ g}}^{2} \frac{180 / \it {v}_{ \rm cr}}{\nu }$| is considered to be spun up to velocity >180 |$\, \text{km s}^{-1}$|⁠. Here, |$M, r_{\mathrm{ g}}, \text{ and }\rm {\mathit{ v}}_{\rm cr}$| are the pre-accretion mass, radius of gyration, and the critical velocity of the star, respectively, and it is assumed that the star is rotating as a solid body. This expression has been adopted from Ghodla et al. (2023), who studied the spin-up efficiency in the secondary component of a binary system due to accretion of mass from a Keplerian disc around the accretor. Since the projection of the VFTS 243 orbit is unknown, we allow all velocities larger than 180 |$\, \text{km s}^{-1}$| (up to critical) in our working.

3 RESULTS

3.1 Posterior kick velocities

In the next section, we will search for the systems that are most similar to VFTS 243, but first, to explore the impact of kicks, we only consider the period and eccentricity criteria. Shenar et al. (2022a) find a period of 10.4 d and an eccentricity of 0.017 |$\pm \, 0.012$|⁠; additionally, the rotation period they infer for the O star is not synchronized with the orbit and therefore tidal interactions (which would circularize the orbit) can be neglected. Consequently, we need to look for models that have a similar orbital configuration after the first SN, without undergoing further binary interactions to reprocess the orbit. Here, we consider all models that have a period between 9 and 11 d and an eccentricity lower than 0.05, after SN kicks have been taken into account.

In Fig. 1, we show the kick velocities of the matched models from numerical simulations performed with a scaled Hobbs kick and the uniform kick distribution. We can see that the flat priors result in preferred kick velocities peaking at the lowest values (below 5 |$\, \text{km s}^{-1}$|⁠) and decreasing steadily – 90 per cent of the distribution is contained below 33 |$\, \text{km s}^{-1}$|⁠. On the other hand, the scaled Hobbs simulations peak between 10 and 20 |$\, \text{km s}^{-1}$|⁠. This is a result of the prior distribution that, with 2000 random kicks, still significantly influences the outcome given the very low weights of small kick velocities (even with the mass scaling). To give a first-order visual representation of the prior distribution in Fig. 1 after mass scaling, we plot a Hobbs distribution scaled by the mean BH remnant mass in the sample on the x-axis, and scaled by an arbitrary constant on the y-axis to make it fit within the frame (the absolute height of the peak is meaningless but the shape is the crucial feature here). We can see that the posterior kick velocities peak ∼20 |$\, \text{km s}^{-1}$| lower than the average prior (a factor of 2 lower), and the best kick velocities from the uniform prior are a factor of ∼30 lower. Consequently, we conclude that VFTS 243-like BH remnants do not receive the same momentum as the ‘typical’ neutron star remnant (assuming that the typical neutron star remnant kick is the Hobbs distribution – see Section 4 for a discussion). This is most probably due to the different nature of the result of the core collapse and explosions yielding the different types of compact remnants.

Figure 1.

Posterior kick distributions of star + BH systems with post-SN periods between 9 and 11 d, and eccentricities lower than 0.05 for a BPASS fiducial population with a metallicity (Z) of 0.006. This is retrieved with 2000 random kicks per SN and remnant masses calculated for weak SNe with E = 10⁵⁰ erg.

3.2 Genealogy of VFTS 243

In addition to the orbital parameters, we can further narrow down matching BPASS systems by using the mass constraints derived by Shenar et al. (2022a). Spectroscopic masses were inferred using cmfgen (19.3 ± 5.2 |$\mathrm{M}_\odot$|⁠; Hillier & Lanz 2001,¹) as well as PoWR (22 ± 5 |$\mathrm{M}_\odot$|⁠; Hamann & Gräfener 2003; Sander et al. 2015)² and fastwind (13 |$\mathrm{\, M}_\odot$|⁠; Puls et al. 2005, 2020).³ To be conservative with their interpretation, they state a mass estimate of 25 ± 12 |$\mathrm{\, M}_\odot$|⁠; later in the analysis when constraints on M_* and M_BH⁴ from the orbital configuration are presented (see their supplementary fig. 10), the masses for the O star are plotted as far as ∼27 |$\mathrm{\, M}_\odot$|⁠. We therefore use a mass window between 13 and 27 |$\mathrm{\, M}_\odot$| for the O star. We also include a total mass constraint of <40 |$\mathrm{\, M}_\odot$| to coincide with the total mass estimate of 36|$^{+3.8}_{-5.4}$|⁠, as well as a mass ratio (q = M_*/M_BH) between 1.5 and 3 (note that because of the mass ratio constraint a low mass limit on the total mass is redundant). Finally, we add a constraint on the radius of the O star to be within roughly 3 |$\mathrm{\, R}_\odot$| of the value derived by Shenar et al. (2022a) of ∼10 |$\mathrm{\, R}_\odot$|⁠. Our search criteria are summarized below and the properties of the two best models are shown in Table 1:

• e (post-SN) < 0.05

• P (days, post-SN) ∈ (9, 11)

• M_* + M_BH < 40 |$\mathrm{\, M}_\odot$|

• M_* ∈ (13, 27) |$\mathrm{\, M}_\odot$|

• q ∈ (1.5–3)

• log(R_*) ∈ (0.85, 1.1)

Table 1.

Properties of the BPASS models that most closely resemble VFTS 243.

Model	M1,i	M2,i	Pi	M1,f	MBH	M*	Pf	e	Mej	N
	(⁠|$\mathrm{\, M}_\odot$|⁠)	(⁠|$\mathrm{\, M}_\odot$|⁠)	(d)	(⁠|$\mathrm{\, M}_\odot$|⁠)	(⁠|$\mathrm{\, M}_\odot$|⁠)	(⁠|$\mathrm{\, M}_\odot$|⁠)	(d)	–	(⁠|$\mathrm{\, M}_\odot$|⁠)	per 106|$\mathrm{\, M}_\odot$|
141048	50	25	15.9	13.8	11.6	25.9	9.2	0.037	0.99	0.8
141238	40	24	15.9	11.6	9.5	24.8	9.0	0.048	1.08	0.9

The masses of the O stars in the best BPASS matches are ∼25–26 |$\mathrm{\, M}_\odot$|⁠, which is the same as was found using BONNSAI (Schneider et al. 2014) by Shenar et al. (2022a). The ages, however, differ: the BONNSAI models predicted an age of 3.9 Myr, whereas our systems are between 4.6 and 5.4 Myr when the primary explodes as an SN. This is likely a direct consequence of the fact that BONNSAI does not take into account the effects of binary interactions. On the other hand, the exemplar evolution calculated by Shenar et al. (2022a) using mesa dies after 7.1 Myr. This is the result of the lower initial masses in their models (M₁ = 30.1 |$\mathrm{\, M}_\odot$| and M₂ = 21.9 |$\mathrm{\, M}_\odot$|⁠) compared to ours (see Table 1) leading to longer lifetimes.

Additionally, the two best models presented here are matched immediately after the first SN, meaning that we only take into account the lifetime of the primary. We can look for secondary models that match our criteria; however, we find that these matches have ages ∼9–10 Myr, which is a little too old for the 30 Doradus region. Our two best models presented here on the other hand are very consistent with the peak of star formation as derived by Schneider et al. (2018).

The most striking difference between the exemplar model presented in Shenar et al. (2022a) and the ones we recover here is the type of mass transfer experienced by the systems. The original VFTS 243 study calls for a binary with an initially very low period (∼3 d) that subsequently undergoes Case A (i.e. on the main sequence) stable mass transfer. Our models on the other hand (see Fig. 2) go through a more typical Case B Roche lobe overflow (RLOF). An episode of common envelope (CE)⁵ occurs after the onset of helium burning. Both the separation and primary mass quickly decrease and the envelope collapses again below the CE threshold. The secondary gains very little mass in the process but the ∼20 solar masses lost by the primary in this interaction lead to a mass reversal. Further mass-loss from the primary occurs when it becomes a helium star that leads to an increase in the period by a few days.

Figure 2.

Evolution of the best progenitor candidates of VFTS 243 in BPASS at Z = 0.006. Note that the time-steps are not evenly spaced in time: the time intervals are dynamically chosen by the code to capture the changing behaviour of the star – in phases where the star changes a lot more steps are taken. This is a convenient way to visualize the changing parameters of a system as age alone can compress the information due to the short time-scales of most critical phases of stellar evolution. Each point on the H–R diagram is a time-step in the model. Note that MH1 is the hydrogen envelope mass of the primary.

Analogues for this system before the first SN exist within our Galaxy. The evolution we see here is similar to that found with earlier BPASS v1.0 models for the binary γ²-Velorum (Eldridge 2009). That system has had a binary interaction in a slightly wider system with some mass transfer having occurred to the O star companion that shows indications of enhanced rotation. While VFTS 243 has a narrower orbit, the highly non-conservative post-main-sequence mass transfer event is observed for both systems.

In the models presented here, the death of the primary is a low-energy (10⁵⁰ erg) SN. Often this is associated with electron-capture SNe, but in this mass range it can be thought of as the output of a partially failed SN. The BH remnant and fast rotation of the progenitor (given the binary interactions it underwent) could lead to the formation of a central engine whose jets, if it is not sufficiently long-lived, may be choked by the envelope. The choked jet model has been put forward as a unifying model for long-gamma-ray-burst SNe and Ic-bl SNe by Lazzati et al. (2012) (with some observational support from spectropolarimetry; Stevance et al. 2017), and it is conceivable that sufficient quenching would lead to lower kinetic energy release, although a quantitative approach is beyond the scope of this study.

In our model search, we tried all three SN energy scenarios available in BPASS as it is unclear what type of energy we should expect a priori from the collapse to a BH. In this case, the weak explosion scenario is the most successful at recreating VFTS 243; hypernova models are excluded completely and typical SNe are strongly disfavoured. Furthermore, given that the ejecta mass for the 10⁵⁰ erg explosions is ∼1 |$\mathrm{\, M}_\odot$|⁠, whereas the ejecta mass constraints in Shenar et al. (2022a) are <0.5 |$\mathrm{\, M}_\odot$|⁠, we can conclude that our models would in fact need to undergo an explosion with E_KE < 10⁵⁰ erg. The fall back of ejected material on to the remnant would increase the remnant masses by a few tenths of |$\mathrm{\, M}_\odot$|⁠; our models would still be fully consistent with the observational criteria and fit a wider range of inclinations for VFTS 243 (see supplementary fig. 10 in Shenar et al. 2022a).

Finally, we show the sample of kick velocities for the two best models highlighted in this section (see Fig. 3). As we can see, the uniform prior samples show a very similar behaviour to that presented in Fig. 1 and very low velocities are preferred. The lack of samples beyond v_kick ∼ 10 |$\, \text{km s}^{-1}$| does not necessarily mean that these velocities are forbidden and could be the result of low-number statistics, as here we have N = 1073 matching models compared to over 8000 in Fig. 1. There are very few matching samples with a Hobbs kick distribution (N = 36) and overall they are outweighed by the uniform kick models by a factor of 30 (when taking into account the incidence rates of each system due to the initial mass function). Overall, considering the physical (mass and radius) constraints on VFTS 243 strengthens the conclusion that this system underwent an extremely low kick that is at odds with the traditional kick velocities assumed for black hole remnants.

Figure 3.

Sampled kick velocities for the two best models (including physical as well as orbital matching criteria).

4 DISCUSSION

To further investigate why the BPASS fiducial models do not naturally result in predictions similar to the system presented in Shenar et al. (2022a), we specifically look for models with similar initial parameters. The discovery paper finds their best model to have M₁ = 30.1 |$\mathrm{\, M}_\odot$|⁠, M₂ = 21.9 |$\mathrm{\, M}_\odot$|⁠, and a period of 3.7 d. On the BPASS model grid, the nearest system is Model 141590, with initial M₁ = 30 |$\mathrm{\, M}_\odot$|⁠, M₂ = 21 |$\mathrm{\, M}_\odot$|⁠, and a period of 2.5 d. Its evolution is shown in Fig. 4.

Figure 4.

Closest BPASS genealogy to the exemplar progenitor model of VFTS 243 presented by Shenar et al. (2022a). Each point on the H–R diagram is a time-step in the model.

Like the mesa model, the BPASS model also undergoes RLOF on the main sequence (Case A), which leads to a significant exchange of mass between the primary and the companion, leading the secondary to reach ∼32 |$\mathrm{\, M}_\odot$|⁠, which is already much greater than the mass estimates from the observations shown in the discovery paper. The primary briefly detaches after the onset of helium burning, then RLOF resumes, and a further ∼5 |$\mathrm{\, M}_\odot$| are accreted. The end product is a 6.7 |$\mathrm{\, M}_\odot$| BH around a 37 |$\mathrm{\, M}_\odot$| star. At this stage, one might argue that future binary interaction could lead to sufficient mass-loss in the secondary star (and mass increase of the BH) to result in a mass regime consistent with VFTS 243; however, Shenar et al. (2022a) demonstrated that the O star in the system is not synchronized with the orbit and therefore no significant interaction has taken place between the companion and the BH.

Overall, the BPASS equivalent to the exemplar model in Shenar et al. (2022a) is not a good progenitor of VFTS 243. Stellar evolution models are complex and there may be several factors at play, but the most likely source of this difference between the BPASS and mesa models is the mass transfer efficiency. The question of how much mass the secondary does accrete during an RLOF episode is not settled and in the mesa models the mass transfer efficiency was treated as a free parameter. Their exemplar model requires 64 per cent of the mass transferred from the primary to be accreted by the secondary. In BPASS, the detailed stellar models are calculated using the Cambridge stars code where the mass transfer efficiency is not a free parameter and in the case of Case A mass transfer, which occurs on the nuclear time-scale, the mass transfer efficiency⁶ is very high (∼97 per cent in this particular case). This explains why the BPASS model with very similar initial parameters results in a system where the O star secondary is too massive to match VFTS 243.

In terms of rate, we find that the models presented in Section 3 are expected to occur 0.8 and 0.9 times (respectively) in a 1 million |$\mathrm{\, M}_\odot$| population. The slightly lower primary mass in the mesa analogue presented here leads to a slightly higher rate: ∼1.2 per million |$\mathrm{\, M}_\odot$|⁠. However, the star formation history of 30 Doradus indicates a higher star formation rate around 5 Myr ago than 7 Myr ago by roughly a factor of ∼1.5 (see Schneider et al. 2018). Overall, occurrence rates alone do not provide clear support to one or the other pathway. A better understanding of stable mass transfer efficiency (particularly on the main sequence) will be required to elucidate which channel is better suited.

Another very uncertain phase in stellar evolution model is the CE. The treatment of CE in BPASS differs from the classical approach (Stevance et al. 2023); one important point is that the CE prescription in our stars implementation does not systematically require the entire envelope to be expelled. This can be seen in the leftmost panels of Fig. 2: >0.5 |$\mathrm{\, M}_\odot$| of hydrogen remains in the envelope when the system detaches. Additionally, we conserve angular momentum rather than binding energy that overall makes our CE more efficient at unbinding the envelope, leading to less orbital shrinkage. A different CE prescription from the one we use would likely lead to very different final systems when using the same initial parameters as those reported in Table 1, and a mismatch with VFTS 243.

The fiducial BPASS populations have been able to reproduce a number of observable phenomena particularly as they relate to massive stars. Nevertheless, CE remains a very uncertain phase of evolution and this should be kept in mind whenever specific evolutionary channels are being studied, as opposed to large populations. In particular, we note here that the CE phase for our favoured models matching VFTS 243 occurs very shortly after the main sequence of the primary. Such Case B interactions in the Hertzsprung gap are uncertain; in rapid population synthesis models, the implementation used is varied as it can have impacts on predictions (e.g. Belczynski et al. 2016). In BPASS, as we can see such a CE leads to significant mass-loss and period decrease but with the orbital radius being mainly unaffected. Most of the orbital evolution occurs before the CE due to mass transfer. This evolution is relatively robust against uncertainties in the CE model, because if CE was to have weaker/stronger orbital evolution than we predict here there are always narrower/wider orbits to provide the same eventual endpoint, respectively.

Another very uncertain aspect of the physics of compact object creation is the natal kicks they undergo. Neutron star kicks have empirical values based on runaway pulsar velocities (e.g. Hobbs et al. 2005; Verbunt, Igoshev & Cator 2017) that are often used in compact object population synthesis. The velocity distribution of a bound neutron star is expected to be lower, however, and there is evidence to support this (e.g. Abbott et al. 2017; Tauris et al. 2017; Vigna-Gómez et al. 2018; Stevance et al., in preparation). As for BH natal kicks, past studies of X-ray binaries have demonstrated that their natal kicks are lower than those of neutron star: between 80 and 310 |$\, \text{km s}^{-1}$| for XTE J1118+480 (Fragos et al. 2009); <80 |$\, \text{km s}^{-1}$| for Cygnus X-1 (Nelemans, Tauris & van den Heuvel 1999); and of the order of a few tens of |$\, \text{km s}^{-1}$| for GRO J1655−40 (Willems et al. 2005). Follow-up analysis by Mandel (2016) demonstrated that these observations did not require natal kicks much greater than 80 |$\, \text{km s}^{-1}$| (although higher kicks could not be ruled out). Due to the lack of BH specific prescription, it is common in modelling to adopt a reduced version of the neutron star kicks (see Section 2).

In this study, we showed that the posterior kick velocity of our matching BPASS models peaks at values significantly lower than one would expect from this reduced kick. From numerical simulations using uniform priors, we find that the BH in VFTS 243 likely received a kick <33 |$\, \text{km s}^{-1}$| with 90 per cent confidence. This is consistent with the finding that our models require a low explosion energy (<10⁵⁰ erg). It is worth considering what change in velocity we might expect the BH to undergo purely from the instantaneous mass-loss. The orbital velocity of the less massive star in a binary system is (following the notation of Blaauw 1961)

$$\begin{eqnarray*} S_2 = 30\times \frac{M_1}{M_1+M_2}\times \sqrt{\frac{M_2}{a} + \frac{M_1}{a}} , \end{eqnarray*}$$

(1)

where M₁ and M₂ are the masses of the most and least massive stellar components, respectively, and a is the separation between the two stars in astronomical units. For our best models, we find that the orbital velocity changes by <2 |$\, \text{km s}^{-1}$| after the death of the primary. This is consistent with the peak of the posterior velocity distribution obtained from our simulations performed with the uniform prior (see Fig. 1), although the tail of this distribution extends past values consistent with this ‘pure Blaauw’ kick. Given the very low eccentricity of VFTS 243, however, it is very likely that the asymmetrical kick velocity was indeed very small.

It is important to note that the Maxwellian kicks and ultra-stripped SN kicks mentioned in this work are not the only prescriptions in the literature. For example, Baker et al. (2008) presented numerical simulations investigating the perpendicular kick as a function of η = [q/(1 + q)]², and Mandel & Müller (2020) developed analytical prescriptions based on the CO core and helium shell masses (which also stand out by their probabilistic nature as opposed to being deterministic like most prescriptions). In another approach presented by Bray & Eldridge (2016), the natal kick values are a function of the mass ejected in the SN explosion. Most recently, Richards et al. (2023) refined the α and β parameters involved by comparing to a larger range of observables. In this prescription, the kick velocity is given by

$$\begin{eqnarray*} v_\text{k} = \alpha \left(\frac{m_\text{ej}}{m_\text{rem}}\right)+\beta \left(\frac{1.4\mathrm{\, M}_\odot }{m_\text{rem}}\right), \end{eqnarray*}$$

(2)

where α and β are constants. We can use their quoted best parameters (⁠|$\alpha =115^{+35}_{-50}\, \text{km s}^{-1}$|⁠; |$\beta = 15^{+10}_{-15}\, \text{km s}^{-1}$|⁠), the ejecta mass (<0.5 |$\mathrm{\, M}_\odot$|⁠), and the remnant mass for VFTS 243 (10.1 ± 2 |$\mathrm{\, M}_\odot$|⁠) to calculate the reduced Bray kick (as opposed to a reduced Hobbs kick; see Section 2). We find a kick velocity of 15.8 ± 2.7 km s⁻¹, which is only a slight overprediction of the peak of the posterior kick distribution obtained from the uniform prior in Section 3 (see Fig. 1). We note that the best parameters obtained by Richards et al. (2023) are calibrated against the Galactic double neutron star system population compiled in Vigna-Gómez et al. (2018) and a catalogue of single-star pulsars from Willcox et al. (2021); the comparison between this prescription and VFTS 243 is encouraging as to its potential to apply the BH population as well, but future research will be required to firmly establish its wider applicability. Future discoveries and studies of quiescent BH binaries will therefore be essential for the further refinement of BH natal kicks, which is crucial as they play an important role in the fate of these systems (Dominik et al. 2012).

5 SUMMARY AND CONCLUSIONS

In this work, we conducted a search through the BPASS v2.2.1 of the recently reported VFTS 243 system containing an O7 star around a quiescent BH in a 10 d orbit at low eccentricity. We investigated three prior kick distributions including the classic Hobbs distributions (Maxwellian with σ = 265 |$\, \text{km s}^{-1}$|⁠) and low-kick uniform distributions (0–100 |$\, \text{km s}^{-1}$|⁠), as well as three SN explosion energies (10⁵⁰, 10⁵¹, and 10⁵² erg). Our key findings are as follows:

• In BPASS, for a metallicity (Z) of 0.006, the most likely progenitor systems have a primary mass of 40–50 |$\mathrm{\, M}_\odot$| with secondaries ∼24–25 |$\mathrm{\, M}_\odot$| in ∼15 d orbits.

• RLOF occurs after the end of the main sequence and the mass transfer is highly non-conservative (accretion efficiency ∼5 per cent).

• The death of the BH progenitor resulted in an explosion with E < 10⁵⁰ erg. This could be a good candidate for a failed SN or, given the evidence for binary interaction and rapid rotation in the O star, a very short lived collapsar leading to choked jets and a low-energy explosion and ejecta mass (∼0.5 |$\mathrm{\, M}_\odot$|⁠; Shenar et al. 2022a).

• The numerical simulations performed with the uniform prior revealed that BH–O star systems with similar orbital parameters as VFTS 243 exhibited low SN kicks (<33 |$\, \text{km s}^{-1}$| in the 90 per cent credible interval) with the peak of the posterior distribution consistent with the magnitude of the Blaauw kick (solely the result of instantaneous mass-loss). This low kick velocity requirement is even stronger in the models that best match VFTS 243, with a strong peak at 1–2 |$\, \text{km s}^{-1}$| and v≤ 10 |$\, \text{km s}^{-1}$|⁠. The very low eccentricity reported for VFTS 243 and very low explosion energy required are also consistent with a scenario with a very low SN kick.

• The Hobbs kick priors are strongly disfavoured, which is consistent with the low SN energy evidence. Indeed, the Hobbs kicks are based on the velocities of runaway pulsars, which are very unlikely to be the result of low-energy explosions.

• The theoretical kick prescriptions based on ejecta masses (Bray & Eldridge 2016; Richards et al. 2023) predict a kick velocity of 15.8 ± 2.7 |$\, \text{km s}^{-1}$|⁠, which is a slight overprediction but similar to peak of the posterior distirbutions obtained from uniform priors.

Further study of quiescent BHs will allow further refinement of the BH natal kick distribution and a better understanding of mass transfer in their progenitor systems. Other candidates found in VFTS (Shenar et al. 2022b) and in Gaia Data Release 3 (Gaia Collaboration 2022) will be the object of a future study, and more observations of such systems are crucial to validating and constraining stellar evolution models at a key milestone in their journey towards becoming BH–BH merger progenitor systems.

ACKNOWLEDGEMENTS

HFS and JJE acknowledge the support of the Marsden Fund Council managed through Royal Society Te Apārangi. SMR, SG, and PT acknowledge support from the University of Auckland. MMB acknowledges the support of the RSNZ.

DATA AVAILABILITY

The BPASS v2.2.1 (and above) data products are available freely – see https://bpass.auckland.ac.nz/. The tui numerical simulations for this work and the codes to perform the analysis and create the plots can be requested by email to the first author. The BPASS python toolkit hoki has been made open source (Stevance et al. 2020). See https://heloises.github.io/hoki/quick_start.html for download instructions.

Footnotes

1

http://kookaburra.phyast.pitt.edu/hillier/web/CMFGEN.htm

2

https://www.astro.physik.uni-potsdam.de/~wrh/PoWR/powrgrid1.php

3

https://fys.kuleuven.be/ster/research-projects/equation-folder/codes-folder/fastwind

4

We use M_* to denote the mass of the bright O star and M_BH to denote the mass of the dark companion.

5

Note that the CE criterion here is R_* > a.

6

Calculated a posteriori.