The enigmatic origin of two dormant BH binaries: Gaia BH1 and Gaia BH2

ABSTRACT

The two systems, namely, Gaia BH1 and Gaia BH2, that have been confirmed as dormant (i.e. no X-ray emission detected) black hole (BH) – low-mass star binaries in the latest Gaia mission data release are intriguing in the context of their formation and evolution. Both systems consist of |$\sim 9\, \mathrm{{\rm M}_{\odot }}$| BH and |$\sim 1\, \mathrm{{\rm M}_{\odot }}$| star orbiting each other on a wide, eccentric orbit (⁠|$e\sim 0.5$|⁠). We argue that formation of such Gaia BH-like systems through the isolated binary evolution (IBE) channel, under the standard common envelope assumptions, and from dynamical interactions in young massive and open clusters are equally probable, and that the formation rate of such binaries is of the order of |$10^{-7}\, \mathrm{{\rm M}_{\odot }}^{-1}$| for both channels. We estimate that, according to our models, there are at most |$\sim 900$| detectable Gaia BH-like binaries in the Milky Way thin disc. What plays an important role in formation of Gaia BH-like systems via the IBE channel is the mutual position of the natal kick velocity vector and the binary angular momentum vector. We find that natal kicks with a median magnitude of |$\sim 40$| km s⁻¹ are preferred for the formation of Gaia BH1-like binaries. Approximately 94 per cent of those binaries are formed with the BH spin misaligned to the orbital axis by less than |$40^{\circ }$|⁠. Gaia BH2-like binaries form if the low-velocity natal kick (of median magnitude |$\sim 20$| km s⁻¹) is directed within |$15^{\circ }$| about the orbital plane. In addition to natal kick, we also discuss the influence of tidal interaction and the adopted common envelope |$\lambda _\mathrm{ce}$| parameter prescription on the evolution of Gaia BH-like binaries. We follow the subsequent evolution of the binaries, once formed as Gaia BH1 and Gaia BH2 systems, to investigate their connection with the low-mass X-ray binary population.

methods: numerical, binaries: general, stars: black holes, stars: evolution, stars: kinematics and dynamics, open clusters and associations: general

1 INTRODUCTION

The dormant black hole (BH) binaries are difficult to observe as they do not emit X-ray radiation and can be discovered only by astrometric and/or spectroscopic analysis. Although such systems were predicted to exist already in the 1960s (Guseinov & Zel’dovich 1966), the new generation of the astronomical instruments and new data analysis methods have allowed to identify the first candidates for such objects only recently. To our knowledge, the first strong candidate for a binary system with a heavy companion unseen in electromagnetic spectrum was found in the globular cluster NGC 3201 (Giesers et al. 2018). Within the last two years, another four systems have been discovered: (i) Binary HD 130 298 where the minimum mass of an unseen companion has been estimated to be |$7.7\pm 1.5$||$\mathrm{{\rm M}_{\odot }}$| (Mahy et al. 2022), (ii) VFTS 243 found in large magellanic cloud which is a binary of a massive star and most probably a black hole companion with mass |$M=10.1\pm 2.0$||$\mathrm{{\rm M}_{\odot }}$| (Shenar et al. 2022), and (iii) two systems which were identified in the most recent Gaia catalogue data release 3 for which there is the strong evidence that they host a non-accreting BH and a low-mass companion star (Chakrabarti et al. 2023; El-Badry et al. 2023a, b; Tanikawa et al. 2023b). In our studies, we focus only on the last two systems and, following El-Badry et al. (2023b), we refer to those binaries as Gaia BH1 and GaiaBH2 (for the observational parameters of the systems see Table 1).

Table 1.

The parameters of Gaia BH1 and Gaia BH2: the orbital period |$P_\mathrm{orb}$| in days, the eccentricity e, BH mass |$M_\mathrm{BH}$|⁠, and companion star mass |$M_{\star }$| in solar units, companion evolutionary stage and age in Gyr, adopted from Chakrabarti et al. (2023); El-Badry et al. (2023a, b). The peculiar velocity at birth |$v_\mathrm{pec}$| is taken from Table 3 in Zhao et al. (2023).

	Gaia BH1	Gaia BH2
|$P_\mathrm{orb}$| (d)	|$185.52\pm 0.08$|	|$1276.7\pm 0.6$|
e	|$0.439^{0.004}_{-0.003}$|	|$0.5176\pm 0.0009$|
|$M_\mathrm{BH}$| (⁠|$\mathrm{{\rm M}_{\odot }}$|⁠)	|$9.326^{0.216}_{-0.209}$|	|$8.94\pm 0.34$|
|$M_{\star }$| (⁠|$\mathrm{{\rm M}_{\odot }}$|⁠)	|$0.93\pm 0.05$|	|$1.07\pm 0.19$|
Star type	MS	RG
Age (Gyr)	|$> 4.0$|	|$8.0-12.0$|
|$v_\mathrm{pec}\, (\mathrm{km\,s^{-1}})$|	|$71.3^{+9.2}_{-10.8}$|	|$34.1^{+5.0}_{-5.1}$|

Once Gaia BH1 and Gaia BH2 had been detected, questions and speculations arose about the origin of those systems. Several studies have already been made which point towards dynamical interactions in young stellar clusters as the most probable scenario to create a BH-low-mass star system having a moderately wide and eccentric orbit (Niccolò Di Carlo et al. 2023; Rastello et al. 2023; El-Badry et al. 2023a; Tanikawa et al. 2023a). However, as we show in this work, the isolated binary evolution (IBE) channel may contribute to the total population of dormant BH binaries at a rate that is comparable to that for the dynamical-interaction channel. The important factors that impact the number of the binaries of interest born through the IBE channel is the magnitude and the geometry of the natal kick imparted to the newborn compact object after the supernova (SN) explosion. One can infer about the natal kick from the peculiar velocity¹ and the orbital parameters of the binary system.

There are two components that change the systemic velocity of a binary after an SN event in the model: (i) the Blaauw kick that impacts the velocity of the binary centre of gravity due to symmetric mass ejection during the SN (Blaauw 1961), (ii) the natal kick that is received by the compact object due to asymmetries of mass ejecta and/or neutrino emission during the SN (e.g. Chugai 1984; Janka 2013, 2017). The observations of proper motions of pulsars and X-ray binaries give constrains on the natal kicks received by neutron stars (NSs) or BHs upon their formation. The Maxwellian distribution with dispersion velocity |$\sigma =265$| km s⁻¹, as derived by Hobbs et al. (2005) from observations of Galactic pulsars, is widely used in the binary population synthesis simulations. This distribution is applied also to BHs but with additional assumption that the kick velocity is reduced by the amount of matter that falls back, being pulled by the gravitational field of the BH forming during the SN. However, recent observations indicate that this approach may underestimate the number of lower velocity systems and overestimate the number of the systems in higher velocity range and that it should be revisited (Atri et al. 2019; O’Doherty et al. 2023; Zhao et al. 2023). The magnitude of the natal kicks is only one part of the puzzle, the other is the kick velocity vector orientation. In the codes based on binary-star evolution (BSE) code (Hurley, Pols & Tout 2000), the direction of a natal kick is drawn from an isotropic spherical distribution. The question is whether the specific direction of the natal kick, that a BH receives in an SN explosion, could lead to the formation of Gaia BH1 and Gaia BH2 (under the standard assumptions of isolated binary evolution model). The reason we may attempt to answer this question for the two systems of interest is that Gaia BH1 and Gaia BH2 are not impacted by the orbit changes due to the mass transfer events after the SN. Gaia BH1 has the orbital separation and eccentricity as it was shaped by SN and Gaia BH2 orbit is only modulated by weak tidal forces.

Our work consists of two main parts: in Section 2, we consider the IBE channel and follow the evolution of the binaries which we classify as Gaia BH1-like and Gaia BH2-like and find the probability of formation of such binaries from the IBE. We then run the evolution of two well-matching Gaia BH1 and Gaia BH2 binaries until the Hubble time and investigate what information Gaia BH1 and Gaia BH2 may bring us about the natal kicks. In Section 3, we present and discuss the formation rates of Gaia BH-like binaries from young massive and open star clusters.

2 GAIA BH BINARIES FROM ISOLATED BINARY EVOLUTION

2.1 Method

We performed the population synthesis calculations using StarTrack. The detailed description of the code can be found in Belczynski, Kalogera & Bulik (2002) and Belczynski et al. (2008, 2016). Here, we highlight the assumptions that are most important in the context of this paper.

The adopted initial mass distribution is a power-law function with three exponents Kroupa, Tout & Gilmore (1993) and Kroupa (2001): |$\alpha _1=-1.3$| for stars of masses |$0.08< M< 0.5\, \mathrm{{\rm M}_{\odot }}$|⁠, |$\alpha _2=-2.2$| for |$0.5\le M< 1.0\, \mathrm{{\rm M}_{\odot }}$| and |$\alpha _3=-2.3$| for stars more massive than |$1.0\, \mathrm{{\rm M}_{\odot }}$|⁠. To maximize the effectiveness of calculations, we focused on the systems that may be the progenitors of Gaia BH binaries: the mass of BH-progenitor star was drawn from the mass range |$18.0-200.0\, \mathrm{{\rm M}_{\odot }}$| and the mass of its companion from the range |$0.2-1.5\, \mathrm{{\rm M}_{\odot }}$|⁠. The orbital period (⁠|$P/\mathrm{orb}$|⁠) (in units of days) and eccentricity e were drawn from the distributions: |$f(\log {P/\mathrm{orb}})=(\log {P/\mathrm{orb}})^{-0.55}$| with |$\log {P/\mathrm{orb}}$| in the range [0.15, 4.0] and |$f(e)=e^{-0.42}$| with e in the range [0.0, 0.9], respectively (Sana et al. 2012, 2013).

In our calculations, the common envelope phase is parametrized by ce parameter |$\alpha _\mathrm{{\small CE}}=\alpha _{\mathrm{ce}}\times \lambda _{\mathrm{ce}}$|⁠, where |$\lambda {\mathrm{ce}}$| is a parameter describing the donor central concentration (de Kool 1989) and |$\alpha _{\mathrm{ce}}$| parametrizes the efficiency at which the orbital energy is transferred into the envelope (Webbink 1984). We assume the standard value |$\alpha _{\mathrm{ce}}=1.0$| and for |$\lambda _{\mathrm{ce}}$| we use the fits from Xu & Li (2010a, b) and Wang, Jia & Li (2016) (see Section 2.3.6 for discussion).

The massive stars (of type O/B) stellar winds were described by the equations for radiation driven mass-loss from (Vink, de Koter & Lamers 2001) which take into account the luminous blue variable mass-loss (Belczynski et al. 2010). For Wolf–Rayet stars we took into account the wind clumping (Hamann & Koesterke 1998) and the metallicity dependence of the winds (Vink & de Koter 2005). The winds of less massive stars were calculated according to the formulae of Hurley et al. (2000).

We assumed the delayed supernova engine (Fryer et al. 2012) in our calculations. There are two components that make up the final kick received by the system in SN explosion: the kick from asymmetric neutrino and/or ejecta emission (natal kick) and the Blaauw kick resulting from the symmetric envelope ejection from the exploding star. We tested two distributions of the natal kick velocities: (i) the widely used Maxwellian distribution with |$\sigma =265\, \mathrm{km\, s^{-1}}$| (Hobbs et al. 2005) which for BH is decreased by the fallback of the fraction of the stellar envelope that is initially ejected in SN explosion, this fraction is defined by the fallback parameter |$f_{\mathrm{fb}}$| which takes the values between 0 and 1 (model V1), (ii) the two-component Maxwellian with |$\sigma _1=21.3\, \mathrm{km\, s^{-1}}$| for low-velocity component and |$\sigma _2=106.7\, \mathrm{km\, s^{-1}}$| for high-velocity component derived from the sample of 85 BH and NS X-ray binaries by Zhao et al. (2023) also decreased by the fallback for BHs (model V2). The directions of the kicks are drawn from isotropic distribution. For each system the drawing of natal kick magnitude and orientation was repeated 200 times.

The systemic velocity |$v_\mathrm{sys}$| calculated in StarTrack is the velocity of the binary centre of mass resulting from Blaauw kick and natal kick. It is assumed that the centre of mass is initially at rest. |$v_\mathrm{sys}$| is consistent with the potential peculiar velocity at birth |$v_\mathrm{p,z=0}$| (Atri et al. 2019), that is the velocity that the binary gained from SN which caused its motion out of the Galactic plane. The potential peculiar velocity is inferred from the observational present-day peculiar velocity (for details see Atri et al. 2019 and Zhao et al. 2023).

We position the orbital plane of a binary before the first SN on x–y plane. All three vectors, the orbital angular momentum and the spins of a star and BH, are aligned and pointing in z direction. Setting binary in such position does not influence our results while it makes the analysis of the post-SN relative orientation of those vectors much easier. The position of the exploding star on the orbit is found by means of mean anomaly M, eccentric anomaly E, eccentricity e, and the time of periapsis passage |$\tau$|⁠.

Taking into account the age constrains derived from the observations of the luminous component in Gaia BH1 and Gaia BH2 we assumed that both binaries are |$\sim 4$| Gyr old (El-Badry et al. 2023a, b) and, therefore, were formed in the environment of metallicity |$\mathit{ Z}=0.7Z_{\odot }$| (see fig. 3 in Olejak et al. 2020) where |$\mathit{ Z}_{\odot }=0.014$| (Asplund et al. 2009).

It is hard to expect that one may get many systems of the parameters exactly corresponding to those inferred from observations of Gaia BH1 and Gaia BH2 in the simulations. Therefore, by putting constrains on the orbital separation a, the eccentricity e, the BH mass |$M_\mathrm{BH}$|⁠, the companion star mass |$M_2$|⁠, and the companion star evolutionary stage, we defined three groups of binaries which are of our interest: Gaia BH1-like, Gaia BH2-like, and Gaia BH-like. Gaia BH1- and Gaia BH2-like are the systems assumed to correspond to observed Gaia BH1 and Gaia BH2 while Gaia BH-like is the general group of all systems which may be classified as dormant BH binaries. To investigate the interconnection between the natal kicks and the post-SN quantities: a, e, |$v_\mathrm{sys}$|⁠, orbital angular momentum, and BH spin vectors orientation we defined the fourth group – the ‘optimistic’ Gaia BH-like binaries – for which we allow the wider range of a, e, and both components masses. The classification criteria for each group are given in Table 2.

Table 2.

The parameters of the binaries that we classify as Gaia BH1-, Gaia BH2-, and Gaia BH-like in our standard models V1 and V2. We take into account both components masses (⁠|$M_\mathrm{BH}$| and |$M_2$|⁠), the evolutionary stage of the companion [main-sequence (MS)/red giant (RG))] the orbital separation a, the eccentricity e, and the systemic velocity |$v_\mathrm{sys}$| (km s⁻¹).

	Gaia BH1-like	Gaia BH2-like	Standard Gaia BH-like	Optimistic Gaia BH-like
|$M_\mathrm{BH}$| (⁠|$\mathrm{{\rm M}_{\odot }}$|⁠)	8.0–12.0	7.0–11.0	<15.0	<15.0
|$M_2$| (⁠|$\mathrm{{\rm M}_{\odot }}$|⁠)	0.8–1.2	0.8–1.2	0.8–1.2	<1.5
Companion star	MS	RG	MS/RG	MS/RG
a (⁠|$\mathrm{R_{\odot }}$|⁠)	200–400	950–1150	200–1200	100–1500
e	0.4–0.6	0.4–0.6	0.3–0.7	0.1–0.9
|$v_\mathrm{sys}\, [\mathrm{km\,s^{-1}}]$|	no limit	no limit	no limit	no limit

We define two angles which help to complement our analysis: (i) |$\gamma$| - a misalignment between the binary angular momentum vector and the natal kick velocity vector, and (ii) |$\theta$| - a misalignment between the binary angular momentum vector and post-SN BH spin vector.

2.2 Results

2.2.1 The formation rate of Gaia BH-like systems from IBE

The formation rate of Gaia BH-like binaries depends on the criteria according to which one classifies the binary as Gaia BH-like and on the adopted natal kick velocity distribution and it may differ as much as three orders of magnitude as can be seen from Table 3. The formation rate for model V1 is an order of magnitude smaller than formation rates for model V2. There is also one order-of-magnitude difference between the standard and optimistic groups in both models: for optimistic and standard Gaia BH-like groups in model V2 the formation rate is |$\sim 1.4\times 10^{-7}\, \mathrm{{\rm M}_{\odot }}^{-1}$| and |$\sim 1.0\times 10^{-8}\, \mathrm{{\rm M}_{\odot }}^{-1}$|⁠, respectively and in model V1 it is |$\sim 1.6\times 10^{-9}\, \mathrm{{\rm M}_{\odot }}^{-1}$| and |$\sim 1.1\times 10^{-8}\, \mathrm{{\rm M}_{\odot }}^{-1}$|⁠, respectively.

Table 3.

The formation rates of Gaia BH-like binaries for different classification criteria and natal kick velocity distributions V1 and V2. The criteria adopted to calculated the formation rate for case ‘V2 low-mass comp’ correspond to the population of BH-star binaries with low-mass companion (⁠|$M_{\rm comp}< 3\, \mathrm{{\rm M}_{\odot }}$|⁠) from young massive and open clusters described in Section 3.2. We present the formation rates of Gaia BH1- and Gaia BH2-like binaries in model V1 and V2 in the last three rows.

Case	Rate (⁠|$\mathrm{{\rm M}_{\odot }}^{-1}$|⁠)
V1 standard	|$\sim 1.6\times 10^{-9}$|
V1 optimistic	|$\sim 1.1\times 10^{-8}$|
V2 standard	|$\sim 1.0\times 10^{-8}$|
V2 optimistic	|$\sim 1.4\times 10^{-7}$|
V2 low-mass comp	|$\sim 2.2\times 10^{-6}$|
V1 standard	|$\sim 1.6\times 10^{-9}$|
V1 |$\&$| V2 Gaia BH1	|$\sim 3.2\times 10^{-10}$|
V1 Gaia BH2	|$\sim 1.4\times 10^{-10}$|
V2 Gaia BH2	|$\sim 9.6\times 10^{-10}$|

The formation rates of Gaia BH1- and Gaia BH2-like systems, which are the two binary categories with the most stringent classification criteria among the binary types defined in Table 2, are on order of |$10^{-10}\, \mathrm{{\rm M}_{\odot }}^{-1}$|⁠. The formation rate of Gaia BH2-like systems is almost an order of magnitude higher in model V2 than in model V1, while the formation rate of Gaia BH1-like binaries does not depend on the natal kick distribution adopted (model V1 versus V2).

Additionally, to make the comparison with the results of dynamical channel described in Section 3.2, we calculated the general population of BH-star binaries with the only constrain imposed on the companion star mass (⁠|$M_{\rm comp}< 3.0\, \mathrm{{\rm M}_{\odot }}$|⁠) and we found that the formation rate of such binaries from IBE channel is |$2.2\times 10^{-6}\, \mathrm{{\rm M}_{\odot }}^{-1}$|⁠.

Among Gaia BH1-like an Gaia BH2-like binaries, we present two which are our closets match to the observational systems in Table 4.

Table 4.

The parameters of the two Gaia BH1-like and Gaia BH2-like binaries that we find closely matching the observed systems: the orbital period |$P_\mathrm{orb}$| in days, the eccentricity e, BH mass |$M_\mathrm{BH}$|⁠, and companion star mass |$M_{\star }$| in solar units and the systemic velocity |$v_\mathrm{sys}$|⁠.

	Gaia BH1-like	Gaia BH2-like
|$P_\mathrm{orb}$| (d)	181.07	1591.63
e	0.48	0.56
|$M_\mathrm{BH}$| (⁠|$\mathrm{{\rm M}_{\odot }}$|⁠)	8.22	8.55
|$M_{\star }$| (⁠|$\mathrm{{\rm M}_{\odot }}$|⁠)	1.38	1.14
|$v_\mathrm{pec}\, [\mathrm{km\,s^{-1}}]$|	57.73	35.57

2.2.2 The evolution of Gaia BH1- and Gaia BH2-like systems

The conditions that we set on BH masses for Gaia BH1- and Gaia BH2-like binaries limit the zero-age main sequence (ZAMS) masses of BH progenitors to the range |$M_\mathrm{ZAMS}\sim 34.0-40.0\, \mathrm{{\rm M}_{\odot }}$| (for |$\mathit{ Z}=0.01$|⁠). We find that the initial orbital separations may reach up to |$\sim 36\, 000\, \mathrm{R_{\odot }}$| for very high-initial eccentricities (⁠|$e\sim 0.9$|⁠). For the moderate initial eccentricities |$e< 0.4$| the initial orbital separations are in the range |$\sim 6000-9000\, \mathrm{R}_{\odot }$|⁠.

Whatever their starting parameters are, all Gaia BH1-like binaries follow the same evolutionary path which we describe on one example also shown on Fig. 1. The evolution starts with two main-sequence stars with masses |$39.7$| and |$1.2\, \mathrm{{\rm M}_{\odot }}$| on a wide (⁠|$a\sim 7\, 941\, \mathrm{R_{\odot }}$|⁠), moderately eccentric (⁠|$e=0.31$|⁠) orbit. In the course of the nuclear evolution of the more massive star (the primary) the orbit widens due to the strong stellar winds and it reaches |$a\sim 8\, 542\, \mathrm{R_{\odot }}$| at the time when the star leaves the Hertzsprung gap and starts to burn helium in its core. The still expanding radius of the core helium burning star is large enough (⁠|$\sim 1550\, \mathrm{R_{\odot }}$|⁠) that the tides start to act between the two stars causing the orbit to gradually tighten and circularize. By the time the orbit becomes circular the orbital separation shrinks to |$a\sim 5\, 841\, \mathrm{R_{\odot }}$|⁠. But the tidal interactions are still at work trying to speed up the slowly rotating massive star and synchronize it with its companion which is spinning three orders of magnitude faster than the giant. The synchronization cannot be reached and Darwin instability sets in leading to further orbital decay until the giant star overflows its Roche lobe. Up to this moment the evolved star has lost |$\sim 57~{{\ \rm per\ cent}}$| of its initial mass (its mass is now |$17.2\, \mathrm{{\rm M}_{\odot }}$|⁠), the mass of MS star has remained almost unchanged and the orbital separation between the stars decreased to |$a\sim 3\, 256\, \mathrm{R_{\odot }}$|⁠. The mass transfer that follows is dynamically unstable even if one considers the more restrictive conditions for unstable mass transfer than in our standard model (see Pavlovskii et al. 2017; Olejak & Belczynski 2021 for details). At this point the evolved donor is a red supergiant with the convective envelope and low-effective temperature |$T_\mathrm{eff}\sim 3.55\times 10^3$| K. The sufficiently low-binding energy of such an envelope leads to its successful ejection during the common envelope (ce) phase.

Figure 1.

The evolution of Gaia BH1-like binary from the isolated binary evolution channel. The chosen example has the parameters close to the observational ones.

The giant star loses its envelope becoming a Wolf–Rayet star on a close orbit (⁠|$a\sim 140\, \mathrm{R_{\odot }}$|⁠) with MS star. In the subsequent evolution that lasts |$\sim 0.22$| Myr, the strong winds drive further mass-loss from the primary star and the orbit expands to |$a\sim 164\, \mathrm{R_{\odot }}$| just prior the SN explosion. The SN explodes |$\sim 5.7$| Myr from the onset of the binary evolution. The primary star becomes |$9.7\, \mathrm{{\rm M}_{\odot }}$| BH and the effective kick (the composition of Blaauw and natal kicks) imparted to the system increases both its orbital separation and the eccentricity to |$a\sim 319\, \mathrm{R_{\odot }}$| and |$e=0.49$| leading to the formation of Gaia BH1-like binary.

Our synthetic Gaia BH1 binary evolution proceeds on the nuclear time-scale of the secondary which slowly evolves towards the end of hydrogen burning in its core. It takes |$\sim 5$| Gyr for the star to leave the main sequence and become a red giant what marks the end of Gaia BH1 stage of the binary evolution. As the radius of the red giant expands the tides act to circularize the orbit. About |$\sim 600$| Myr later the red giant radius reaches the radius of the star Roche lobe and the mass transfer on to BH begins. The system turns into low-mass X-ray binary (LMXB) with the orbital period |$P_\mathrm{orb}\sim 140$| d (⁠|$a\sim 251\, \mathrm{R_{\odot }}$|⁠) and donor mass |$M_2=1.2\, \mathrm{{\rm M}_{\odot }}$|⁠. The mass transfer proceeds from less to more massive binary component, therefore, the orbit widens and so does the radius of the secondary star Roche lobe. When the Roche lobe extends to the point that the star fits inside it again, the mass transfer ceases. It happens when the orbital separation is |$a\sim 892$| |$\mathrm{R_{\odot }}$|⁠, |$M_\mathrm{BH}\sim 10.2\, \mathrm{{\rm M}_{\odot }}$|⁠, and |$M_2\sim 0.6\, \mathrm{{\rm M}_{\odot }}$|⁠. The LMXB phase lasts for |$\sim 4$| Myr. From there, the secondary star continues its nuclear evolution until it transforms into a carbon–oxygen white dwarf (WD) of mass |$M_\mathrm{WD}=0.5\, \mathrm{{\rm M}_{\odot }}$|⁠. At the Hubble time |$t_\mathrm{Hub}=13.7$| Gyr, the distance between the two compact objects in the system is |$a\sim 1244\, \mathrm{R_{\odot }}$|⁠. For the rest of its life BH–WD binary slowly spirals in due to the gravitational radiation. The system will not merge within Hubble time.

To illustrate the evolution of Gaia BH2 binary, we choose the binary which starts as the system consisting of massive MS star of mass |$22.6\, \mathrm{{\rm M}_{\odot }}$| and MS star of mass |$1.1\, \mathrm{{\rm M}_{\odot }}$| orbiting each other on the moderately eccentric (⁠|$e=0.66$|⁠) wide orbit (⁠|$a\sim 7\, 475\, \mathrm{R_{\odot }}$|⁠) (see Fig. 2). As in the case of Gaia BH1, the system first evolves on the more massive (primary) star nuclear evolution time. At first the orbital separation increases due to the stellar wind from the primary but when the primary becomes a core helium burning giant and its radius reaches |$R_2\sim 520.4\, \mathrm{R_{\odot }}$| the tidal force starts to control the orbital evolution. At this point the orbital separation (⁠|$a\sim 7\, 795\, \mathrm{R_{\odot }}$|⁠) starts to shrink. It takes |$\sim 0.32$| Myr for the tidal force to circularize the orbit. While the orbit keeps shrinking the primary evolves on to the early asymptotic giant branch (AGB). The mass of the primary just before it explodes as SN has decreased to |$M_1=18.5\, \mathrm{{\rm M}_{\odot }}$| and the orbital separation has contracted by over a half of its initial size, down to |$a\sim 3043\, \mathrm{R_{\odot }}$|⁠. The massive star becomes a |$9.2\, \mathrm{{\rm M}_{\odot }}$| BH after |$\sim 9.53$| Myr from the start of the evolution. The large amount of mass (⁠|$\sim 50~{{\ \rm per\ cent}}$|⁠) that it loses in the explosion results in high eccentricity (⁠|$e=0.77$|⁠) that the orbit gains mainly due to the Blaauw kick but at the same time the systemic velocity remains low |$v_{\mathrm{sys}}\sim 32.8\, \mathrm{km\, s^{-1}}$| due to the low-natal kick velocity (⁠|$v_{\mathrm{kick}}\sim 34.8\, \mathrm{km\, s^{-1}}$|⁠).

Figure 2.

The main stages of Gaia BH2-like binary evolution from the isolated binary evolution channel. The details are described in Section 2.2.2.

The system consisting of BH–RG forms |$\sim 6.2$| Gyr after the onset of the binary evolution. While RG expands the tides act to tighten and circularize the orbit. When the orbital separation and the eccentricity decrease to |$a\sim 1175\, \mathrm{R_{\odot }}$| and |$e=0.59$|⁠, respectively, they fall into the ranges which define the binary as Gaia BH2-like (see Table 2). This phase of the binary evolution lasts only for |$\sim 2$| Myr. In next |$\sim 4$| Myr the eccentricity drops to zero and the orbital separation reaches its minimum, i.e. |$a\sim 758\, \mathrm{R_{\odot }}$|⁠. At this point the wind from the expanding RG comes into play: the orbit widens again, and BH accrets part of the matter lost by the evolved companion. By the time He is ignited in the core (⁠|$\sim 6.7$| Myr after orbit circularization) the star mass drops to |$M_2=0.75\, \mathrm{{\rm M}_{\odot }}$|⁠, BH mass slightly increases to |$M_{\mathrm{BH}}=9.4\, \mathrm{{\rm M}_{\odot }}$| and orbit widens to |$a\sim 1165\, \mathrm{R_{\odot }}$|⁠. The exhaustion of He in the core occurs after |$\sim 112$| Myr when the companion becomes an early AGB star and after next |$\sim 5.7$| Myr it evolves to become the thermally pulsing AGB star. The AGB stage, when the star loses roughly 25 per cent of its mass and the orbit widens by almost a factor of two reaching |$a\sim 2217\, \mathrm{R_{\odot }}$|⁠, is followed by the formation of carbon oxygen WD (CO WD) of mass |$M_\mathrm{WD}=0.57\, \mathrm{{\rm M}_{\odot }}$| orbiting BH of mass |$M_\mathrm{BH}=9.5\, \mathrm{{\rm M}_{\odot }}$| at a distance |$a\sim 2247\, \mathrm{R_{\odot }}$|⁠.

2.2.3 Systemic velocities and natal kicks

We present |$v_\mathrm{sys}$| distributions of standard Gaia BH-like binaries in models V1 (Hobbs et al. 2005; natal kick velocity distribution) and V2 (Zhao et al. 2023; natal kick velocity distribution) in green in Fig. 3. On top of them, we show the distribution of |$v_\mathrm{sys}$| for subsets of binaries which we identify as Gaia BH1-like (magenta) and Gaia BH2-like (blue), the systems which have |$v_\mathrm{sys}$| exactly in the observational range for those systems are marked in pattern. The observed |$v_\mathrm{sys}$| of Gaia BH1 and Gaia BH2 are within the ranges predicted by isolated binary evolution for Gaia BH1-like and Gaia BH2-like binaries.

Figure 3.

The distribution of the systemic velocity of all standard Gaia BH-like binaries (the dominant distribution in green) in model V1 (top) and model V2 (bottom). The magenta distribution stands for the subset of Gaia BH1-like binaries and the distribution in blue (the least numerous) is for Gaia BH2-like binaries. The binaries in both groups which have |$v_\mathrm{sys}$| exactly in the observational ranges are shown as boxes with pattern. The criteria that we adopt for each Gaia BH-like group are listed in Table 2.

For all standard Gaia BH-like binaries the range of achievable |$v_\mathrm{sys}$| is wide (up to |$\sim 300$| km s⁻¹ in model V1 and |$\sim 200$| km s⁻¹ in model V2) but the median values are well below 100 km s⁻¹: it is |$v_\mathrm{sys,med,V1}=72.8$| km s⁻¹ for all Gaia BH-like binaries in model V1 and |$v_\mathrm{sys,med,V2}=45.3$| km s⁻¹ in model V2. For Gaia BH1-like binaries |$v_\mathrm{sys,med,V1}\sim 54.8$| km s⁻¹ in model V1 and |$v_\mathrm{sys,med,V2}\sim 35.4$| km s⁻¹ in model V2, for Gaia BH2-like binaries |$v_\mathrm{sys,med,V1}\sim 19.2$| km s⁻¹ in model V1 and |$v_\mathrm{med,sys,V2}\sim 16.7$| km s⁻¹ in model V2. For Gaia BH2-like binaries the median systemic velocity is comparable between the two natal kick velocity distributions while Gaia BH1-like binaries move on average by |$\sim 55~{{\ \rm per\ cent}}$| faster in model V1 than in model V2.

Overall there is 25 per cent more Gaia BH-like binaries in model V2 than in model V1. Considering the two specific subsets of binaries it appears that the number of Gaia BH1-like systems is almost the same in both models (⁠|$\pm 3$| systems) but Gaia BH2 are almost seven times more abundant in model V2. It is a straightforward consequence of the lower mean velocity of the natal kick velocities distribution of Zhao et al. (2023) and low velocities of natal kicks allow the wide and eccentric binaries like Gaia BH2 to survive SN.

The corresponding median magnitude of the natal kick is not very different from |$v_\mathrm{sys}$| due to the low mass of the companion star (⁠|$\sim 1\, \mathrm{{\rm M}_{\odot }}$|⁠) and for standard Gaia BH-like binaries it is |$v_\mathrm{kick,med,V1}=82.9$| km s⁻¹ in model V1 and |$v_\mathrm{kick,med,V2}=51.0$| km s⁻¹ in model V2.

The distributions of the natal kick velocities received by Gaia BH1-like and Gaia BH2-like binaries for model V2 are shown on top and bottom panels of Fig. 4 correspondingly. The median natal kick magnitude is |$v_\mathrm{kick}\sim 39.3$| km s⁻¹ for Gaia BH1-like systems and |$v_\mathrm{kick}\sim 18.7$| km s⁻¹ for Gaia BH2-like systems. We mark the distributions of natal kick velocities for Gaia BH1- and Gaia BH2-like binaries which have the systemic velocities exactly in the observational limits (Zhao et al. 2023) in dark-blue.

Figure 4.

The distribution of the natal kick velocity received by the synthetic Gaia BH1-like (top) and Gaia BH2-like (bottom) in model V2. The dark-blue distribution stands for the subset of the binaries that have |$v_\mathrm{sys}$| in the range corresponding to the observational estimations: |$v_\mathrm{sys}\in [60.0-80.0]$| km s⁻¹ for Gaia BH1-like and |$v_\mathrm{sys}\in [30.0-40.0]$| km s⁻¹ for Gaia BH2-like.

In our further analysis in Section 2.3.4, we focus on model V2 as it is well motivated by observations, it gives us better statistics of Gaia BH2 binaries and the conclusions of Section 2.3.4 do not depend on the choice of |$v_\mathrm{kick}$| distribution.

2.3 Discussion

2.3.1 The formation rate of Gaia BH-like binaries

The total stellar mass of one IBE model is |$M_\mathrm{sim}=3.917\times 10^{11}\, \mathrm{{\rm M}_{\odot }}$|⁠. It comes from our choice to evolve only the binaries which may be the progenitors of Gaia BH-like binaries, that is we take primary ZAMS mass to be |$M_\mathrm{1,ZAMS}> 18\, \mathrm{{\rm M}_{\odot }}$| (so that it becomes BH) and the second star ZAMS mass to be |$M_\mathrm{2,ZAMS}< 1.5 \mathrm{{\rm M}_{\odot }}$|⁠. We then evolve |$N=10^6$| such binaries. To calculate |$M_\mathrm{sim}$| we draw the binaries from the whole initial mass function (IMF) and sum up their masses and we count the number of the subset of drawn binaries that have |$M_\mathrm{1,ZAMS}> 18 \mathrm{{\rm M}_{\odot }}$| and |$M_\mathrm{2,ZAMS}< 1.5 \mathrm{{\rm M}_{\odot }}$|⁠. We draw the binaries (summing up their masses) until we have a subset of N (⁠|$10^6$|⁠) binaries. We also take into account the median number of natal kick draws that are needed to get a Gaia BH-like binary. Therefore, the calculated formation rates are corrected for using truncated IMF for both of the binary members and for the number of kick draws.

If we scale down IBE simulations to match the total simulated initial cluster mass of star clusters simulations |$M_{\rm cl,tot}$||$=3.6\times 10^6\, \mathrm{{\rm M}_{\odot }}$| (see Section 3.2), we find one Gaia BH-like binary per simulation. This corresponds to the result of the dynamical interaction channel where also one Gaia BH-like binary ejected from the cluster was found in a given model cluster set. In conclusion, both channels give the same formation rate that is |$\sim 2.7\times 10^{-7}\, \mathrm{{\rm M}_{\odot }}^{-1}$|⁠.

2.3.2 The estimated number of Gaia BH-like binaries

A considerable amount of work has so far been done to estimate the number of dormant BHs that can be present in the Milky Way within Gaia detection capability (e.g. Breivik, Chatterjee & Larson 2017; Mashian & Loeb 2017; Yamaguchi et al. 2018; Wiktorowicz et al. 2020; Chawla et al. 2022; Shikauchi, Tanikawa & Kawanaka 2022; Shikauchi et al. 2023). Typically, the resulting number of binaries of interest varies from few tens to few hundreds (e.g. Yamaguchi et al. 2018; Wiktorowicz et al. 2020; Chawla et al. 2022). However, the number can reach as high as |$10^3-\sim 10^5$| (e.g. Breivik et al. 2017; Mashian & Loeb 2017) depending on the physics adopted in the evolutionary models, on the accuracy with which the Milky Way is modelled (where for example the realistic stellar distribution, e.g. Wiktorowicz et al. 2020; Chawla et al. 2022), and/or how the position-dependent reddening (Chawla et al. 2022) is taken into account.

Although the main focus of our paper is the comparison of the formation rates of Gaia BH-like binaries for the isolated-binary and the dynamical formation channels, we complement this work with an approximate estimation of the expected number of Gaia BH-like binaries in the Milky Way within the scope of Gaia detectability, according to the two models. We limit our considerations to the Galactic thin disc population, since our simulations of the isolated binary evolution are done for a metallicity that is close to solar.

As an intrinsic synthetic population of Gaia BH-like binaries, we choose that from the model ‘V2 low-mass comp’ (see 2.2.1), where we assume that the natal kicks distribution is a double peaked Maxwellian and require the mass of the companion to be smaller than |$3.0\, \mathrm{{\rm M}_{\odot }}$| for the binary to be classified as Gaia BH-like. We calibrate the mass of the simulations to match the stellar mass of the Galactic thin disc, where, for the latter, we adopt |$M_\mathrm{d,thin}=2.15\times 10^{10}\, \mathrm{{\rm M}_{\odot }}$| following Breivik et al. (2017) who assumed a constant star formation rate of |$2.15\, \mathrm{{\rm M}_{\odot }}\mathrm{yr}^{-1}$| over 10 Gyr. We draw the present-day population of BH–low-mass star binaries from this intrinsic binary population, taking into account the estimated age of Gaia BH1 (⁠|$\sim 8$| Gyr) and the estimated duration of the red giant phase of the companion star in Gaia BH2 (⁠|$\sim 100$| Myr).

In the next step, the resulting binary population is spatially distributed at random, following the thin disc mass distribution described by the function from Yu & Jeffery (2015), i.e.

$$\begin{eqnarray} \rho (R,z)=\frac{M_\mathrm{d,thin}}{4\pi h_\mathit{ R}\, h_z}e^{R/h_\mathit{ R}}\, \mathrm{ sech}^2(-z/h_z), \end{eqnarray}$$

(1)

where |$h_R=2.5$| kpc is the scale length of the disc, |$h_z=0.352$| kpc is its scale height, and R and z are the cylindrical coordinates of the axisymmetric disc.

To calculate the apparent magnitude |$m_V$| of all considered present-day Gaia BH-like binaries, we first calculate the absolute magnitude and the bolometric correction (Torres 2010) for each binary system in the sample using the luminosity, |$L_2$|⁠, and the effective temperature, |$T_\mathrm{eff,2}$|⁠, of the companion star. For the interstellar extinction in the V-band, |$A_V$|⁠, we follow the approach of Yamaguchi et al. (2018) who assumed that |$A_V(d)=d$|⁠. Here, d in the distance from Earth in kpc that we calculate from the binary’s position in the thin disc. We set a lower limit of |$m_V< 20$| mag for the apparent magnitude of the binary, if it is to be detected by Gaia.

The astrometric signature, |$\alpha$|⁠, is calculated following Breivik et al. (2017):

$$\begin{eqnarray} \alpha =\biggl (\frac{a_\mathrm{proj}}{\mathrm{au}}\biggr)\biggl (\frac{d}{\mathrm{pc}}\biggr)^{-1}\, \, \, \mathrm{arcsec}, \end{eqnarray}$$

(2)

where d is the distance to the binary and |$a_\mathrm{proj}$| is the projected semimajor axis of the luminous companion.

We calculate |$a_\mathrm{proj}$| for each system, using the exact position of the star on its orbit around the compact object and the orbital elements i, |$\Omega$|⁠, and |$\omega$| at a randomly chosen time elapsed from the star’s periastron passage. We consider the binary as detectable if |$\alpha > 3\sigma _\mathrm{G}$|⁠, where |$\sigma _\mathrm{G}$| is Gaia’s astrometric precision calculated following Gaia Collaboration (2016) (see equation 4 therein).

The final estimated number of BH–low-mass star binaries expected in the Galactic thin disc from the model of isolated binary evolution, assuming the most optimistic classification criteria, is of the order of |$1.4\times 10^3$| and, among them, |$\sim 900$| binaries are within Gaia’s detectability limit. In Fig. 5, we present the distributions of the orbital periods and eccentricities of those binaries. It can be seen that binaries with orbital period shorter than 100.0 d (⁠|$\sim 66~{{\ \rm per\ cent}}$| of all binaries) and eccentricity larger than 0.5 (⁠|$\sim 68~{{\ \rm per\ cent}}$| of all binaries) dominate the distributions. About half of the binaries (⁠|$\sim 44~{{\ \rm per\ cent}}$|⁠) have both |$P_{\rm orb}$| and e in these ranges.

Figure 5.

The distribution of the orbital periods |$P_\mathrm{orb}$| (top) and eccentricity e (bottom) of the predicted population of Gaia BH-like binaries evolved through IBE detectable by Gaia during 5 yr observational period.

2.3.3 The evolution Gaia BH-like binaries

In the following paragraphs we focus on the three phases of the binary evolution that we find important to discuss. The thorough studies of the impact of the choice of SN model (rapid versus delayed, with or without fallback), or of changing the value of the common envelope efficiency parameter |$\alpha _\mathrm{ce}$| on the formation of detached BH–star binaries can be found for example in Chawla et al. (2022) and Shikauchi et al. (2023). Although our results differ from conclusions of Shikauchi et al. (2023) on |$\alpha _\mathrm{ce}$| (see Section 2.3.6) it is instructive to see the influence of |$\alpha _\mathrm{ce}> 1$| on the dormant BH populations.

2.3.3.1 The evolution Gaia BH-like binaries: the impact of the natal kicks

The natal kick velocity |$v_\mathrm{kick}$| is a vector that plays an important if not the decisive role in the Gaia BH-like systems formation. What counts is not only the magnitude of the kick but also its direction and relative position of the stars on their mutual orbit when SN takes place.

The natal kick of a given magnitude determines |$v_\mathrm{sys}$| but not the eccentricity of the binary, what is demonstrated in Fig. 6. However, there is a correlation between the eccentricity, post- and pre-SN orbital separation (see e.g. Kalogera 1996) which is followed by our population of Gaia BH-like binaries as is shown in Fig. 7. For Gaia BH1, where the companion star is still on its main sequence, the evolution has not have enough time to change the binary orbit from its immediate post-SN configuration. The system may have widen or shrunken due to the natal kick but the total amount of the orbital separation change is tightly connected to the system post-SN eccentricity. If the natal kick causes binary orbit to widen the orbital separation may increase no more than twice (⁠|$a_\mathrm{post-SN}/a_\mathrm{pre-SN}\sim 2$|⁠) for binaries that have eccentricities |$e\in [0.4;0.6]$|⁠. If the orbit shrinks due to the natal kick then it can shrink by no more than 40 per cent (⁠|$0.6< a_\mathrm{post-SN}/a_\mathrm{pre-SN}< 1.0$|⁠) (see Fig. 7) for such binaries.

Figure 6.

The relation between the binary eccentricity e, systemic velocity |$v_\mathrm{sys}$| (km s⁻¹), and natal kick |$v_\mathrm{kick}$| (km s⁻¹) magnitudes in model V2 for ‘optimistic’ Gaia BH-like binaries. We note that a given |$v_\mathrm{kick}$| may result in a whole range of eccentricities.

Figure 7.

The relation between the change of the binary orbital separation after and just before SN (⁠|$a_\mathrm{post-SN}$|/|$a_\mathrm{pre-SN}$|⁠) and the eccentricity (e). The circles are Gaia BH1-like binaries and the crosses are Gaia BH2-like binaries. For binaries which tend to form Gaia BH1-like systems the fraction by which the orbit may change during first SN is strictly limited.

There is no similar dependence for Gaia BH2-like binaries. They harbour an already evolved star (RG) and the tidal forces between the binary components influence the orbital separation and eccentricity during Gaia BH2-like phase. In the simulations, we catch most of Gaia BH2-like binaries when they reach the upper limits for a and e that we set for Gaia BH2-like classification during their post-SN evolution. The relation between |$a_\mathrm{post-SN}$|⁠, |$a_\mathrm{pre-SN}$|⁠, and e for those systems is not straightforward.

Whether the binary orbit shrinks or expands due to the natal kick depends on the position of the exploding star on the orbit around its companion at the moment of the SN explosion and on the natal kick vector, where both the position and the natal kick vector are drawn from the distributions described in Section 2.1.

Apart from the shrinkage/widening of the orbit the natal kick changes the angle between the direction of BH spin and the direction of the binary angular momentum (⁠|$\theta$|⁠). We show the relation between angles |$\theta$|⁠, |$\gamma$|⁠, and |$v_\mathrm{sys}$| (top) and |$a_\mathrm{pre-SN}$|/|$a_\mathrm{post-SN}$| (bottom) for Gaia BH1-like and Gaia BH2-like binaries in Fig. 8. On both plots Gaia BH1-like binaries are marked with stars and Gaia BH2-like binaries are marked with filled circles. The formation of Gaia BH1-like binaries depends differently on |$\gamma$| than the formation of Gaia BH2-like binaries.

Figure 8.

The relation between the misalignment angle between the binary angular momentum vector and BH spin vector after SN (⁠|$\theta$|⁠), the angle between the natal kick velocity vector and orbital angular momentum of the system |$\gamma$|⁠, and natal kick velocity magnitude (top) or the the ratio of pre-SN to post-SN orbital separation (bottom). Both angles are in degrees units and |$v_\mathrm{nk}$| is in km s⁻¹. Gaia BH1-like systems are marked with stars and Gaia BH2-like systems with filled circles.

Most of Gaia BH1-like binaries form when the natal kick is directed out of the binary orbital plane. If |$\gamma \sim 90^{\circ }$| only a few (⁠|$< 20$|⁠) Gaia BH1-like binaries form and only with BH spin almost aligned with binary angular momentum vector (⁠|$\theta < 2^{\circ }$|⁠). The higher is the natal kick velocity the further from the binary orbital plane the kick has to be oriented and the more BH spin is misaligned from binary angular momentum in Gaia BH1-like binary, specifically |$\sim 94~{{\ \rm per\ cent}}$| of Gaia BH1-like population forms with |$\theta < 40^{\circ }$| and the median is |$\theta _\mathrm{med}\sim 10.5^{\circ }$| which can be seen in Fig. 9.

Figure 9.

The distributions of |$theta$| (top) and |$\gamma$| (bottom) for Gaia BH1-like systems. In dark colour is marked the part of the distribution corresponding to Gaia BH1-like binaries with |$v_\mathrm{sys}$| in observational ranges.

In the case of Gaia BH2-like binaries |$\gamma$| is a major factor determining their population. Most of them (⁠|$\sim 95~{{\ \rm per\ cent}}$|⁠) form if natal kick direction is close to the orbital plane (⁠|$\gamma \sim 90^{\circ }\pm 15$|⁠) and |$v_\mathrm{kick}< 41.6$| km s⁻¹. The binaries that survive SN and form Gaia BH2-like systems have |$\theta$| in the whole range of values. It can be clearly seen from the |$\gamma$| and |$\theta$| distributions shown in Fig. 10: |$\theta$| distribution is almost flat and |$\gamma$| distribution is centred around |$90^{\circ }$|⁠.

Figure 10.

The distributions of |$theta$| (top) and |$\gamma$| (bottom) for Gaia BH2-like systems. In dark colour is marked the part of the distribution corresponding to Gaia BH2-like binaries with |$v_\mathrm{sys}$| in observational ranges.

There is a dependence between the natal kick direction and the tightening/expanding of binary orbit after SN. All systems marked in black in lower plot of Fig. 8 have expanded due to the natal kick and the majority of them are Gaia BH1-like binaries. The only Gaia BH1-like binaries that have formed through tightening of the orbit after SN are those that form a pronounced narrow blue strip in bottom plot of Fig. 8. There is almost one-to-one relation between |$\gamma$| and |$\theta$| which allows the binary to shrink after SN and to form Gaia BH1-like system.

The natal kick has to tighten the orbit if the binary is to form Gaia BH2-like system. The wider is the pre-SN binary the lower velocity has to have the kick and the closer to the orbital plane it has to be oriented to become Gaia BH2-like.

The distribution of angles |$\theta$| and |$\gamma$| for Gaia BH1-like and Gaia BH2-like binaries are presented on top and bottom plots in Figs 9 and 10, respectively. The subsets of Gaia BH1-like and Gaia BH2-like binaries that have |$v_\mathrm{sys}$| matching the observations are marked in dark colours. There is a specific range in |$\gamma$| distribution which results in Gaia BH1-like systems having observed |$v_\mathrm{sys}$|⁠: most of them form if the natal kick acts at the angle |$30^{\circ }-45^{\circ }$| to the orbital plane. In general the systems with observational |$v_\mathrm{sys}$| fit into the distributions of |$\theta$| and |$\gamma$| for Gaia BH1-like and Gaia BH2-like populations.

2.3.3.2 The evolution Gaia BH-like binaries: the tidal interactions

To calculate the evolution of the orbital separation and eccentricity due to tidal interaction between the binary components, the convective damping in equilibrium tide and radiative damping in the dynamical tide (Zahn 1977) are taken into account following the formalism of Hut (1981). The discrepancies between the tidal circularization periods derived from observations of coeval binary stars populations (Meibom & Mathieu 2005) and from synthetic binary populations have led to the introduction of an additional scaling factor, |$F_\mathrm{tid}$|⁠, in the tidal evolution equations of orbital separation and eccentricity in StarTrack. It appeared that to match the observed cut off period in the open cluster M67 with simulations, the scaling factor has to be set to |$F_\mathrm{tid}=50$| in the case of convective damping (Belczynski et al. 2008). Although different approaches has been proposed to determine the circularization periods of binary populations (e.g. Zanazzi 2022; Bashi, Mazeh & Faigler 2023) and there are also new approaches to the theory of tides (e.g. Terquem 2021; Terquem & Martin 2021), it seems that the need of enhancement of convective damping is maintained (Mirouh et al. 2023).

To evaluate the impact of |$F_\mathrm{tid}$| on our results on Gaia BH-like binary formation, we tested three values of the scaling factor: |$F_\mathrm{tid}=1.0$|⁠, 50.0 and 100.0. We find that setting |$F_\mathrm{tid}=1.0$| (the standard assumption in most binary population synthesis calculations) reduces the number of Gaia BH-like binaries in the ‘V2 standard’ model by a factor of 3.2 (from 5550 to 1757 binaries), while |$F_\mathrm{tid}=100$| increases that number only by 5 per cent (from 5550 to 5815). Changing the efficiency of tidal interaction in the convective case impacts mainly the binaries that we classify as Gaia BH1-like: their number drops from 219 for |$F_\mathrm{tid}=100.0$| to 51 for |$F_\mathrm{tid}=1.0$|⁠. However, it has no such effect on Gaia BH2-like binaries (166 binaries for |$F_\mathrm{tid}=1.0$|⁠, 73 for |$F_\mathrm{tid}=50.0$| and 96 for |$F_\mathrm{tid}=100.0$|⁠). The reason for the reduction in the number of Gaia BH1-like binaries can be seen from their evolutionary path. At first, the orbit of the system widens due to the wind mass-loss from the massive primary star and, after a few Myr, the star evolves into the red supergiant phase when the radius becomes larger than |$1000\, \mathrm{R_{\odot }}$| and the star’s envelope becomes convective. This is when the tidal interaction between the stars sets in and the efficiency of the convective damping impacts the subsequent binary evolution. For |$F_\mathrm{tid}=50.0$|⁠, the efficient convective damping causes the orbital separation to start shrinking and, at the point when the evolved star’s radius reaches almost |$2000\, \mathrm{R_{\odot }}$|⁠, the binary components are close enough to enter a ce phase. On the other hand, for |$F_\mathrm{tid}=1.0$|⁠, the orbit does not shrink at all. The binary widens throughout the evolution of the massive primary due to the latter’s stellar wind, until supernova explosion disrupts the binary. That way, the initial parameter space for the binaries to form Gaia BH1-like system is much more limited for |$F_\mathrm{tid}=1.0$|⁠.

There is no such effect in the case of Gaia BH2-like binaries because they are much wider than Gaia BH1-like binaries and most of them do not need to get close enough to evolve through a ce phase. The formation of Gaia BH2-like binaries is mainly determined by the magnitude and the direction of the natal kick. The difference between the |$F_\mathrm{tid}=1.0$| and the |$F_\mathrm{tid}=100.0$| models is that most of the progenitors of Gaia BH2-like binaries start their evolution with smaller orbital separations (⁠|$a_0=4000-7000\, \mathrm{R_{\odot }}$|⁠) in the former case than in the latter case (⁠|$a_0=6000-9000\, \mathrm{R_{\odot }}$|⁠).

2.3.3.3 The evolution Gaia BH-like binaries: the common envelope parameters

The results from the IBE evolution channel depend on the choice of two parameters that determine the ce evolution of the binary, namely, the envelope binding energy parameter |$\lambda _{\mathrm{ce}}$| and the envelope efficiency parameter |$\alpha _{\mathrm{ce}}$|⁠. In this section, we discuss the choice of |$\lambda _{\mathrm{ce}}$| and |$\alpha _{\mathrm{ce}}$| for our models.

As was mentioned in Section 2.1, for the calculations done for the purpose of this paper we adopt the |$\lambda _{\mathrm{ce}}$| formulae by Wang et al. (2016) for stars more massive than |$18.0\, \mathrm{{\rm M}_{\odot }}$|⁠, which we find to be more complete compared to the older prescriptions.The reason is that to calculate the binding energy of an envelope, Wang et al. (2016) take into account not only the contributions from the gravitational binding energy, internal energy of the envelope consisting of the thermal energy, radiation energy, and ionization and dissociation energies but also the enthalpy (see Ivanova & Chaichenets 2011, and discussion therein). The inclusion of enthalpy decreases the envelope’s binding energy with respect to the cases where it is not taken into account. This approach differs from the commonly used |$\lambda _{\mathrm{ce}}$| prescription given by Claeys et al. (2014), where only gravitation and ionization are considered as the energy sources that unbind the donor’s envelope. To see how different choices of |$\lambda _{\mathrm{ce}}$| impact the envelope binding energy for stars of different masses and metallicities, we refer to fig. A1 in the appendix of Iorio et al. (2023), although it does not present the case of the Wang et al. (2016) formulae. From that figure, it can be seen that, for a |$M_\mathrm{ZAMS}\sim 30\, \mathrm{{\rm M}_{\odot }}$| star of metallicity |$Z=0.01$|⁠, the envelope binding energy during the core He-burning phase is |$\sim 2$| orders of magnitude smaller for the |$\lambda _{\mathrm{ce}}$| of Claeys et al. (2014) than for the |$\lambda _{\mathrm{ce}}$| of Xu & Li (2010b). Consequently, the envelope ejection during a ce phase is easier in the former than in the latter case, for the same value of |$\alpha _{\mathrm{ce}}$| in both cases. This, in turn, may result in different post-CE orbital separations in models where, apart from |$\lambda _{\mathrm{ce}}$|⁠, all other assumptions are identical.

The standard version of StarTrack has been using |$\lambda _{\mathrm{ce}}$| formulae from Xu & Li (2010a, b). However, for the stars with ZAMS masses in the range of |$\sim 35-45\, \mathrm{{\rm M}_{\odot }}$|⁠, which are the progenitors of |$\sim 8-10\, \mathrm{{\rm M}_{\odot }}$| BHs, those formulae give |$\lambda _{\mathrm{ce}}$| values as high as 4.0 for the core He burning stars because they underestimate their envelope binding energy. The reason is that the formulae in Xu & Li (2010a, b) were fitted to the models calculated with the old Eddington’s code ev only for stars with |$M_\mathrm{ZAMS}< 20\, \, \mathrm{{\rm M}_{\odot }}$|⁠. Those formulae were updated using Modules for Experiments in Stellar Astrophysics code (mesa; Paxton et al. 2011) in 2016 by Wang et al. (2016) and should be used in their most recent form (Xiang-Dong Li private communication). We further motivate our choice of the Wang et al. (2016) |$\lambda _{\mathrm{ce}}$| by the fact that the formulae were fitted using a dense grid of stellar models calculated with mesa. Furthermore, they assumed a self-consistent core-envelope boundary located in the hydrogen-burning shell of the maximum compression prior to the ce, following Ivanova (2011; see, however, Klencki et al. 2021). The formulae from Wang et al. (2016), that we use, were calculated using the prescription for the wind mass-loss rates from Hurley et al. (2000) and Vink et al. (2001) (Wind1 in Wang et al. 2016) and we take |$\lambda _{\mathrm{ce}}$| as the average of |$\lambda _\mathrm{g}$| (corresponding to gravitational binding energy only), |$\lambda _\mathrm{b}$| (corresponding to the total energy, comprising both gravitational and internal energy) and |$\lambda _\mathrm{h}$| (which, in addition to the total energy, takes into account the enthalpy of the stellar envelope).

Note that our results are obtained with the standard value of |$\alpha _\mathrm{ce}$|⁠, i.e. |$\alpha _\mathrm{ce}=1.0$|⁠. There are analytical estimates for the final orbital separation, after the ce phase in pre-Gaia BH1 systems (e.g. El-Badry et al. 2023a). However, they are based on the assumption that the donor star has to lose a mass of |$20\, \mathrm{{\rm M}_{\odot }}$| in envelope, during the ce. This is an overestimation, considering that the BH progenitor loses up to half of its mass in stellar winds during its evolution, before the ce phase starts (see Fig. 1 and its description in Section 2.2.2). Consequently, the envelope mass, that has to be unbound during the ce phase, has a mass of only |$\sim 3\, \mathrm{{\rm M}_{\odot }}$|⁠. If |$M_\mathrm{env}=3.0\, \mathrm{{\rm M}_{\odot }}$| is considered in the calculations instead, the post-CE orbital separation of the binary turns out to be |$> 100\, \mathrm{R_{\odot }}$| instead of a few |$\mathrm{R_{\odot }}$|⁠.

Taking into account the uncertainties of the |$\lambda _{\mathrm{ce}}$| prescription and the correction to the post-CE orbital separation calculations as mentioned above, we conclude that invoking |$\alpha _{\mathrm{ce}}> 1$| should not be necessary to form Gaia BH-like binaries via the IBE channel. Also, it is so far unclear what (if any) physical mechanism could lead to |$\alpha _{\mathrm{ce}}> 1$|⁠. Therefore, we assume |$\alpha _{\mathrm{ce}}=1$|⁠, which still seems to be an optimistic assumption as it implies that all orbital energy goes into unbinding the envelope without any energy loss.

For the analyses of the impact of different |$\alpha _{\mathrm{ce}}$| values, of different SN engines (e.g. rapid and delayed), and of the reduction of SN natal kicks due to matter fallback on the isolated-binary formation of BH–star binaries, we refer to the works of Yamaguchi et al. (2018) and Chawla et al. (2022).

2.3.4 Gaia BH1 as low-mass X-ray binary

There is no X-ray emission detected from Gaia BH1 and Gaia BH2. However, it is found that during the post-dormant BH binary phase, both of these binaries may evolve through the symbiotic X-ray binary phase to the phase of stable Roche-lobe overflow (RLOF) Rodriguez et al. (2023). At this point, the systems convert into LMXBs. The population of LMXBs, that could arise from Gaia BH-like systems, would have long orbital periods of |$P_\mathrm{orb}> 100$| d.

Of the two binaries that we consider as the best examples of Gaia BH1 and Gaia BH2 in our IBE simulations and which evolution we describe in Section 2.2.2 only Gaia BH1-like system undergoes RLOF becoming LMXB in its subsequent evolution while Gaia BH2-like binary remains detached.

As of today, all known LMXBs with BH accretors are transient sources, that is they show periodic rapid increase in their luminosity up to |$L_\mathrm{X}\sim 10^{39}$| erg s⁻¹ lasting from days to years, followed by long periods (tens of years) of quiescence when their luminosity drops below |$L_\mathrm{X}\sim 10^{33}$| erg s⁻¹. The outbursts are triggered by an instability in accretion disc formed around BH. The instability arises from the changes in opacities when ionization/recombination sets in, what in turn entails the changes in the disc viscosity (for details of disc instability model see e.g. Lasota 2001; Hameury 2020). According to the model the system is persistently bright (has high level of luminosity with no outbursts) if the mass transfer rate |$\dot{M}_\mathrm{tr}$| is higher than critical value of mass accretion rate at the outer disc radius (⁠|$\dot{M}^{+}_\mathrm{crit}(R_\mathrm{d})$|⁠) and the system is transient (undergoes outbursts) if |$\dot{M}_\mathrm{tr}< \dot{M}^{+}_\mathrm{crit}(R_\mathrm{d})$| and |$\dot{M}_\mathrm{tr}> \dot{M}^{-}_\mathrm{crit}(R_\mathrm{in})$| where |$\dot{M}^{-}_\mathrm{crit}(R_\mathrm{in})$| is the critical mass accretion rate at the inner disc radius. The formulae for |$\dot{M}_\mathrm{crit}(R)$| are taken from Lasota, Dubus & Kruk (2008):

$$\begin{eqnarray} \begin{split} & \dot{M}^{+}_\mathrm{crit}(R)=8.07\times 10^{15}\alpha _{0.1}^{-0.01}\, R_{10}^{2.64}\, M_1^{-0.89}\, & \mathrm{g\, s^{-1}} \\ & \dot{M}^{-}_\mathrm{crit}(R)=2.64\times 10^{15}\alpha _{0.1}^{0.01}\, R_{10}^{2.58}\, M_1^{-0.85}\, &\mathrm{g\, s^{-1}} \end{split} \end{eqnarray}$$

(3)

where |$\alpha _{0.1}=0.1\alpha$| is the viscosity parameter in units 0.1, |$R_{10}$| is the radius in units of |$10^{10}$| cm, and |$M_1$| is the primary mass in solar masses.

In following calculations, we assumed that the disc outer radius extends to 90 per cent of BH Roche lobe radius in the outburst and that the inner disc radius |$R_\mathrm{in}$| does not coincide with the last stable orbit in the quiescence but instead it is truncated due to the transition from optically thick to optically thin radiatively inefficient flow (e.g. Lasota, Narayan & Yi 1996; De Marco et al. 2021). We further assumed that the truncation happens at |$3.4\times 10^4 R_\mathrm{g}$| where |$R_\mathrm{g}$| is gravitational radius |$R_\mathrm{g}=GM_\mathrm{BH}/c^2$|⁠. |$\dot{M}^{-}_\mathrm{crit}(R)$| decreases with R and the assumption about |$R_\mathrm{in}$| is the upper limit for the inner disc truncation that is justified by the observations of V404 Cyg (a transient LMXB; Bernardini et al. 2016). The viscosity parameter in the hot disc is taken to be |$\alpha _{0.1}=1.0$| and in the cold disc |$\alpha _{0.1}=0.2$| (Lasota 2001). We calculated the critical mass accretion rates at the outer and inner disc radii for parameters that we found for exemplary Gaia BH1-like systems described in Section 2.2.2 at the onset of LMXB phase:

$$\begin{eqnarray} \begin{split} & \dot{M}^{+}_\mathrm{crit}(R_\mathrm{d})=7.95\times 10^{-6}\, & \mathrm{M_\odot \, yr^{-1}} \\ & \dot{M}^{-}_\mathrm{crit}(R_\mathrm{in})=3.93\times 10^{-10}\, & \mathrm{M_\odot \, yr^{-1}} \end{split} \end{eqnarray}$$

(4)

for |$R_\mathrm{d}=20.44$||$\mathrm{R_{\odot }}$| and |$R_\mathrm{in}=0.71$||$\mathrm{R_{\odot }}$|⁠.

The mass transfer rate for the considered binary in its LMXB state is |$\sim 3.53\times 10^{-8}\, \mathrm{{\rm M}_{\odot }}\, \mathrm{yr}^{-1}$|⁠, therefore, according to the disc instability model, the system should become a transient LMXB. If LMXBs that formed from Gaia BH1-like binaries exist, they would be observed through their outbursts. It is of interest to predict what would be the outburst recurrence time of such systems. The majority of the disc mass is accreted on to the compact object during the long outburst. The disc depleted from matter needs to refill so that another outburst may develop. Following Hameury & Lasota (2020), the refilling time corresponding to the outburst recurrence time is:

$$\begin{eqnarray} t_\mathrm{recc}=\frac{\Delta M}{\dot{M}_\mathrm{tr}}, \end{eqnarray}$$

(5)

where |$\Delta M$| is the amount of the disc mass accreted during outburst and |$\dot{M}_\mathrm{tr}$| is the mass transfer rate.

We assumed that the amount of mass retained in the disc at the end of outburst is 10 per cent of the maximum disc mass (just before the outburst is triggered) and for the purpose of this calculations we took |$\Delta M=0.9M_\mathrm{d}$|⁠, where the disc mass is:

$$\begin{eqnarray} M_\mathrm{d}=3.3\times 10^{26}\psi \alpha _{0.2}^{-0.8}\dot{M}_{19}^{0.7}M_1^{0.25}R_{12}^{1.25}\, \, \mathrm{g}, \end{eqnarray}$$

(6)

where |$\psi$| is a factor connected to opacity (for details see Hameury & Lasota 2020) and we used |$\psi =1.3$| following the authors, |$\alpha$| is the disc viscosity parameter in units 0.2 (Lasota 2001), |$\dot{M}$| is the accretion rate in units of |$10^{19}\, \mathrm{g\,s^{-1}}$|⁠, |$M_1$| is the compact object mass in solar units, and |$R_{12}$| is the outer disc radius in outburst in units of |$10^{12}$| cm.

The mass transfer rate (here equal to mass accretion rate)² changes in time with the nuclear evolution of the donor and so does the Roche lobe radius of the accretor. Therefore, we estimate two limiting values of |$t_\mathrm{recc}$|⁠: at the onset (minimum |$t_\mathrm{recc}$|⁠) and at the end (maximum |$t_\mathrm{recc}$|⁠) of LMXB phase of the considered Gaia BH1-like system. The parameters used for calculations are given in Table A1 in the Appendix. The obtained recurrence times are |$\sim 47.3$| and |$\sim 159.8$| yr, respectively. This suggests that there may be a population of long orbital period transient LMXBs that originated from dormant BH binaries and that are hard to observe because their luminosity rises above the detection threshold only once every few tens of years.

3 GAIA BH BINARIES FROM YOUNG MASSIVE AND OPEN STAR CLUSTERS

3.1 Computed model clusters

It would be worth considering and comparing our isolated-binary formation rate of Gaia BH1- and Gaia BH2-like binaries with those from star clusters. Dynamical interactions among stars, stellar binaries, and BHs in star clusters, particularly in low-mass open clusters, is widely proposed to be a channel for assembling detached BH-star binaries (Banerjee 2018b; Niccolò Di Carlo et al. 2023; Rastello et al. 2023; El-Badry et al. 2023a; Tanikawa et al. 2023a). In this work, we consider the evolutionary model star cluster set of Banerjee (2021a, hereafter Ba21) to estimate the formation rate of dynamical BH-star binaries in the galactic field.

We waive a detailed description and discussion of the computed models in Ba21, which are available in the Ba21 paper itself and are discussed further in several follow-up studies (e.g. Banerjee 2020, 2021b; Belczynski et al. 2022). To summarize, the Ba21 model set comprise 72 model star clusters that evolve from a young massive phase to an old, open cluster-like phase. The clusters’ initial mass cover the range of |$10^4\mathrm{{\rm M}_{\odot }}\le M_{\rm cl}\le 10^5\mathrm{{\rm M}_{\odot }}$| and initial size (half-mass radius) range between |$1-3$| pc. About half of the clusters have an initial overall primordial binary fraction between |$5~{{\ \rm per\ cent}}\, \mathrm{ and}\, 10~{{\ \rm per\ cent}}$|⁠; however, the clusters’ massive-star members, above |$m_\ast > 16\mathrm{{\rm M}_{\odot }}$|⁠, are paired separately among themselves with an initial binary fraction of 100 per cent. The clusters cover the metallicity range |$0.0001 \le Z \le 0.02$|⁠. All models are subjected to a solar-neighbourhood-like galactic tidal field.

The cluster models are computed with the direct, star-by-star N-body integration code nbody7 (Aarseth 2012), with various updates as described in Ba21. In particular, stellar and binary evolution and remnant formation take place in the models via the coupling of nbody7 with an updated version of the fast binary evolution code bse (Hurley, Tout & Pols 2002; Banerjee et al. 2020). In most models, stellar remnants (NSs and BHs) form via the ‘rapid’ remnant mass model of Fryer et al. (2012), as implemented in Banerjee et al. (2020).

As seen in the above mentioned followup works of Ba21, these cluster models produce gravitational waves (GW) mergers with properties and rate consistent with those inferred from GW observations. The main motivation for choosing the Ba21 models for this study is that they are evolved either at least until all the in-cluster BHs are practically depleted or until |$\approx 11$| Gyr (see appendix C of Ba21). Therefore, the population of escaped BH-star binaries (and other escaped BH-containing systems) from these clusters can be considered to be the complete yield, allowing for robust rate estimates for such field objects from clusters. This is unlike any other recent study addressing the dynamical formation of detached BH-star binaries (see above). Therefore, in the following, we shall focus on only those BH-star binaries that have genuinely escaped from the clusters, i.e. those which have crossed their parent cluster’s tidal radius and have become tidal-tail members. This would allow for a fair comparison with the yield from isolated binaries (see above) and with Gaia BH1 and Gaia BH2, both of which are field objects.

3.2 Results: formation rate of field BH-star binaries

Fig. 11 shows the orbital period, |$P_{\rm orb}$|⁠, eccentricity, e, and companion mass, |$M_{\rm comp}$|⁠, against BH mass, |$M_{\rm BH}$|⁠, of all the BH-star binaries that have escaped into the galactic field from the model clusters of Ba21. In all panels, the individual data points are colour coded according to the velocity of ejection, |$v_{\rm ej}$|⁠, which is the velocity of the escaping binary’s centre of mass, at ejection. Unless otherwise stated, all quoted or plotted values of a BH-star binary are those when the escaping binary crosses the instantaneous tidal radius of its parent cluster; i.e. at the moment of its ‘formal’ escape form the cluster. These confirmed escaped binaries and their parameters are recovered from the nbody7 simulation logs, with the help of the ‘BINARY ESCAPE’ records. In Fig. 11, the four rows distinguish between assembly mechanism (dynamically or primordially assembled) and the stellar companion’s evolutionary status (main sequence or evolved) of the binary.

Figure 11.

Demographics of the BH-star binaries that have escaped into the galactic field from the evolutionary model star clusters of Ba21. The filled circles in the panels show the mass of the BH member (⁠|$M_{\rm BH}$|⁠; X-axes) against the orbital period (⁠|$P_{\rm orb}$|⁠; left column), eccentricity (middle column), and companion-star mass (⁠|$M_{\rm comp}$|⁠; right column) of these model BH-star binaries. The binaries are plotted separately, based on the evolutionary stage of the companion star (MS or evolved) and the assembly channel (primordial or dynamical), as indicated in each row’s title. The data points are colour-coded according to the velocity of ejection of the binary (⁠|$v_{\rm ej}$|⁠; colour bar). All of the plotted values correspond to the event of the BH-star binary crossing the instantaneous tidal radius of its parent cluster. For comparison, the observed BH-star binary candidates Gaia BH1, that contains a main-sequence companion, and Gaia BH2, that contains an evolved companion, are shown in the upper and lower pair of rows, respectively (empty squares). For these observed data points, the colour coding is not followed and their error bars are generally invisible due to the scale of the figure axes. (The size of the squares are chosen for legibility and does not represent measurement uncertainties.)

The figure demonstrates that escaped BH-star binaries with a low-mass companion (⁠|$M_{\rm comp}\lesssim 3\mathrm{{\rm M}_{\odot }}$|⁠) are most likely dynamically assembled, as anticipated by El-Badry et al. (2023a). In the Ba21 clusters, such pairing between a low-mass star and a BH can happen only via dynamical channel since, initially, massive, BH-progenitor stars are always primordially paired among themselves (and never with a low-mass star). The dynamically assembled BH–MS binaries generally cover the |$P_{\rm orb}$|⁠, eccentricity, |$M_{\rm comp}$|⁠, and peculiar velocity range of Gaia BH1. However, the escaped BH-evolved binaries (both dynamically and primordially paired) struggle to meet the observed properties of Gaia BH2. As such, there are only a handful of BH-evolved binaries that have escaped from the entire Ba21 set.

If all orbital parameters are compared together, there is only one ejected BH-star binary that is somewhat comparable to Gaia BH1: this binary has |$M_{\rm comp}=1.1\mathrm{{\rm M}_{\odot }}$| (main sequence), |$M_{\rm BH}=17.2\mathrm{{\rm M}_{\odot }}$|⁠, |$P_{\rm orb}=290$| d, |$e=0.31$|⁠. The total computed initial star cluster mass of the Ba21 set being |$M_{\rm cl,tot}=3.6\times 10^6\,\mathrm{{\rm M}_{\odot }}$|⁠, this gives an unweighted (plain) formation rate of Gaia BH1-like system, per unit initial cluster mass, of |$2.7\times 10^{-7}\,\mathrm{{\rm M}_{\odot }}^{-1}$|⁠. This particular binary is formed via a dynamical exchange encounter in one of the |$M_{\rm cl}=2\times 10^4\,\mathrm{{\rm M}_{\odot }}$| cluster and is also ejected dynamically from the cluster. In the Ba21 set, there is no ejected BH-star system that is or can evolve (after the ejection) anywhere similar to Gaia BH2.

Table 5 (middle column) shows the plain formation rates for escaped BH-star systems with a low-mass companion (⁠|$M_{\rm comp}\le 3\mathrm{{\rm M}_{\odot }}$|⁠), with a solar mass companion (⁠|$0.8\mathrm{{\rm M}_{\odot }}\le M_{\rm comp}\le 1.1\mathrm{{\rm M}_{\odot }}$|⁠), and of Gaia BH1-type, as obtained from the Ba21 set. The table (right column) also states the corresponding weighted rates when the initial masses of the clusters are assumed to follow a power-law distribution of index |$-2$|⁠, as motivated by observations (e.g. Portegies Zwart, McMillan & Gieles 2010; Bastian et al. 2012). The weighted rates are similar to the plain rates since there are only a handful of field low-mass-companion systems (21, 5, and 1 for the three cases, respectively) from a number of model clusters spanning over a wide initial mass range (see Section 3.1). Fig. 12 shows the mass, orbital parameter, and |$v_{\rm ej}$| distributions of the escaped BH-star binary population in Ba21, with low mass (orange line and circles) and solar mass (blue line and circles) companion. The distributions are constructed assuming a |$\propto M_{\rm cl}^{-2}$| initial cluster mass distribution. Note that there is only one specimen that is comparable to a Gaia-detected BH-star binary (see above). The distributions are nevertheless presented keeping in mind the overall diverse demographics of field BH-star binary population from these computed star clusters (c.f. Fig. 11), that may address such binaries discovered in the future.

Figure 12.

Normalized distributions of BH mass (⁠|$M_{\rm BH}$|⁠), orbital period (⁠|$P_{\rm orb}$|⁠), eccentricity, companion star mass (⁠|$M_{\rm comp}$|⁠), and ejection velocity (⁠|$v_{\rm ej}$|⁠) of BH-star binaries that have escaped into the galactic field from the evolutionary model star clusters of Ba21 (upper to lower panels; all values correspond to the event of the BH-star binary crossing the instantaneous tidal radius of its parent cluster). The distributions of binaries with a low-mass companion (⁠|$M_{\rm comp}\lt =3\,\mathrm{{\rm M}_{\odot }}$|⁠; orange (light) line) and with a solar mass companion (⁠|$0.8\,\mathrm{{\rm M}_{\odot }}\lt =M_{\rm comp}\lt =1.1\mathrm{{\rm M}_{\odot }}$|⁠; blue (dark) line) are shown separately. These distributions incorporate a weighted distribution of cluster birth mass, |$M_{\rm cl}$|⁠, according to the power law |$\propto M_{\rm cl}^{-2}$|⁠. Furthermore, the two distributions on each panel are normalized according to the relative number of BH-low mass and BH-solar mass binaries ejected from the weighted cluster population. The observed Gaia BH1 and Gaia BH2 binaries are indicated in the panels; in the final panel, peculiar velocities of the Gaia BH binaries are indicated.

Table 5.

Formation rates of field BH-star binaries from young massive and open star clusters. The rates in the first three rows are based on the evolutionary model star cluster set of Ba21. The Gaia BH2-like systems’ rate in the final row is based on a newly computed cluster model grid (see text).

Case	Rate (plain) |$[\mathrm{{\rm M}_{\odot }}^{-1}]$|	Rate (weighted) |$[\mathrm{{\rm M}_{\odot }}^{-1}]$|
Low mass companion	|$5.8\times 10^{-6}$|	|$5.3\times 10^{-6}$|
Solar mass companion	|$1.4\times 10^{-6}$|	|$1.4\times 10^{-6}$|
Gaia BH1-like	|$2.7\times 10^{-7}$|	|$3.5\times 10^{-7}$|
Gaia BH2-like	|$7.4\times 10^{-7}$|	–

As evident from Figs 11 and 12, the vast majority of dynamically assembled field BH-star systems have |$M_{\rm comp}< 3\mathrm{{\rm M}_{\odot }}$|⁠. Nearly half of them have zero or small eccentricity, an MS companion, and |$P_{\rm orb}< 1$| d – these very tight binaries are outcomes of tidal capture between the BH and the companion star inside the parent cluster and the binary’s subsequent tidal circularization.³ Inside clusters, such close binaries behave like single, massive objects most of the time and are ejected dynamically in ways similar to those of ejection of single BHs (Kulkarni, Hut & McMillan 1993). These systems are potentially tidally interacting or mass transferring and are not compatible with the Gaia detected, detached BH-star binaries. There are also a few wide, |$P_{\rm orb}\gtrsim 100$| day, MS-companion binaries – these systems are assembled in the cluster via dynamical exchange encounters and are also ejected via dynamical scatterings (see above), making them eccentric in general.

Inside a cluster, dynamical formation of wide and eccentric BH-star binaries with a low-mass stellar companion is rather common (Banerjee 2018b). After formation, such a binary would undergo dynamical hardening whilst inside the cluster. However, the BH-star pair is also vulnerable inside a cluster containing a dynamically active BH subsystem, since the low-mass stellar companion can easily be displaced by a much more massive BH via a close exchange encounter. In other words, the in-cluster BH-star phase is at best a short-lived, intermediate phase (Banerjee 2018a). This explains the general dearth of field BH-low mass star systems from the computed models.

Dynamical ejection of wide BH–RG systems is further suppressed: owing to the large envelope size of the RG, an in-cluster close interaction with such a binary almost inevitably results in tidal interaction and/or geometrical interception with the giant’s envelope, in turn resulting in a merger, ce, or stable mass transfer with the RG. Indeed, for the Ba21 clusters, all dynamically paired field BH–RG systems are stable-mass-transferring. Inspection reveals that the interactions, that ejected these binaries, are also responsible for putting them in a mass transferring state. However, it is possible that an escaped, wide enough (⁠|$P_{\rm orb}\sim 10^3$| d), eccentric BH–MS binary evolves into a detached, Gaia BH2-like BH–RG binary.

While no such case is recovered from the Ba21 set, one such field BH–MS binary has been dynamically produced by a newly computed model cluster set. This set comprises a grid of model clusters with homogeneous initial conditions and properties: mass |$M_{\rm cl}=3\times 10^4\,\mathrm{{\rm M}_{\odot }}$|⁠, size 2 pc, ‘delayed’ remnant mass model (Fryer et al. 2012), metallicities |$Z=0.0002$|⁠, 0.001, 0.005, 0.01, 0.02, and overall and massive-star primordial binary fractions (see above) of 10 per cent and 100 per cent, respectively. The distributions of the primordial binaries are similar to those of Ba21. For each Z, model clusters are evolved by placing them on circular orbits on the equatorial plane of an axisymmetric, static, Milky Way-like, central mass-disc-halo galactic potential (Allen & Santillan 1991), at different galactocentric distances, namely, |$r_{\rm g}=1.0$|⁠, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 6.0, and 8.5 kpc (hence, the grid contains 45 evolutionary models in total). All of these clusters are evolved (using nbody7) until they are completely dissolved into the external galactic field (otherwise, at least until all BHs are depleted from them). A detailed account of this model set (hereafter nicknamed as ‘|$r_{\rm g}$|-grid’) and its outcomes will be described in a forthcoming dedicated paper.

As for the BH-star yield, |$r_{\rm g}$|-grid has produced one field BH–MS binary with |$M_{\rm BH}=14.8\mathrm{{\rm M}_{\odot }}$|⁠, |$M_{\rm comp}=1.1\mathrm{{\rm M}_{\odot }}$|⁠, |$P_{\rm orb}=1200$| d, |$e=0.56$|⁠. This binary is qualitatively similar to the binary demonstrated in Fig. 2 in the latter’s BH–MS phase, and would evolve into a Gaia BH2-type detached (non-interacting) binary. The total simulated initial mass of |$r_{\rm g}$|-grid being |$M_{\rm cl,tot}=1.35\times 10^6\mathrm{{\rm M}_{\odot }}$|⁠, the unweighted formation rate of field Gaia BH2-type systems from the set is |$7.4\times 10^{-7}\mathrm{{\rm M}_{\odot }}^{-1}$|⁠. These clusters have not produced any field BH-star system that is compatible with Gaia BH1.

Notably, |$v_{\rm ej}$|s of the field BH-star binaries with low-mass companion are typically lower than the observed peculiar velocities of the Gaia BH binaries (see Fig. 12, lowermost panel). However, the |$v_{\rm ej}$| distributions have long tails that extend well beyond the velocities of both Gaia BH1 and Gaia BH2. In the scenario of dynamical formation and ejection from moderately massive clusters such as those in Ba21, a generally low |$v_{\rm ej}$| is expected due to the clusters’ low-velocity dispersion and escape velocity. Still, high-velocity ejections can rarely occur due to chaotic interactions inside the clusters, invoving compact subsystems. The cluster Gaia BH1 and Gaia BH2 rates quoted in Table 5 do not include any velocity criterion.

Since the model clusters produce only one Gaia BH1- and one Gaia BH2-like binary (see above), we do not repeat the source count analysis in Section 2.3.2 for the cluster channel. We note that the relative contribution of star clusters to the field Gaia-BH binaries depends on cluster formation efficiency (hereafter CFE), i.e. the relative fractions of the star forming mass that go into star clusters and isolated binaries. An appropriate CFE for a Milky Way-like spiral galaxy is largely an open question (see e.g. Krumholz, McKee & Bland-Hawthorn 2019 for a review). The recent study by Niccolò Di Carlo et al. (2023) adopted a CFE of 10 per cent. Notably, the number, |$\sim 900$|⁠, of isolated-binary Gaia BH systems, as obtained in Section 2.3.2, corresponds to the highest formation efficiency from the isolated-binary channel, which formation efficiency is similar to that for the cluster channel (c.f. Tables 3 and 5). Hence, the 5-yr count of Gaia BH-like systems from clusters would simply proportionate with the CFE, with respect to the above isolated-binary-channel count. A detailed study of CFE in the appropriate cluster mass range is beyond the scope of this work.

3.3 Discussions: on spin-orbit inclinations of dynamical BH-star binaries

As for the spin-orbit alignment, it can be expected that wide, dynamically assembled field BH-star binaries (resembling Gaia BH1 and BH2; see above) would have an isotropic spin-orbit misalignment w.r.t. the BH member, since the exchanging BH is generally uncorrelated to the stellar binary (in a similar way as for dynamical binary black hole mergers; see e.g. Banerjee, Olejak & Belczynski 2023). In case the pre-exchange (in-cluster) star-star binary is a primordial binary, it may be spin-orbit aligned to some extent. However, after the exchange, the exchanged BH’s spin would still be uncorrelated to and hence isotropically misaligned w.r.t. the binary’s orbital angular momentum and the stellar companion’s spin. Inside the cluster, close dynamical encounters would impart additional, random, common spin-orbit tilts to the (BH-star) binary’s members (Trani et al. 2021), so that the isotropic distribution is maintained until the binary’s escape from the cluster.

On the other hand, for the very tight field BH-star systems that are formed via in-cluster, dynamically mediated tidal interactions (see above; which are unlike Gaia BH systems), the tidal interaction and circularization may align the stellar member’s spin w.r.t. the orbital angular momentum. However, the BH member would still remain misaligned, presumably isotropically, at least initially. In the case of a mass transfer forming an accretion disc around the BH in such a binary, the system can exhibit low-frequency quasi-periodic oscillations (QPO) related to Lense–Thirring precession of the disc (Ingram, Done & Fragile 2009, and references therein). In other words, stellar clusters may dynamically produce a population of field stellar mass BH candidates that exhibit potential misalignment signatures such as low-frequency QPOs. Although interesting, a detailed study of this is beyond the scope of this work.

A drawback of the present estimates from clusters is that the criteria for selecting Gaia BH1 and Gaia BH2 systems could not be matched with those for the isolated-binary counterpart. This is simply due to the sparsity of escaped Gaia BH systems (see above), which makes applying well defined selection criteria nearly impossible. Note that, in this work, we have simply compared the two formation channels and have not attempted to combine them, unlike in some recent works (Niccolò Di Carlo et al. 2023).

4 CONCLUSIONS

In this paper, we have investigated the two possible channels of the formation of dormant BH binaries Gaia BH1 and Gaia BH2: the isolated binary evolution (Section 2) and the dynamical interactions in young massive and open star clusters (Section 3). We used the population synthesis code StarTrack for IBE calculations (Section 2.1) and the updated star-by-star N-body integration code nbody7 for cluster models (Section 3.1). The two systems from IBE Gaia-BH like binaries population, that match well the parameters of Gaia BH1 and Gaia BH2, were evolved till the Hubble time to find the characteristics of the low-mass X-ray binaries that possibly form from such systems (Section 2.3.7). We have also analysed what are the constrains on magnitude and direction of the natal kick imparted to BH set by the orbital parameters and systemic velocities of Gaia BH1-like and Gaia BH2-like binaries (Section 2.2.3 and Section 2.3.4).

As seen in Section 3.2, medium mass (⁠|$\sim 10^4\, \mathrm{{\rm M}_{\odot }}$|⁠), pc-scale stellar clusters struggle to produce dynamically paired Gaia BH-type binaries in the field. Nevertheless, the overall population of field BH-star binaries from such clusters does present interesting and diverse demographics. In the case when no other constrains (⁠|$P_{\rm orb}$|⁠, e, |$M_\mathrm{BH}$|⁠) were imposed on the population of binaries except that they should consist of BH and a low-mass star, the formation rate of field BH-star binaries is |$1.4-5.8\times 10^{-6}\, \mathrm{{\rm M}_{\odot }}^{-1}$| depending on the mass of the companion. This implies that further types of BH-star systems can potentially be discovered in the future. Of course, internally, a cluster continues to assemble BH-star systems (as long as its BH population is not fully depleted), potentially at a much higher rate. The differences between internal and field BH-star population from clusters will be explored in a future work.

In contrast, the IBE simulations result in a large Gaia BH-like binaries population, but this is the effect of pre-choosing possible Gaia BH-like progenitors as described in Section 2.3.1. Among Gaia BH-like binaries that are produced through IBE there are systems with the parameters matching closely those inferred from observations of Gaia BH1 and Gaia BH2.

The results from IBE simulations were obtained using the standard value of parameter |$\alpha _\mathrm{ce}$| which describes the fraction of the orbital energy that is used to unbind the envelope in ce phase, that is |$\alpha _\mathrm{ce}=1.0$| (see Section 2.3.6).

While we acknowledge the uncertainties that the stellar and binary evolution modelling is subject to, our results indicate that isolated binary evolution should be considered as the scenario equivalent to dynamical interactions in open clusters for Gaia BH1 and Gaia BH2 binaries formation.

The formation rate of Gaia BH-like binaries, as inferred from the present simulations, is sensitive to the criteria adopted to classify a binary as Gaia BH-like (Section 2.2.1) and it can differ by three orders of magnitude for the IBE channel (see Table 3). Ensuring that the common criteria are set and the same total stellar mass per simulation is adopted for IBE and dynamical interactions in young clusters, we find that the formation rate of Gaia BH-like binaries for both channels is equal |$2.7\times 10^{-7}\, \mathrm{{\rm M}_{\odot }}^{-1}$|⁠. If we consider the general population of field BH-low-mass star (⁠|$< 3\, \mathrm{{\rm M}_{\odot }}$|⁠) binaries from young massive and open clusters and corresponding BH-low mass star population from IBE the formation rate of such binaries in both channels is on order of |$\sim 10^{-6}\, \mathrm{{\rm M}_{\odot }}^{-1}$|⁠.

Our approximate estimate is that |$\sim 900$| dormant BH binaries, formed through the IBE channel, would be detectable by Gaia in the Milky Way thin disc.

Considering the IBE channel, the formation of Gaia BH1- and Gaia BH2-like binaries depends not only on the natal kick magnitude but also on two angles: the angle between the binary angular momentum and the natal kick vectors (⁠|$\gamma$|⁠), and the misalignment between the binary angular momentum and the BH-spin vectors (⁠|$\theta$|⁠). We summarize the results that we inferred from the model Gaia BH1- and Gaia BH2-like population as follows:

• All Gaia BH1- and Gaia BH2-like binaries have systemic velocities |$< 80$| km s⁻¹ and the peculiar velocities at birth (which correspond to |$v_\mathrm{sys}$| in the model) of the observed systems fit well into the |$v_\mathrm{sys}$| distributions of Gaia BH1- and Gaia BH2-like binaries.

• The natal kick velocity has to be lower than 100 km s⁻¹ for Gaia BH1- and Gaia BH2-like binaries to form.

• The population of Gaia BH1-like binaries is not vulnerable to the natal kick velocity distribution ^adopted but a seven fold increase in the number of Gaia BH2-like binaries in model V2 compared to model V1 favours the two-component Maxwellian distribution with low- and high-velocity components (⁠|$\sigma \approx 21$| and |$\sigma \approx 107$| km s⁻¹), as suggested by Zhao et al. (2023).

• The median value of natal kick velocity imparted to the systems, that become Gaia BH1-like and Gaia BH2-like binaries, is |$\sim 39$| km s⁻¹ for the former and |$\sim 19$| km s⁻¹ for the latter (assuming the natal kick distribution from Zhao et al. 2023).

• |$\sim 94~{{\ \rm per\ cent}}$| of Gaia BH1-like binaries that form in IBE channel have |$\theta < 40^{\circ }$| and the median value is |$\theta _\mathrm{med}\sim 10^{\circ }$|⁠.

• For Gaia BH1-like binary to be formed, the kicks should be directed out of the binary orbital plane: the higher is the natal kick velocity the further from orbital plane should the kick be directed.

• |$\sim 95~{{\ \rm per\ cent}}$| of the binaries that become Gaia BH2-like experience the natal kick directed within |$\pm 15^{\circ }$| from the orbital plane.

We find that, in addition to the natal kicks, the treatment of the tidal interactions in binaries should not be overlooked in the context of the formation of Gaia BH-like binaries. The enhancement of the convective damping efficiency relative to the classical prescription Hut (1981), introduced as the additional factor |$F_\mathrm{tid}$| in the equations, increases the number of the systems like Gaia-BH1 in the IBE channel. However, increasing |$F_\mathrm{tid}$| above the certain value, which was originally calibrated to match the observations (⁠|$F_\mathrm{tid}=50$|⁠, Belczynski et al. 2008), does not lead to further significant increase in the Gaia BH1-like binaries number (see Section 2.3.5). The results show that the prescriptions for tidal interactions that are used in the codes play a decisive role in the formation of Gaia BH1-like binaries via the IBE channel and that more detailed investigation on the subject is needed.

Tracking the future evolution of Gaia BH-like binaries shows that those systems may be the progenitors of LMXBs. Those LMXBs would be transient systems with the recurrence times of order of tens to over hundred years according to our calculations.

ACKNOWLEDGEMENTS

The co-author, Prof. Krzysztof Belczynski, passed away on 2024 January 13. He contributed to this paper with his ideas and discussions.

IK would like to thank Prof. Jean-Pierre Lasota for his comments that helped to improve the manuscript and to Ataru Tanikawa and Aleksey Generozov for helpful discussions.

We thank the referee for their constructive report that helped to improve the manuscript.

IK and KB were supported by the Polish National Science Center (NCN) grant Maestro (2018/30/A/ST9/00050). SB acknowledges funding for this work by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the project ‘The dynamics of stellar-mass black holes in dense stellar systems and their role in gravitational wave generation’ (project number 405620641; PI: S. Banerjee). Special thanks go to tens of thousands of citizen-science project ‘Universe@home’ (universeathome.pl) enthusiasts that help to develop StarTrack, the population synthesis code used in this study.

DATA AVAILABILITY

The data presented in this work can be made available based on the request to the corresponding author.

Footnotes

1

The peculiar velocity corresponds to systemic velocity in the model, as is explained in Section 2.1.

2

The mass accretion rate on compact object is assumed to be equal to the mass transfer rate during the stable Roche lobe overflow in StarTrack, but in general this is not true for the discs of transient LMXBs. However, it is a sufficiently good approximation for our purpose of the recurrence time estimation.

3

In the Ba21 runs, tidal energy loss and analytical tidal circularization were enabled (option |$27=1$|⁠; Aarseth 2003).

APPENDIX A:

Table A1.

The parameters used to calculate the outburst recurrence times at the onset and finish of LMXB phase.

	Onset	Finish
|$M_\mathrm{BH}$| (⁠|$\mathrm{{\rm M}_{\odot }}$|⁠)	9.7	10.2
|$\dot{M}_\mathrm{tr}$| (g s−1)	|$4.45\times 10^{18}$|	|$6.07\times 10^{18}$|
|$R_\mathrm{d}$| (cm)	|$9.61\times 10^{12}$|	|$2.73\times 10^{13}$|

APPENDIX B: GAIA BH1-LIKE BINARIES AS LISA SOURCES

We checked if Gaia BH1-like binaries become the Laser Interferometer Space Antenna (LISA) sources at some point of their late evolution using the parameters of a binary described in Section 2.2.2. We calculate the GW strain amplitude h taking into account that Gaia BH1 will become BH–WD binary according to our results. We start from the formula for the gravitational wave strain derived by Nelemans, Yungelson & Portegies Zwart (2001):

$$\begin{eqnarray} h=1.0\times 10^{-21}\frac{\sqrt{g(n,e)}}{n}\left(\frac{M_\mathrm{Ch}}{\mathrm{M}_\odot }\right)^{5/3}\left(\frac{P_\mathrm{orb}}{1\, \mathrm{hr}}\right)^{-2/3}\left(\frac{d}{1\, \mathrm{kpc}}\right)^{-1}, \\ \end{eqnarray}$$

(B1)

where |$g(n,e)$| is the Fourier decomposition of GW signal to n harmonics, given by the combination of Bessel functions of first kind (see Barack & Cutler 2004, equation 10), |$M_\mathrm{Ch}$| is a chirp mass in solar masses defined as |$M_\mathrm{Ch}=(M_\mathrm{BH}M_2)^{3/5}(M_\mathrm{BH}+M_2)^{-1/5}$|⁠, |$P_\mathrm{orb}$| is the orbital period in hours and d is the distance to the binary in kilo-parsecs.

For circular binaries the factor |$g(n,e)$| is non-zero only when |$n=2$| (⁠|$g(2,0)=1.0$|⁠), therefore, the frequency of GW emitted by circular binaries is always twice their orbital frequency. We calculate the GW strain |$h(n,e)$| for |$n=2$|⁠, |$e=0.0$|⁠, |$M_\mathrm{BH}=10.2\,$||$\mathrm{{\rm M}_{\odot }}$|⁠, |$M_2=0.5$||$\mathrm{{\rm M}_{\odot }}$|⁠, |$P_\mathrm{orb}=979.6$| d and |$d=0.48$| kpc, that is when Gaia BH1 becomes BH–WD binary:

$$\begin{eqnarray} h(2,0)=2.03\times 10^{-23} \, \, \, \, \, \, \, f_\mathrm{GW}=4.99\times 10^{-7} \mathrm{Hz} \end{eqnarray}$$

(B2)

Comparing the result with the expected LISA sensitivity (Robson, Cornish & Liu 2019) it is clear that Gaia BH1 will not contribute to LISA sources.

The Gaia BH2 parameters, when it becomes BH–WD binary, are comparable to Gaia BH1 except for the orbital separation which is |$\sim 2.5$| times wider. It is apparent that GW emitted from Gaia BH2 will be undetectable by LISA likewise, therefore we skip the calculations of GW strain for Gaia BH2.

Author notes

Deceased