Populations of stellar mass Black holes from binary systems

TL;DR

This paper models black hole populations from binary systems in galaxies, predicting their characteristics and abundance, and providing data for future observational studies and surveys.

Contribution

It introduces a realistic evolutionary model for black hole populations from binaries, including effects of initial parameters and metallicity, and offers a publicly available database.

Findings

Black hole populations are dominated by single black holes from binary disruptions.

Approximately 26 microlensing events per year are expected due to black holes in the Galactic Bulge.

The population characteristics are significantly influenced by initial parameters and metallicity.

Abstract

In large and complicated stellar systems like galaxies it is difficult to predict the number and characteristics of a black hole population. Such populations may be modelled as an aggregation of homogeneous (i.e. having uniform star formation history and the same initial chemical composition) stellar populations. Using realistic evolutionary models we predict the abundances and properties of black holes formed from binaries in these environments. We show that the black hole population will be dominated by single black holes originating from binary disruptions and stellar mergers. Furthermore, we discuss how black hole populations are influenced by such factors as initial parameters, metallicity, initial mass function, and natal kick models. As an example application of our results, we estimate that about 26 microlensing events to happen every year in the direction of the Galactic Bulge due to black holes in a survey like OGLE-IV. Our results may be used to perform in-depth studies related to realistic black hole populations, e.g. observational predictions for space survey missions like Gaia, or Einstein Probe. We prepared a publicly available database with the raw data from our simulations to be used for more in-depth studies.

Tables

Table 1. Table 1 Summary of models

Model	difference in respect to standard model
STD	standard (reference) model:
	solar metallicity $Z_{\odot }$
	distribution of initial periods
	$P\left(\log P\right)\sim {\left(\log P\right)}^{-0.55}$
	distribution of initial eccentricities
	$P\left(e\right)\sim e^{-0.42}$
	BH/NS natal kicks are drawn from
	Maxwellian distribution with $\sigma =265$ km s-1
	BH natal kicks reduced due to fallback
	moderate slope for high-mass end of the IMF
	$\Gamma =-2.3$
mid-$Z$	metallicity equal to $10\%Z_{\odot }$
low-$Z$	metallicity equal to $1\%Z_{\odot }$
SS0	distribution of initial separations $\log \left(a\right)\sim 1$
	distribution of initial eccentricities $P\left(e\right)\sim e$
NKR	BH NKs are inversely proportional to the BH’s mass
	(Eq. 2)
NKBE	BH/NS natal kick proportional to
	ratio of ejecta mass and remnant mass (Eq. 3)
flat IMF	flat slope for high-mass end of the IMF ($\Gamma =-1.9$)
steep IMF	steep slope for high-mass end of the IMF ($\Gamma =-2.7$)

Table 2. Table 2 Number of BHs in the Milky Way Equivalent Galaxy

Model	$N_{\mathrm{BH},\mathrm{tot}}$	$N_{\mathrm{dSBH}}$		$N_{\mathrm{BHB}}$		$N_{\mathrm{BH},\mathrm{DCO}}$		$N_{\mathrm{mSBH}}$
STD	$1.1\times {10}^8$	$5.4\times {10}^7$	($49.7\%$)	$1.5\times {10}^6$	($1.4\%$)	$5.1\times {10}^6$	($4.7\%$)	$4.8\times {10}^7$	($44.2\%$)
mid-$Z$	$1.1\times {10}^8$	$3.4\times {10}^7$	($30.6\%$)	$2.0\times {10}^6$	($1.8\%$)	$1.2\times {10}^7$	($10.8\%$)	$6.3\times {10}^7$	($56.8\%$)
low-$Z$	$1.2\times {10}^8$	$3.7\times {10}^7$	($30.0\%$)	$1.5\times {10}^6$	($1.2\%$)	$1.9\times {10}^7$	($15.4\%$)	$6.6\times {10}^7$	($53.4\%$)
SS0	$1.1\times {10}^8$	$6.6\times {10}^7$	($58.8\%$)	$3.4\times {10}^6$	($3.0\%$)	$6.9\times {10}^6$	($6.1\%$)	$3.6\times {10}^7$	($32.1\%$)
NKR	$1.1\times {10}^8$	$6.4\times {10}^7$	($58.1\%$)	$5.9\times {10}^4$	($0.1\%$)	$1.1\times {10}^4$	($0.0\%$)	$4.6\times {10}^7$	($41.8\%$)
NKBE	$1.1\times {10}^8$	$6.5\times {10}^7$	($58.5\%$)	$5.4\times {10}^4$	($0.0\%$)	$2.7\times {10}^3$	($0.0\%$)	$4.6\times {10}^7$	($41.4\%$)
flat IMF	$2.7\times {10}^8$	$1.4\times {10}^8$	($52.6\%$)	$3.2\times {10}^6$	($1.2\%$)	$1.3\times {10}^7$	($4.9\%$)	$1.1\times {10}^8$	($41.3\%$)
steep IMF	$3.4\times {10}^7$	$1.6\times {10}^7$	($46.8\%$)	$5.6\times {10}^5$	($1.6\%$)	$1.6\times {10}^6$	($4.7\%$)	$1.6\times {10}^7$	($46.8\%$)

Table 3. Table 3 Typical formation routes of BHs

Route	evolutionary route
Single BHs from disrupted binaries
${ℛ}_{\mathrm{dSBH},1}$	MT1(1/2/4/5-1/2/4/5) SN1 Disruption SN2
${ℛ}_{\mathrm{dSBH},2}$	SN1 SN2 Disruption
BH XRBs
${ℛ}_{\mathrm{MTBHB},1}$	CE1(4/5-1;7/8-1) SN1 MT2(14-1/2/3/4/5/6)
${ℛ}_{\mathrm{MTBHB},2}$	CE1(4/5-1;7/8-1) SN1 MT2(13-1/2/3) AICBH1
	MT2(14-1/2/3)
DCOs
${ℛ}_{\mathrm{DCO},1}$	SN1 SN2
${ℛ}_{\mathrm{DCO},2}$	MT1(1/2/4/5-1/2/4) SN1
	CE2(13/14-4/5;13/14-7/8) SN2
Mergers
${ℛ}_{\mathrm{mSBH}}$	CE1(2/3/4-1;7-1) MT1(7-1) Merger

Table 4. Table 4 Single BHs from disrupted binaries

Model	Route	Number	$a_{\mathrm{ZAMS}}\left[R_{\odot }\right]$
STD	${ℛ}_{\mathrm{dSBH},1}$	$3.6\times {10}^7\left(67\%\right)$	$28÷400$
	${ℛ}_{\mathrm{dSBH},2}$	$1.7\times {10}^7\left(32\%\right)$	$>4400$
mid-$Z$	${ℛ}_{\mathrm{dSBH},1}$	$2.0\times {10}^7\left(58\%\right)$	$20÷800$
	${ℛ}_{\mathrm{dSBH},2}$	$1.2\times {10}^7\left(33\%\right)$	$>7700$
low-$Z$	${ℛ}_{\mathrm{dSBH},1}$	$1.5\times {10}^7\left(40\%\right)$	$17÷180$
	${ℛ}_{\mathrm{dSBH},2}$	$2.0\times {10}^7\left(55\%\right)$	$>1600$
SS0	${ℛ}_{\mathrm{dSBH},1}$	$4.0\times {10}^7\left(61\%\right)$	$28÷1000$
	${ℛ}_{\mathrm{dSBH},2}$	$2.4\times {10}^7\left(36\%\right)$	$>6200$
NKR	${ℛ}_{\mathrm{dSBH},1}$	$4.0\times {10}^7\left(62\%\right)$	$28÷400$
	${ℛ}_{\mathrm{dSBH},2}$	$2.2\times {10}^7\left(35\%\right)$	$>5200$
NKBE	${ℛ}_{\mathrm{dSBH},1}$	$4.1\times {10}^7\left(63\%\right)$	$28÷400$
	${ℛ}_{\mathrm{dSBH},2}$	$2.2\times {10}^7\left(34\%\right)$	$>5200$
flat IMF	${ℛ}_{\mathrm{dSBH},1}$	$9.4\times {10}^7\left(66\%\right)$	$28÷430$
	${ℛ}_{\mathrm{dSBH},2}$	$4.7\times {10}^7\left(33\%\right)$	$>5400$
steep IMF	${ℛ}_{\mathrm{dSBH},1}$	$1.1\times {10}^7\left(68\%\right)$	$28÷400$
	${ℛ}_{\mathrm{dSBH},2}$	$5.1\times {10}^6\left(31\%\right)$	$>5000$

Table 5. Table 5 Double compact objects containing black holes

Model	$N_{\mathrm{BH},\mathrm{DCO}}$	BH+BH		${ℛ}_{\mathrm{DCO},1}$		${ℛ}_{\mathrm{DCO},2}$
STD	$5.1\times {10}^6$	$4.9\times {10}^6$	(95.9%)	$2.4\times {10}^6$	(90.1%)	$1.6\times {10}^4$	(0.6%)	$2.5\times {10}^4$	(0.9%)	$1.7\times {10}^4$	(0.6%)
mid-$Z$	$1.2\times {10}^7$	$1.2\times {10}^7$	(96.6%)	$2.6\times {10}^6$	(42.6%)	$4.0\times {10}^5$	(3.4%)	$1.7\times {10}^6$	(27.3%)	$8.1\times {10}^5$	(13.0%)
low-$Z$	$1.9\times {10}^7$	$1.9\times {10}^7$	(97.5%)	$5.7\times {10}^6$	(57.8%)	$1.8\times {10}^6$	(18.1%)	$2.7\times {10}^6$	(27.0%)	$1.9\times {10}^6$	(19.0%)
SS0	$6.9\times {10}^6$	$6.6\times {10}^6$	(95.6%)	$3.4\times {10}^6$	(93.1%)	$2.4\times {10}^4$	(0.7%)	$3.0\times {10}^4$	(0.8%)	$2.3\times {10}^4$	(0.6%)
NKR	$1.1\times {10}^4$	$5.1\times {10}^3$	(45.5%)	$8.8\times {10}^2$	(10.3%)	$1.6\times {10}^3$	(19.1%)	$6.8\times {10}^3$	(79.4%)	$3.2\times {10}^3$	(36.8%)
NKBE	$2.7\times {10}^3$	$2.5\times {10}^2$	(9.5%)	–		–		$2.5\times {10}^3$	(100%)	$1.3\times {10}^2$	(5.0%)
flat IMF	$1.3\times {10}^7$	$1.3\times {10}^7$	(96.6%)	$6.2\times {10}^6$	(90.6%)	$3.8\times {10}^4$	(0.6%)	$6.0\times {10}^4$	(0.9%)	$3.8\times {10}^4$	(0.6%)
steep IMF	$1.6\times {10}^6$	$1.5\times {10}^6$	(95.3%)	$7.7\times {10}^5$	(90.7%)	$3.7\times {10}^3$	(0.4%)	$5.7\times {10}^3$	(0.7%)	$3.7\times {10}^3$	(0.4%)

Table 6. Table 6 Wide and mass-transferring binaries

Model	Wide BHB			MTBHB
	$N$	$a\left[R_{\odot }\right]$	Typ. comp.	$N$	$a\left[R_{\odot }\right]$	Typ. comp.
STD	$1.5\times {10}^6$	$\gtrsim 10,000$	CO WD ($86\%$)	$8.5\times {10}^1$	$\lesssim 115$	MS ($75\%$)
mid-$Z$	$7.5\times {10}^5$	$\gtrsim 8,300$	CO WD ($79\%$)	$4.2\times {10}^3$	$\lesssim 110$	MS ($85\%$)
low-$Z$	$1.2\times {10}^6$	$\gtrsim 5,800$	CO WD ($85\%$)	$4.4\times {10}^2$	$\lesssim 84$	MS ($91\%$)
SS0	$3.3\times {10}^6$	$\gtrsim 8,300$	CO WD ($74\%$)	$5.7\times {10}^2$	$\lesssim 108$	MS ($73\%$)
NKR	$3.1\times {10}^3$	$\gtrsim 4,800$	CO WD ($89\%$)	$1.7\times {10}^2$	$\lesssim 90$	MS ($86\%$)
NKBE	$1.0\times {10}^{-4}$	$\gtrsim 14,000$	CHeB ($100\%$)	$1.8\times {10}^2$	$\lesssim 75$	MS ($82\%$)
flat IMF	$3.1\times {10}^6$	$\gtrsim 10,000$	CO WD ($86\%$)	$1.3\times {10}^2$	$\lesssim 113$	MS ($73\%$)
steep IMF	$5.5\times {10}^5$	$\gtrsim 10,000$	CO WD ($86\%$)	$3.2\times {10}^1$	$\lesssim 116$	MS ($70\%$)

Table 7. Table 7 Number of BHs in the Milky Way Equivalent Galaxy (model A)

Model	$N_{\mathrm{BH},\mathrm{tot}}$	$N_{\mathrm{dSBH}}$		$N_{\mathrm{BHB}}$		$N_{\mathrm{BH},\mathrm{DCO}}$		$N_{\mathrm{mSBH}}$
STD	$1.0\times {10}^8$	$5.4\times {10}^7$	($54.2\%$)	$1.5\times {10}^6$	($1.5\%$)	$6.1\times {10}^6$	($6.1\%$)	$3.8\times {10}^7$	($38.2\%$)
mid-$Z$	$1.2\times {10}^8$	$3.5\times {10}^7$	($29.9\%$)	$2.0\times {10}^6$	($1.7\%$)	$2.2\times {10}^7$	($18.8\%$)	$5.8\times {10}^7$	($49.6\%$)
low-$Z$	$1.2\times {10}^8$	$3.7\times {10}^7$	($30.7\%$)	$1.5\times {10}^6$	($1.2\%$)	$2.1\times {10}^7$	($17.4\%$)	$6.1\times {10}^7$	($50.6\%$)
SS0	$1.0\times {10}^8$	$6.6\times {10}^7$	($65.0\%$)	$3.4\times {10}^6$	($3.3\%$)	$8.1\times {10}^6$	($8.0\%$)	$2.4\times {10}^7$	($23.6\%$)
NKR	$1.0\times {10}^8$	$6.4\times {10}^7$	($63.8\%$)	$6.0\times {10}^4$	($0.1\%$)	$2.1\times {10}^5$	($0.2\%$)	$3.6\times {10}^7$	($35.9\%$)
NKBE	$1.0\times {10}^8$	$6.5\times {10}^7$	($64.3\%$)	$5.5\times {10}^4$	($0.1\%$)	$8.7\times {10}^4$	($0.1\%$)	$3.6\times {10}^7$	($35.6\%$)
flat IMF	$2.5\times {10}^8$	$1.4\times {10}^8$	($56.0\%$)	$3.2\times {10}^6$	($1.3\%$)	$1.6\times {10}^7$	($6.4\%$)	$9.1\times {10}^7$	($36.4\%$)
steep IMF	$3.1\times {10}^7$	$1.6\times {10}^7$	($50.9\%$)	$5.6\times {10}^5$	($1.8\%$)	$1.9\times {10}^6$	($6.0\%$)	$1.3\times {10}^7$	($41.3\%$)

Equations

$$V_{\mathrm{kick}}=V_{\mathrm{kick},\mathrm{Maxwell}}(1-f_{\mathrm{FB}}),$$

$$\begin{array}[]{lcl}V_{\rm kick,NS}&=&V_{\rm kick,Maxwell},\\ V_{\rm kick,BH}&=&V_{\rm kick,Maxwell}\frac{M_{\mathrm{max,NS}}}{M_{\rm BH}},\end{array}$$

$$V_{\mathrm{kick}}=60\frac{\mathrm{km}}{\mathrm{s}}\left(\frac{M_{\mathrm{ejecta}}}{M_{NS/BH}}\right)+130\frac{\mathrm{km}}{\mathrm{s}},$$

$$f=\frac{{\int }_{M_{\mathrm{NS},\mathrm{ZAMS},\max }}^{150M_{\odot }}{P}_{\mathrm{BH}+\mathrm{NS}}(M_1)\mathrm{IMF}(M_1)d{M}_1}{{\int }_{M_{\mathrm{NS},\mathrm{ZAMS},\max }}^{150M_{\odot }}{P}_{\mathrm{BH}+\mathrm{BH}}(M_1)\mathrm{IMF}(M_1)d{M}_1},$$

$$P_{\mathrm{BH}+\mathrm{NS}}(M_1)=\frac{M_{\mathrm{NS},\mathrm{ZAMS},\max }-M_{\mathrm{NS},\mathrm{ZAMS},\min }}{M_1-M_{\mathrm{NS},\mathrm{ZAMS},\min }}$$

$$P_{\mathrm{BH}+\mathrm{BH}}(M_1)=\frac{M_1-M_{\mathrm{NS},\mathrm{ZAMS},\max }}{M_1-M_{\mathrm{NS},\mathrm{ZAMS},\min }},$$

$${\rho }_{\mathrm{star}}\propto e^{\frac{|z|}{0.3\mathrm{kpc}}-\frac{R}{2.6\mathrm{kpc}}},$$

$${\theta }_{\mathrm{E}}=\sqrt{\kappa M_{\mathrm{L}}{\pi }_{\mathrm{rel}}},{\pi }_{\mathrm{rel}}=\frac{1}{D_{\mathrm{L}}}-\frac{1}{D_{\mathrm{S}}},$$

$${\mu }_{\mathrm{rel},y/z}=\frac{v_{\mathrm{rel},y/z}}{x_{\odot }-x},$$

$$\begin{gathered} v_{\mathrm{rel},y/z}=v_{\mathrm{BH},y/z}-v_{\mathrm{Earth},y/z} \\ -(v_{\mathrm{source},y/z}-v_{\mathrm{Earth},y/z})(1-\frac{x}{x_{\odot }}), \end{gathered}$$

$$\mathrm{E}(N_{\mathrm{ML}})=\frac{{\sum }_i(2{\theta }_{\mathrm{E},\mathrm{i}}\frac{v_{\mathrm{rel},\mathrm{i}}}{x_{\text{sun}}-x_{\mathrm{i}}}t+\pi {\theta }_{\mathrm{E},\mathrm{i}}^2)}{{\Omega }_{\mathrm{bulge}}}{N}_{\star ,bulge},$$

Full text

1. Introduction

A BH is defined as a region in space from which nothing, even light, can escape (for a recent review see Bambi, 2017). BHs may be detected when interacting with other objects (e.g. in X-ray binaries; Remillard & McClintock, 2006; Corral-Santana et al., 2016), as gravitational wave sources (e.g., Abbott et al., 2016a, 2017a), when they interfere with radiation (e.g. as gravitational lenses; Wyrzykowski et al., 2016), or in non-interacting binaries by observations of their companions (e.g. Thompson et al., 2018; Mashian & Loeb, 2017). It was also proposed that SBHs may be observed as X-ray novae (e.g. Matsumoto et al., 2018).

The so called stellar-mass BHs, i.e. stellar-origin BHs and BHs originating from mergers of stars and/or stellar-origin compact objects, with masses from $\sim 5$ to possibly a few $\times 100\,M_{\odot}$ are final stages of massive stars evolution (e.g. Neugebauer, 2003) and are distinguished from other subgroups such as super-massive BHs ($M_{\rm BH}\gtrsim 10^{6}\,M_{\odot}$; e.g. Ferrarese & Ford, 2005), intermediate-mass BHs ($\sim 10^{2}<M_{\rm BH}<\sim 10^{5}\,M_{\odot}$; e.g. Mezcua, 2017; Koliopanos, 2018) and hypothetical primordial BHs (e.g. Chapline, 1975; Khlopov, 2010). Hereafter, we focus exclusively on stellar-mass BHs.

Corral-Santana et al. (2016) provided a list of 59 BHs in transient X-ray binaries (XRB). 22 of them have dynamically confirmed mass estimates (see Casares et al., 2017, for a recent list). 5 BHs were detected in High-mass XRBs (defined as heaving a donor mass above $10\,M_{\odot}$). Recently, gravitational waves made it possible to discover first double BH merger events (Abbott et al., 2016a, b, 2017a, 2017b), and a double NS merger GW170817, which may have formed a low-mass BH (Pooley et al., 2017). Interestingly, Adams et al. (2017) observed a $\sim 25\,M_{\odot}$ star to disappear after a short brightening, which may be interpreted as a BH-formation event. Up to now, none SBHs were confirmed (e.g. Tsuna et al., 2018).

Previous studies of BH populations are usually outdated and do not take into account recent progress in our understanding of massive star evolution (e.g., Langer, 2012; Vink et al., 2015; Vink, 2015). The earliest estimations predicted $\sim 100$ million BHs in the Milky Way (MW) galaxy (Shapiro & Teukolsky, 1983) with as much as $45\%$ of supernovae occurring in close binaries (Tutukov et al., 1992). Timmes et al. (1996) provided an upper limit of $1.4\times 10^{9}$ and Samland (1998) estimated $1.8\times 10^{8}$ BHs in the galaxy taking into account changes in star formation rate (SFR). Belczynski et al. (2002a, 2004) used the StarTrack population synthesis code in its earlier version to investigate BHs formed in star formation (SF) bursts. Their results indicate that most of BHs in the MW are actually single objects (not in binaries), whereas these which remain bound will have mostly main-sequence (MS) companions, provided that the stellar population is not very old.

Recent studies do not provide the estimations for parameter distributions of the entire BH population. For example, Elbert et al. (2018) predicted $\sim 10^{8}$ BHs in a MW-type galaxy, however, their approach does not take into account binary interactions directly and focuses on binary BHs (BH+BH) only. Lamberts et al. (2018) performed a detailed cosmological simulation of the MW evolution (including e.g. changes in metallicity) and predicted, that $\sim 1$ million of BH binaries have already merged in the galaxy, whereas $\sim 3$ millions are still present. However, their study does not account for BHs with lower-mass companions, or disrupted systems, which may actually constitute a bulk of the BH population (Belczynski et al., 2004).

Although the current population of BHs is predicted to be mainly single, their progenitors did not have to be single stars. Due to a NK (e.g. Lyne & Lorimer, 1994; Janka, 2013), a binary (especially with large orbital separation, or a low-mass companion) may be disrupted (Iben & Tutukov, 1997). Recent observational studies (see e.g. Sana, 2017, for a review) suggest that even $\gtrsim 90\%$ of massive stars (BH progenitors) are born in binaries (or multiple systems). The presence of a companion may significantly affect the evolution of a progenitor and, therefore, final properties of a BH. Some studies used a simplified approach. E.g., Elbert et al. (2018) encapsulated all the interactions in just two parameters. In this study we focused on a detailed consideration of binary interactions, however, we neglected any higher-order systems (e.g. Toonen et al., 2016).

A serious problem in obtaining estimations of BH populations is the lack of knowledge about the stellar environment. Taking the MW as an example, not even the total stellar mass is well constrained, let alone the SFH, or the distribution and motions of stars. What is more, a galaxy with complicated evolution and structure cannot be modelled with a uniform population of objects evolved from a single SF burst. Therefore, in our approach we take into account a range of evolutionary models with different initial parameters. Such results can be joined together according to SFH and chemical evolution of a stellar population in order to obtain more realistic distributions.

Specifically, we propose a step-by-step approach in which we use results from modelling of homogeneous and scalable BH populations in order to build a complicated (non-homogeneous) stellar systems like galaxies. We provide the estimates for a simplified MW galaxy as an exemple (for detailed modeling of the MW galaxy see Olejak et al. in prep.). However, our results may be used to simulate virtually any large stellar population providing the dynamical interactions between stars may be neglected (it is not true e.g. for globular clusters).

The two most important factors that influence BH populations are the metallicity, and NKs. The metallicity of a star significantly affects the mass loss in stellar wind and, therefore, the final compact object’s mass (Fryer & Kalogera, 2001). Belczynski et al. (2010a) showed that the maximal BH mass for solar metallicity environment is $\sim 15\,M_{\odot}$, whereas for lower metallicity environments may reach $\sim 80\,M_{\odot}$. However, their study did not consider the pair-instability supernova and pair-instability pulsation supernova (e.g. Woosley, 2017, 2019) and was performed for single-stars only. Both, pair-instability supernovae and pair-instability pulsation supernova are believed to produce the second ”mass gap” in the distribution of BH masses between $\sim 50\,\hbox{--}\,135\,M_{\odot}$ (Belczynski et al., 2016c; Spera & Mapelli, 2017; Marchant et al., 2018; Woosley, 2019). However, in binary systems the masses of BHs can fill this range due to mass transfer (MT) phases and mergers (e.g. Spera et al., 2019).

Uncertainties in the modelling of common envelope (CE) phase add additional significant error to the predictions of binary evolution (for a review see Ivanova et al., 2013). The CE seems to be essential for the formation of XRBs where a strong reduction of initial separation is necessary to produce a Roche lobe overflow (RLOF), or significant accretion from wind. However, it was shown (Wiktorowicz et al., 2014), that different prescriptions for the CE phase, give similar results with respect to the population of XRBs. In contrast, the formation of close double compact objects (DCO), with time-to-merger smaller than the Hubble time, which also typically involves the CE phase, is highly influenced by the adopted CE model (e.g. Dominik et al., 2012). Nonetheless, as we show in this paper, only a small part of the total BH population resides in XRBs and close DCOs. Most stellar mergers of isolated binaries occur as a result of failed CE ejection, therefore, CE model may strongly influence the population of BHs originating form merger products. However, even though outcomes of mergers may constitute a significant fraction of a BH population, the poorly understood physics of stellar mergers and post-merger evolution make it impossible to reliably describe the population of these objects. In this paper, therefore, we do not include the analysis of CE models.

The fraction of binaries that remain bound after a supernova (SN) explosion as well as peculiar velocities of BHs are particularly sensitive to the assumptions about the NKs. The former impacts the predicted ratio of the number of SBHs to those found in binaries. The latter influences the spacial distribution of BHs in stellar systems (e.g. galaxies). Although the SN mechanism is not well understood, NKs are usually being connected to asymmetries arising in this process. Main candidates are: mass ejection (e.g. Wongwathanarat et al., 2013), gravitational waves (Bonnell & Pringle, 1995), and neutrinos (e.g. Fryer & Kusenko, 2006). Miller-Jones (2014) derived peculiar velocities of several BH binaries (BHB, i.e. binaries composed of a BH and a non-compact companion) obtaining values between $19\div 144\,\mathrm{km}\,\mathrm{s}^{-1}$. However, it is usually difficult to assess NKs from present-day motions (Fragos et al., 2009; Repetto et al., 2012; Repetto & Nelemans, 2015). Recent analysis performed by Repetto et al. (2017) has shown that at least some of the BHs in the MW should have obtained high NKs ($\sim 100$ km/s), comparable to these of NSs. On the other hand, Mandel (2016) showed that velocities higher than $\sim 80$ km/s are not necessary, although cannot be ruled out. Additionally, Jonker & Nelemans (2004) showed that XRBs with BH accretors have similar spacial distribution as theses with NS accretors, what suggests a similar NK magnitude at birth. See Belczynski et al. (2016b) and references therein for a recent discussion on the BH NKs.

In this paper, we performed the first simulation of the entire BH population for different evolutionary models. Our analysis is focused on SBHs from binary disruption events (Sec. 3.1), BHBs, DCOs (Sec. 3.2), and potential BHs originating from stellar mergers (Sec. 3.3). We pay a particular attention to the description of the results for the standard model (STD), which may depict the MW disk field population, and differences in relation to other models.

The main goals of this work are: 1) to describe the general characteristics of the BH population and its dependence on the most important model parameters (e.g., metallicity, NKs), and 2) to provide a reference point for forthcoming in-depth studies focusing on astrophysical problems involving BH populations like: XRBs, microlensing, gravitational wave sources, to name a few. Especially educating are possible comparison studies with the results of actual and future surveys focused on gravitational wave sources (e.g., aLIGO, Einstein Telescope, LISA), X-ray sources (e.g., XMM-Newton, Chandra, NuSTAR), or gravitational microlensing (OGLE, Gaia, LSST), which may result in obtaining better constrains for the evolutionary models.

Particularly, in Sec. 4 we estimate a number of microlensing events toward the Galactic bulge. The microlensing method seems to be a promising way to search for SBHs, which, as we show in this work, constitute a vast majority of all BHs in stellar populations.

In order to allow for easier and more flexible usage of the data introduced by this paper, we prepared a publicly available database where all the results are available for download:

https://universeathome.pl/universe/bhdb.php

2. Methodology

We utilize a recent version of the StarTrack population synthesis code (Belczynski et al., 2002b, 2008, with further updates). The code has been frequently used to study BHs in XRBs and DCOs (e.g. Dominik et al., 2012, 2013, 2015; Wiktorowicz et al., 2014; Belczynski et al., 2016c; Klencki et al., 2017). Recent updates include, but are not limited to, new prescriptions for wind mass loss from massive stars (Vink, 2011), pair-instability supernovae and pair-instability pulsation supernovae (Belczynski et al., 2016a; Woosley, 2017).

Here, we analyse 8 main models (Tab. 1) which differ in parameters that significantly affect the resulting BH population. Additionally, in the web database, we provide a grid of all 54 models which allow to investigate combinations of the main models. For each model we have simulated the evolution of $2\times 10^{6}$ binaries from zero-age main-sequence (ZAMS) to disruption, merger, or reaching the age of $15\,\mathrm{Gyr}$. For mergers that are not DCOs we apply a simple formalism to estimate the endpoint of post-merger single-star evolution (see Sec. 2.2).

Initial primary masses ($M_{\rm a}$) are drawn from a broken power-law distribution (Kroupa et al., 1993) with $\Gamma=-1.3$ for $M_{\rm a}<0.5\,M_{\odot}$ and $\Gamma=-2.2$ for $0.5<M_{\rm a}<1.0\,M_{\odot}$. For initial masses $M_{\rm a}>1.0\,M_{\odot}$ we chose $\Gamma=-2.3$ (STD; Kroupa, 2001), $-2.7$ (steep IMF model; Kroupa & Weidner, 2003), or $-1.9$ (flat IMF model; Schneider et al., 2018). Conclusions of Schneider et al. (2018) were recently revised by Farr & Mandel (2018) who argued in favor of a much steeper $\Gamma\approx-2.11$, or $-2.15$. Being aware of this, we leave the original value in order to emphasize the influence of IMF steepness on the population of BHs. The mass ratio $q=M_{\rm b}/M_{\rm a}$ (where $M_{\rm b}$ is the mass of the secondary – initially less massive star) is drawn from a uniform distribution between $0.08\,M_{\odot}/M_{\rm a}$ and $1$. In STD model, we assume that the initial distribution of periods is $P(\log P)\sim(\log P)^{-0.55}$, and distribution of eccentricities is $P(e)\sim e^{-0.42}$ (Sana et al., 2012). Even though, the results of Sana et al. (2012) are for stars with initial masses between $15\div 60\,M_{\odot}$, we extrapolate them to entire investigated range ($0.08\div 150\,M_{\odot}$). Additionally, in model SS0, we test a distribution of initial separations that is flat in logarithm (Abt, 1983) with maximal value set as $10^{5}\,R_{\odot}$ and a thermal distribution of eccentricities ($P(e)\sim e$; Duquennoy & Mayor, 1991).

The effects of pair-instability supernovae and pair-instability pulsation supernovae (Woosley, 2017) may significantly alter the final BH mass (e.g. Belczynski et al., 2016a). The instability leads to a significant mass loss for massive stars (Woosley et al., 2007), or disruption of the entire star (Heger & Woosley, 2002). Following (Belczynski et al., 2016a), we assume that for stars that form massive helium cores ($M_{\rm He}=45\,\hbox{--}\,65\,M_{\odot}$) all the envelope above the inner $45\,M_{\odot}$ is lost due to pulsations. Furthermore, we assume that stars which form heavier helium cores ($M_{\rm He}=65\,\hbox{--}\,135\,M_{\odot}$) are subject to the pair-instability supernovae and leave no remnant.

The CE phase is important for the evolution of many binaries and may lead to the formation of a much closer system, or a merger. In a general situation, we utilize the simple energy balance (Webbink, 1984). However, the donor type influences significantly the binary survival during this phase (e.g. Belczynski et al., 2007). MS and Hertzsprung gap (HG) donors lack a clear core-envelope structure (no clear entropy jump; e.g. Taam & Sandquist, 2000; Ivanova & Taam, 2004), so the orbital energy is transferred to the entire star, instead of being transferred to the envelope only, what prevents envelop ejection and leads to a merger. It is not clear when the core-envelop boundary emerges (late HG, or red giant phase), thus two models for the treatment of HG donors in CE were introduced (e.g. Dominik et al., 2012). First one (model B), assumes that a CE phases with a HG donor always results in a merger. Second one (model A), producing significantly higher DCO merger rates (e.g. Belczynski et al., 2007; Dominik et al., 2012) allows for survival in such situations. We adopt model B in the present study and analyse the influence of model A in Sec. 5.4.

We test models with 3 different metallicities. Most of the models have solar metallicity ($Z=0.02=\,Z_{\odot}$, where $\,Z_{\odot}$ is the solar metallicity; Villante et al., 2014). We also introduce two models with lower metallicity: mid-$Z$ model ($Z=0.002=10\%\,Z_{\odot}$), and low-$Z$ model ($Z=0.0002=1\%\,Z_{\odot}$). Additionally, three different NK models are included in our simulations. Those models are described below.

2.1. Natal kick models

Hobbs et al. (2005) performed a study of 233 galactic pulsars proper motions and find out that for young (i.e. with characteristic age $\tau_{\rm c}=P/2\dot{P}<3\,\mathrm{yr}$) pulsars the velocity distribution is well fitted with Maxwellian distribution with $\sigma=265$ km/s. Velocities of such a population may quite well resemble the NKs of NSs. However, it is not obvious, if the NKs of BHs follow the same distribution. Especially, if NKs are driven by the asymmetries in ejecta, the post-SN fall-back may significantly lower the BH’s velocity. Proper motion measurements are available only for a few BHs residing in XRBs (Miller-Jones, 2014) and the connection between the current peculiar velocity and the NK is not straightforward (e.g. Fragos et al., 2009). On the other hand, the analysis of the positions of XRBs harbouring BHs may suggest that distributions of NKs of BHs and NSs are similar (e.g., Jonker & Nelemans, 2004; Repetto et al., 2017, however, see Mandel 2016).

In the standard model, we lower the remnant’s NK proportionally to the fall-back of material after a SN explosion. Precisely, the NK for a compact object is equal

[TABLE]

where $V_{\rm kick,Maxwell}$ has a Maxwellian distribution with $\sigma=265$ km/s and $f_{\rm FB}$ indicates what fraction of ejected mass falls back onto a compact object (see Fryer et al., 2012). This model produces both low-velocity and high-velocity BHs which is in agreement with observational estimates for XRBs (e.g. Belczynski et al., 2016b).

If NKs are driven by neutrino based mechanism, it is predicted that the imposed momentum will be the same for BHs and NSs (e.g., Janka, 2013; Rodriguez et al., 2016). Such NKs will be inversely proportional to the mass of a BH ($V_{\rm kick}\approx M_{\rm BH}^{-1}$). In NKR model we assume that NSs’ NK distribution has a Maxwellian shape with $\sigma=265$ km/s, but BHs NK distribution is lowered. Precisely,

[TABLE]

where $M_{\mathrm{max,NS}}=2.5\,M_{\odot}$ is an adopted limit for a NS mass above which it collapses to a BH. It must be noted that, even for low-mass BHs ($\sim 5$–$7\,M_{\odot}$), such a prescription will produce significantly lower NKs, than the NKs of NSs, what stands in contradiction with the results of Repetto et al. (e.g. 2017). Nevertheless, we included this model as a parameters study.

Finally, we incorporated a simple model proposed by Bray & Eldridge (2018), who suggested that the NK distribution can be described as a linear function of a ratio between the mass of the ejecta ($M_{\rm ejecta}$) and the mass of the remnant ($M_{\rm remnant}$). We applied (model NKBE) their fit to the observations provided by Hobbs et al. (2005) according to which the NK in 3D is given by

[TABLE]

where $M_{\rm ejecta}$ is mass ejected from the star during SN and $M_{\rm NS/BH}$ is mass of a resulting NS, or a BH. Although Bray & Eldridge (2018) provided their fit only for NSs, we apply it also to BHs as a parameter study.

2.2. Mergers of binary components

Stellar mergers may constitute a significant fraction of massive stars (e.g. Langer, 2012), and therefore, progenitors of BHs. Head-on collisions may play a significant role in dense stellar systems like globular clusters (e.g. Glebbeek et al., 2013), however, stellar mergers in field populations are more probable due to orbital angular momentum loss (e.g. during CE phase). Podsiadlowski et al. (1992) predicted that $\sim 10\%$ of $8<M<20\,M_{\odot}$ primaries merge with their companions before a SN and de Mink et al. (2014) showed that $\sim 8^{+9}_{-4}\%$ of massive single stars may actually be products of binary mergers. Objects like Red Novae, or Blue Stragglers are thought to be observed during, or after the merging process (e.g., Blagorodnova et al., 2017; Leonard, 1989; Kochanek et al., 2014).

In case of non-DCO mergers (including mergers of a compact object with a non-compact object), a merger process is not well understood. Especially, the amount of mass ejected from the system, which is an important factor for post-merger evolution, is poorly constrained. A product of a merger of two MS stars may evolve similarly to a star which was always single, but if evolved stars are involved, the post-merger evolution is more complicated (Glebbeek et al., 2013). Especially, a merger outcome may evolve unlike any single star (Vigna-Gómez et al., 2019).

We assume that non-DCO mergers occur in the following situations:

• •

failed envelope ejection during a CE event (e.g. Justham et al., 2014),

• •

donor in a CE phase haven’t developed a clear boundary between core and the envelope (e.g. MS and HG donors; Ivanova & Taam, 2004),

• •

donor’s radius exceeds two times its Roche lobe radius ($R_{\rm donor}>2R_{\rm donor,RL}$) during RLOF.

Due to the fact that the merger physics and post-merger evolution is poorly understood, we provide information on the binary parameters just before the merger to allow for different approaches to this conundrum.

As far as DCO mergers are concerned, many studies were devoted to an in-depth analysis of this phenomenon in stellar populations (e.g., Lipunov et al., 1997a; Dominik et al., 2012, 2013, 2015; Mennekens & Vanbeveren, 2014; Mandel & de Mink, 2016; Belczynski et al., 2016a). DCO mergers affect the BH populations, especially the BH mass distribution. As a part of our results, in the public database we provide information about DCO formation and estimate the time to merger using the formula of Peters (1964). In order to give predictions and estimate the fraction of BHs which come from DCO mergers, we assumed that the mass loss during a DCO merger is negligible.

In this work, in order to estimate the population of SBH originating from binary mergers, we adopt a simplistic approach of Olejak et al. (in prep.). We include only main channels responsible for $>95\%$ of all the mergers. Specifically, for different evolutionary types of stars we use the following procedure:

• •

MS+MS - Outcome is a MS star. We assume that half of the mass of the lighter component is lost in the process.

• •

MS+HeS - Outcome is a helium star (HeS) star. We assume that half of the mass of the MS component is lost in the process.

• •

HeS+HeS - Outcome is a HeS star. We assume that half of the mass of the lighter component is lost in the process.

• •

NS+MS/HeS - We assume that half of the mass of the MS/HeS star is lost in the process and that it becomes a NS/BH if its final mass is lower/higher than $M_{\rm max,NS}=2.5\,M_{\odot}$.

• •

BH+MS/HeS - We assume that half of the mass of the MS/HeS star is lost in the process and that it becomes a BH.

• •

BH+NS/BH - The outcome is a BH with a mass equal to the total mass of the binary before the merger.

We neglect other types of mergers as we found that they constitute only a small fraction ($\lesssim 5\%$) of all mergers and predictions for their outcomes are even less certain. For post-merger MS and HeSs we calculate their further evolution assuming that they are on ZAMS, or zero-age helium main-sequence, respectively, in order to find BH predecessors and calculate compact object’s mass. We note that presented approach is only quantitative, but helps to estimate the importance of stellar merger for BH populations.

2.3. General simulation properties

In this work, we concentrate on a description of BH populations originating from massive binaries, which we define as binaries with a primary’s (i.e., heavier star on ZAMS) initial mass $M_{\rm ZAMS,a}>10\,M_{\odot}$. Our analysis does not include BHs in triples and higher-order systems (e.g. Antonini et al., 2017), as well as BHs in binaries formed due to stellar encounters in dense stellar systems (e.g. Banerjee, 2017).

The total simulated stellar mass is $4.8\times 10^{8}\,M_{\odot}$ for the STD model and $1.1\times 10^{9}\,M_{\odot}$ ($2.4\times 10^{8}\,M_{\odot}$) for the steep IMF(flat IMF) model. Throughout this paper, we assume $50\%$ binary fraction on ZAMS for low-mass stars ($M_{\rm ZAMS}<10\,M_{\odot}$) and $100\%$ binary fraction for heavier stars ($M_{\rm ZAMS}\geq 10\,M_{\odot}$). The simulated stellar mass corresponds to about $\sim 0.8\%$ ($1.9\%$ for steep IMF model, or $0.4\%$ for flat IMF model) of the stellar mass of the MW equivalent galaxy (MWEG; $M_{\rm MW}\approx 6\times 10^{10}\,M_{\odot}$; e.g. Licquia & Newman, 2015).

For the purpose of presentation, all results are scaled to the stellar mass of a MWEG ($M_{\rm MW}\approx 6\times 10^{10}\,M_{\odot}$; Licquia & Newman, 2015) and the SFH is chosen to be constant through the last $10\,\mathrm{Gyr}$ (${\rm SFR}=6\,M_{\odot}\,yr^{-1}$). Such a model, although simple, allows us to draw general conclusions. Raw results may be scaled by using any realistic total stellar mass and any SFH may be applied. In general, the changes in SFH are not as important as the total mass of formed stars, as there is only a small delay time between the formation of a binary (ZAMS) and formation of a BH (typically less than $a\;few\times 10\,\mathrm{Myr}$), which are, typically, single. We note, that for a part of the BH population the SFH may, actually, be important, e.g. for DCOs, which have steep time-to-merger distribution ($t_{\rm merge}\propto t^{-1}$; e.g. Dominik et al., 2012). Additionally, if the spacial distribution of BHs is concerned, the change of SFH may play a role, because it affects the time which BHs have to spread throughout the space. What is more, BHs formed earlier have more time to be ejected from the stellar systems due to their NKs, or as a result of dynamical interactions. Deeper analysis of changes in spacial distributions, SFH, and chemical evolution are left for separate studies (e.g. Olejak et al. in prep.).

3. Results

The main results of the simulations are summarized in Tab. 2. Most of the BHs are predicted to exist as SBHs either as a result of a binary disruption (dSBH), or a stellar merger (mSBH; with a BH involved, or producing a star massive enough to form a BH). The population of BHs in binaries (both in DCOs and in BHBs) is about an order of magnitude smaller and in the case of models with increased BH NKs (NKR and NKBE) even two orders of magnitude smaller. Only a very small fraction of the BHBs are the interacting ones (see Sec. 3.2.3).

The initial separation influences significantly the fate of a massive binary. In general picture, if the separation is very large ($a\gtrsim 3,000\,R_{\odot}$) the binary will mostly disrupt and produce SBHs (Sec. 3.1). A low separation ($a_{\rm ZAMS}\lesssim 30\,R_{\odot}$) will lead frequently to a non-DCO merger (Sec. 3.3), which may be massive enough to form a SBH. Objects with medium initial separations have much more uncertain fate and may merge, become disrupted, or form a bound binary with at least one BH inside (Sec. 3.2). In the case of a binary harbouring a BH, the system after some time may still merge (e.g. as a double compact object merger, Sec. 3.2.1), or be disrupted during the formation of a second compact object (e.g., $\mathcal{R}_{\rm dSBH,2}$ formation route, Sec. 3.1.1).

3.1. Single BHs from disrupted binaries

If a system hosting a BH progenitor becomes disrupted, what may occur before, after, or during the BH formation, a compact object becomes an dSBH (we use name dSBH to represent a SBH from a disrupted binary, as SBH can also form due to mergers; see Sec. 3.3). Binaries may become disrupted due to a NK which one component obtains after a SN explosion, due to a Blaauw kick (e.g. because of a significant loss of mass ($\gtrsim 50\%\;M_{\rm tot}$) from a binary in SN explosion; Blaauw, 1961), or interaction with a third star. Although in our simulations we included only isolated binaries, we note that interactions, even in sparse stellar systems like the Galactic disk, may significantly alter the binary evolution (Kaib & Raymond, 2014; Klencki et al., 2017).

We found out that the tested initial parameters distributions have little effect on the properties of the resulting population of dSBHs. The steepness of the IMF changes number of resulting dSBH by $\Delta N_{\rm dSBH}\approx 70\%/160\%$ due to the decreased/increased relative number of BH progenitors for steep IMF/flat IMF models. The effect of metallicity and NK models is not as strong ($\Delta N_{\rm dSBH}\approx 40\%$ and $20\%$, respectively).

3.1.1 Typical evolutionary routes

We found that two evolutionary routes (see Tab. 3) dominate most of the dSBH production and are common for all models. Number of dSBHs produced through these routes and for different models is presented in Tab. 4. Both channels involve a binary initially composed of two massive stars that will undergo SN explosions, or collapse directly into BHs. The former case may lead to a binary disruption (due to NK and/or mass loss). The main difference between the routes comes from the initial separation. In the case of $\mathcal{R}_{\rm dSBH,1}$, the initial distance between stars is lower than $\sim 1000\,R_{\odot}$, whereas for $\mathcal{R}_{\rm dSBH,2}$ it is usually larger then $\sim 5000\,R_{\odot}$. Consequently, the former route leads through an interaction phase (typically the MT), whereas a binary evolving through the latter one experiences no interactions.

For a typical binary evolving through $\mathcal{R}_{\rm dSBH,1}$ route (Fig. 2, upper plot), the primary is about $26\,M_{\odot}$ on ZAMS, whereas secondary is slightly lighter ($\sim 21\,M_{\odot}$). The separation is modest ($\sim 300\,R_{\odot}$). The primary evolves faster and after $\sim 6.7\,\mathrm{Myr}$, while expanding as a HG star, fills its Roche lobe and a MT begins onto the secondary, which is still on its MS. Although the MT is relatively short ($\sim 4\,\mathrm{kyr}$), the primary loses its hydrogen envelope and becomes a $\sim 7\,M_{\odot}$ helium star. Half of the expelled envelope is accreted by the secondary which grows to $\sim 29\,M_{\odot}$ still being on its MS. After about $1\,\mathrm{Myr}$ primary goes through SNIb/c and becomes a $\sim 1.5\,M_{\odot}$ NS. A NK disrupts the system. Due to its significant mass, the velocity which the secondary obtains is very low ($\sim 5\,\mathrm{km}\,\mathrm{s}^{-1}$). The massive secondary needs only $\sim 2\,\mathrm{Myr}$ to become a BH with a mass of $\sim 8\,M_{\odot}$. The NK in a case of such a BH is negligible, therefore, its velocity doesn’t change significantly. It is noteworthy that most of the BHs formed through $\mathcal{R}_{\rm dSBH,1}$ route originate from secondary stars, i.e. less massive ones on ZAMS.

The $\mathcal{R}_{\rm dSBH,1}$ channel is similar to the one proposed for the formation of single massive stars (potential BH progenitors) by Renzo et al. (2019). Using different population synthesis code, they found that $86^{+11}_{-9}\%$ (depending on the model parameters) of binaries evolving through this channel will become disrupted during the first SN.

In the other typical route ($\mathcal{R}_{\rm dSBH,2}$; Fig. 2, bottom plot) the primary is about $25\,M_{\odot}$ on ZAMS, whereas the secondary is $\sim 10\,M_{\odot}$. The initial separation is as large as $\sim 6000\,R_{\odot}$. As a consequence, their Roche lobes are huge and no interaction is possible throughout their evolution. The heavier star evolves faster and in $\sim 8\,\mathrm{Myr}$ forms a $\sim 7.5\,M_{\odot}$ BH with a small NK. The separation grows to $\sim 10,000\,R_{\odot}$ what results from a loss of mass in stellar wind from both stars. After additional $\sim 17\,\mathrm{Myr}$ secondary explodes as a SN and forms a NS. A strong NK disrupts the system. Due to a large mass ratio between compact objects ($q\approx 7$) the post-disruption velocity of the BH is relatively small ($\sim 5\,\mathrm{km}\,\mathrm{s}^{-1}$).

3.1.2 Mass distribution of dSBHs

The dSBHs’ mass distribution is presented in Fig. 3. The shape is similar to the distribution of BH masses in single star evolution (compare e.g. $M_{\mathrm{remnant}}(M_{\mathrm{ZAMS}})$ relations in Belczynski et al., 2010a). The main peak ($\sim 7$–$8\,M_{\odot}$ for solar metallicity (STD model; $Z=\,Z_{\odot}$), or $\sim 15$–$30\,M_{\odot}$ for lower metallicities (models mid-$Z$ and low-$Z$; $Z=10\%\,Z_{\odot}$ and $Z=1\%\,Z_{\odot}$, respectively)) relates to BHs formed from $\sim 20$–$35\,M_{\odot}$ progenitors through failed SN explosion (e.g. Fryer et al., 2012; Belczynski et al., 2012). The second peak ($\sim 15\,M_{\odot}$ for STD, or $\sim 40\,M_{\odot}$ for mid-$Z$ model) comes from BHs originating from most massive stars ($\gtrsim 100\,M_{\odot}$) on ZAMS, which lose a large part of their mass in stellar wind. In the case of low metallicity environments (model low-$Z$) the peak results from the pair-pulsation SN (e.g. Belczynski et al., 2016c), which prevents the formation of heavier BHs in single star evolution. The third peak ($\sim 22\,M_{\odot}$ for STD, or the tail extending to $\sim 45\,M_{\odot}$ for mid-$Z$ and low-$Z$) is formed through binary interactions and is not present in distributions for single star evolution (compare Belczynski et al., 2010a, 2016c).

The similarity of the mass distributions of BHs for single stars and binaries is a direct consequence of the fact that massive binaries, which later on became disrupted, usually interact only little (mass of BH progenitor is not changed significantly; $\mathcal{R}_{\rm dSBH,1}$), or not at all ($\mathcal{R}_{\rm dSBH,2}$), so in most cases the binary components evolve as in isolation.

A notable difference is the presence of the high-mass peak and tail in the mass distributions. It is most prominent for the models with solar metallicity where the tail extends up to $\sim 22\,M_{\odot}$, whereas the heaviest BHs forming at that metallicity from single stars are only $\sim 15\,M_{\odot}$ (Belczynski et al., 2010a). These massive SBHs constitute only a tiny fraction of the population ($<1\%$). Typically, the initial total mass of a system in which such over-massive BHs are formed is above $200\,M_{\odot}$ and a mass ratio is $\sim 1$. The primary, being slightly heavier, is the first to expand and fill its Roche lobe while being on a HG. At the moment of interaction both stars have masses of about $60\,M_{\odot}$ due to strong mass losses in stellar wind. When the mass transfer ceases, the primary becomes a $\sim 25\,M_{\odot}$ core helium burning (CHeB) star. The secondary being now much heavier ($\sim 80\,M_{\odot}$), appears rejuvenated (e.g. Tout et al., 1997), evolves much faster and after $\sim 400\,\mathrm{kyr}$ forms a BH first with a mass of $\sim 21\,M_{\odot}$. The primary needs additional $100\,\mathrm{kyr}$ to become a $7.3\,M_{\odot}$ BH. During the second SN, NK easily disrupts a wide ($\sim 15,000\,R_{\odot}$) binary.

3.1.3 Significance of metallicity

The relation between the number of dSBH and metallcity was found to be non-monotonic. The largest number of dSBHs in respect to metallicity is produced in the STD model ($5.4\times 10^{7}$), whereas the smallest is in the mid-$Z$ model ($3.4\times 10^{7}$). For low-$Z$ model the number of dSBHs is slightly higher ($3.7\times 10^{7}$) than in the mid-$Z$ model.

The number of dSBHs produced through the $\mathcal{R}_{\rm dSBH,1}$ route increases monotonically with metallicity. The higher is the metallicity, the stronger is the line-driven stellar wind. Only a fraction of this mass may be accreted by a companion, thus bulk of it leaves the system what results in orbital expansion. Consequently, the stronger is the stellar wind, the wider the binary becomes. Simultaneously, the stronger mass loss also makes the companion lighter during the NS formation. Wider orbits and lighter companions make the system easier to disrupt, what as a result, produces more SBHs.

A different mechanism results in the non-monotonic relation between the metallicity and the number of BHs produced through the $\mathcal{R}_{\rm dSBH,2}$ route. Metallicity affects the evolutionary expansion of the stars which happens to be the strongest in mid-$Z$ environment and the smallest in the low-$Z$ environment (Fig. 4). Only binaries with separations high enough to avoid Roche lobe filling may evolve without interactions. The number of such systems will be smaller if the nuclear expansion is on average higher as in the case of mid-$Z$ model.

3.1.4 Initial BH velocities and high velocity dSBHs

When a SBH forms (after disruption of the binary, or BH formation, whatever happens later) its velocity is not only a result of the motion in the gravitational potential (e.g. of a galaxy), but also preceding evolutionary processes, like NK, Blaauw kick, and binary disruption. The velocity which is the result of the latter processes calculated at the moment of the binary disruption, or the formation of the BH (whatever happens later) we call here an initial BH’s velocity ($v_{\mathrm{BH,0}}$). Such a definition in which we do not involve the motion in the gravitation potential allows for application of our results to any gravitational potential (e.g. different galaxies, or different models of the MW galaxy). In this work, we provide the magnitudes of 3D velocities (i.e. lengths of the velocity vectors) and assume that the distribution of NKs is isotropic. We note that in a realistic situation, the velocity of a dSBH will be modified by the gravitational potential of a stellar system and interactions with other stars. Although these effects may also affect the binary evolution, the short dSBH formation timescale ($\lesssim 50\,\mathrm{Myr}$) allow us to assume that the effect is usually negligible.

The highest values of $v_{\mathrm{BH,0}}$ are obtained for BHs with the highest NKs. We found that majority of dSBHs with the highest initial velocities ($v_{\mathrm{BH,0}}>300$ km/s) within the highest $10\%$ were formed through $\mathcal{R}_{\rm dSBH,1}$ route and involve lowest BH masses and the shortest pre-disruption periods. We note, that Renzo et al. (2019) found that the majority of the secondaries ($\sim 90\%$) after the disruption will have low velocities ($\lesssim 20$ km/s) what agrees with our predictions for standard NK model.

$v_{\mathrm{BH,0}}$ distributions for all tested NK models are presented in Fig. 5. The maximal $v_{\mathrm{BH,0}}$ of a dSBH may reach $\sim 500\,\mathrm{km}\,\mathrm{s}^{-1}$ for STD model and nearly $\sim 1000\,\mathrm{km}\,\mathrm{s}^{-1}$ for low-$Z$ model. The higher velocities for lower metallicities are the result of a formation of closer binaries through $\mathcal{R}_{\rm dSBH,1}$ route, what together with higher masses of BHs, give larger orbital speeds (see also Renzo et al., 2019).

A distribution of $v_{\mathrm{BH,0}}$ for standard NK prescription (e.g. STD model; Hobbs et al., 2005) has two distinctive parts (Fig. 5). The low-velocity peak is connected with BHs formed with no NK due to a heavy fallback, which, as we assume, decreases the NK, or in a direct collapse associated with no NK. The wide high-velocity component comprises BHs which are much lighter, thus formed with no significant fallback. Models NKR and NKBE lead to a comparable number of dSBHs ($6.4\times 10^{7}$ and $6.5\times 10^{7}$) as the STD model ($5.4\times 10^{7}$), however, the shape of $v_{\mathrm{BH,0}}$ distribution is significantly different. NKR model lacks a low-velocity peak as the NK is independent of the fallback. Therefore, the average NK is much higher than in the STD model and typically equals $\sim 125$ km/s for a typical BH mass of $7.5\,M_{\odot}$. For STD model there is typically full fallback associated with the formation of BHs with masses of $\sim 7$–$8\,M_{\odot}$ and no NK is applied. Even more dissimilar is the $v_{\mathrm{BH,0}}$ distribution for the NKBE model, where NKs are proportional to the ratio of ejected mass to compact object mass ($v_{\rm NK}\sim M_{\rm ej}/M_{\rm BH}$; eq. 3) and typically, for $7.5\,M_{\odot}$ BHs, are $\sim 130$ km/s.

3.2. Bound systems

The number of BHs which are bound in binaries is between $5.7\times 10^{4}$ – $2.1\times 10^{7}$ (depending on the model, see Tab. 2) per MWEG. The fraction of BHs which reside in binaries for most of the tested models is $\lesssim 17\%$. For models NKR and NKBE, the fraction of BHs in binaries is $\lesssim 0.1\%$, therefore, the model of NKs highly influences the survival of BH binaries.

The maximum mass for a BH is $\sim 22\,M_{\odot}$ (or $\sim 45$ for lower metallicity models). The binary with such a massive BH exists for $\sim 100\,\mathrm{kyr}$ before it disrupts after the second SN. Afterwards, these BHs become the most heavy dSBHs (Sec. 3.1.2).

Most of the binaries harbouring BHs are wide, i.e. formed without any interactions between stars. BHs reside predominantly in DCOs (Sec. 3.2.1). The rest is mostly accompanied by WDs (such binaries we treat as non-DCOs, i.e. BHBs; Sec. 3.2.2). A small fraction of BHBs could be observed during their MT phase, although this fraction is strongly dependent on the adopted SFH (Sec. 3.2.3).

3.2.1 Double compact objects

In this study, by DCOs we understand only these harbouring at least one BH (i.e. BH+BH and BH+NS), totally neglecting NS+NS systems. BH+WD systems are included in non-DCO binaries (BHBs; see Sec. 3.2.2). The fact that the majority of binaries harbouring BHs are actually DCOs is a direct consequence of the adopted uniform distribution of mass ratios on ZAMS - if the primary is a BH progenitor ($M_{\rm ZAMS}\gtrsim 20\,M_{\odot}$) the most likely companions are NS and BH progenitors ($M_{\rm ZAMS}\gtrsim 8\,M_{\odot}$). NSs generally receive higher NKs than BHs and many potential BH+NS binaries get disrupted during the NS formation (e.g. Dominik et al., 2012), therefore, BH+BH systems typically dominate ($\sim 95$–$98\%$; see Tab. 5) over BH+NS binaries in most tested models. The exception are models with higher average NKs (NKR and NKBE) where most of DCOs are compact BH+NS systems. Many groups studied the formation of DCOs harbouring BHs over the years (e.g., Lipunov et al., 1997b; Nelemans et al., 2001; Belczynski et al., 2002b; Voss & Tauris, 2003; Dominik et al., 2012; Belczynski et al., 2016a, c; Eldridge & Stanway, 2016; Stevenson et al., 2017). However, the main focus of those studies was lied on the merger rates and properties of gravitational wave sources and, hence, they performed no deeper analysis of wider DCOs. Here, we present the analysis of the entire DCO population as a part of the total BH population. Therefore, our results may be applied not only to the study of merging DCOs, but also to other phenomena, e.g. microlensing by DCOs.

In Tab. 5 we present detailed results for the main models. The number of BHs in DCOs is typically between $1.5\times 10^{6}$ – $1.9\times 10^{7}$ (note that BH+BH counts as two BHs), although for NKR and NKBE models, the number drops to $5.1\times 10^{3}$ and $2.5\times 10^{2}$, respectively. The low number of DCOs in the models with higher average NKs is connected to the fact that most of binaries with massive components which are potential progenitors of DCOs are wide (route $\mathcal{R}_{\rm DCO,1}$; see below), thus easily disrupted by even moderate NKs.

There are two main evolutionary routes leading to the formation of DCOs harbouring BHs (Tab. 3). The main feature that distinguish these routes is the initial separation. It is connected with the lack ($\mathcal{R}_{\rm DCO,1}$), or presence ($\mathcal{R}_{\rm DCO,2}$) of the pre-DCO formation interactions (MT and CE), which results in the formation of wide ($a\gtrsim 10^{4}$), or close ($a\lesssim 10^{3}\,R_{\odot}$) DCOs, respectively (see Fig. 6). Among close DCOs (mainly route $\mathcal{R}_{\rm DCO,2}$) are potential DCO merger progenitors (e.g. Abbott et al., 2016a) and their number and relative fraction strongly depends on metallicity and NKs as previously noted by Chruslinska et al. (2019).

Route $\mathcal{R}_{\rm DCO,1}$ is the main formation route of DCOs if standard NK model is concerned. In this formation channel none of strong interactions are present (Fig. 7). In a typical situations the stars on ZAMS are not very massive ($\sim 25$ and $\sim 22\,M_{\odot}$) and the separation is large ($\sim 12,300\,R_{\odot}$). After about $7.7\,\mathrm{Myr}$ the primary explodes and forms a $\sim 7.5\,M_{\odot}$ BH. In additional $\sim 1\,\mathrm{Myr}$, the secondary forms a $\sim 7.5\,\mathrm{Myr}$ BH in a SN explosion. The final separation is $\sim 37,500\,R_{\odot}$.

The other formation route ($\mathcal{R}_{\rm DCO,2}$) is responsible for much smaller fraction of DCOs ($\lesssim 21\%$), but systems which form through this route have much shorter separations and, thus, may merger after less than $10\,\mathrm{Gyr}$ due to the emission of gravitational waves. In a typical case, the binary is massive on ZAMS ($M_{\rm ZAMS,a}\approx 57\,M_{\odot}$ and $M_{\rm ZAMS,b}\approx 38\,M_{\odot}$) and the separation is small ($a_{\rm ZAMS}\approx 600\,R_{\odot}$). The primary evolves fast and after $\sim 4.1\,\mathrm{Myr}$ fills its Roche lobe while expanding on the HG. The MT commences and results in the loss of nearly entire hydrogen envelope by the primary, which shortly after becomes a $\sim 16\,M_{\odot}$ HeS. Half of the primaries envelope is accreted by the secondary which grows to $\sim 49\,M_{\odot}$ while still on its MS. After additional $\sim 600\,\mathrm{kyr}$, the primary forms a $\sim 8\,M_{\odot}$ BH in a SN explosion. The secondary starts to expand after leaving the MS and while being a CHeB star fills its Roche lobe and a CE commences. As a result the separation shrinks from $\sim 4000\,R_{\odot}$ to about $5.2\,R_{\odot}$ and the secondary loses its entire envelope and becomes a $\sim 15\,M_{\odot}$ HeS. A second SN occurs $\sim 600\,\mathrm{kyr}$ later and $\sim 7.7\,M_{\odot}$ BH forms. The NK is not very strong, thus the separation remains small ($\sim 7.3\,R_{\odot}$). Such a binary is estimated to merge after $\sim 420\,\mathrm{Myr}$. This channel was previously thoroughly analysed in the context of DCO mergers by e.g. Belczynski et al. (2002b); Dominik et al. (2012); Belczynski et al. (2016a); Woosley (2016); Kruckow et al. (2018).

BH+NS formation through $\mathcal{R}_{\rm DCO,1}$ and $\mathcal{R}_{\rm DCO,2}$ routes is qualitatively similar to the formation of BH+BH systems, although the ZAMS secondary masses are significantly lower ($\sim 3\times$ for $\mathcal{R}_{\rm DCO,1}$ and $\sim 1.5\times$ for $\mathcal{R}_{\rm DCO,2}$). Additionally, the progenitors are more easily disrupted during the second SN and have a smaller chance of CE survival due to more extreme mass ratio. This results in a significantly lower fraction of BH+NS systems among DCOs than BH+BH systems, in spite of a higher abundance of potential progenitors in initial (ZAMS) populations111The ratio of potential BH+NS progenitors to potential BH+BH progenitors (i.e. binaries with components’ masses high enough to form a NS, or BH in single star evolution) may be calculated as

$f=\frac{\int_{M_{\rm NS,ZAMS,max}}^{150\,M_{\odot}}P_{\rm BH+NS}(M_{1})\,{\rm IMF}(M_{1})\,dM_{1}}{\int_{M_{\rm NS,ZAMS,max}}^{150\,M_{\odot}}P_{\rm BH+BH}(M_{1})\,{\rm IMF}(M_{1})\,dM_{1}},$

where $P_{\rm BH+NS/BH+BH}(M_{1})$ is the probability that a companion of a BH progenitor with mass $M_{1}$ is a NS/BH progenitor. Assuming flat mass ratio distribution for companions ($P(q)=1/(M_{1}-M_{\rm NS,ZAMS,min})$),

$P_{\rm BH+NS}(M_{1})=\frac{M_{\rm NS,ZAMS,max}-M_{\rm NS,ZAMS,min}}{M_{1}-M_{\rm NS,ZAMS,min}}$

and

$P_{\rm BH+BH}(M_{1})=\frac{M_{1}-M_{\rm NS,ZAMS,max}}{M_{1}-M_{\rm NS,ZAMS,min}},$

where $M_{\rm NS,ZAMS,min/max}=8/22\,M_{\odot}$ are the minimal/maximal ZAMS masses for a NS progenitor (assuming no interactions with the companion) for solar metallicity ($\,Z_{\odot}$). ${\rm IMF}(M_{1})\propto M_{1}^{\Gamma}$ is the initial mass function. $150\,M_{\odot}$ is the upper limit for the ZAMS mass adopted in our calculations. Assuming $\Gamma=-2.3$ and solar metallicity (STD model) $f\approx 1.1$ meaning more BH+NS than BH+BH progenitors in initial populations. This value is slightly affected by metallicity, as $Z=1\%\,Z_{\odot}$ (low-$Z$ model; where $M_{NS,ZAMS,max}\approx 19\,M_{\odot}$) gives $f\approx 1$, and significantly affect by the steepness of the IMF, as $\Gamma=-1.9$ (flat IMF model) gives $f\approx 0.9$, i.e. more BH+BH progenitors.. There are typically (except NKR and NKBE models) $\gtrsim 20$ times more BHs in BH+BH binaries than in BH+NS. We observe no significant difference between the shapes of BH+BH and BH+NS distributions of orbital separations at the moment of of DCO formation (Fig. 6).

We have found that the number of BHs in DCOs is $3$–$4$ times higher in lower metallicity environments than in the STD model. The number of wide DCOs ($\mathcal{R}_{\rm DCO,1}$) is only slightly affected by metallicity (by a factor of $\lesssim 2.5$), whereas the number of close DCOs ($\mathcal{R}_{\rm DCO,2}$) is significantly affected (by more than two orders of magnitude). The latter is an effect of higher BH masses (typically $\sim 15$) in low metallicity environments in comparison to $\,Z_{\odot}$ environments (typically $\sim 7$–$8\,M_{\odot}$). Heavier BHs in the standard NK model usually obtain lower NKs due to fallback, or no NK when formed in a direct collapse. Therefore, DCOs in lower metallicity environments, where BHs form with higher masses, are less frequently disrupted (e.g. Belczynski et al., 2010b).

Significantly different situation concerns the relation for close systems ($\mathcal{R}_{\rm DCO,2}$), which in lower metallicity environments are about two orders of magnitude more numerous than in the STD model. This results from an interplay between mass loss in stellar wind and stellar expansion (Fig. 4). The former is proportional to metallicity and leads to orbital widening due to mass loss from the system. On the other hand, solar metallicity stars do not expand as significantly as these in mid-$Z$ model. This means that massive binaries in high-metallicity environment (e.g. STD model) more frequently have separations too large to go through a CE phase which is an essential step in the $\mathcal{R}_{\rm DCO,2}$ route. Consequently, lower-metalicity models have more close DCOs, than high metallicity ones. A similar result was recently obtained by Spera et al. (2019) with the use of SEVN population synthesis code.

Distributions of component masses are shown in Fig. 8. For $\mathcal{R}_{\rm DCO,1}$ formation channel (wide DCOs) distributions for primaries and secondaries are joined together, except BH+NS systems for which only BH masses are shown. The results resemble these for BHs from disrupted binaries (Fig. 3). For $\mathcal{R}_{\rm DCO,2}$ channel (close DCOs) the total mass of a binary is provided. Such information is more important from the point of view of microlensing surveys, because a close binary will most probably act as a single lens. For lower metallicites the total mass may reach $\sim 70$–$80\,M_{\odot}$ what matches the current range of observed mass in BH+BH mergers (The LIGO Scientific Collaboration et al., 2018a, see also Belczynski et al. 2016a).

Within our models up to $\sim 15\%$ of BHs are found in DCOs (Tab. 2) which are predominantly wide (Voss & Tauris, 2003). For the standard NK model, up to $19\%$ of them have merger times ($t_{\rm merge}$) smaller than $10\,\mathrm{Gyr}$, what mainly depends on metallicity (see Tab. 5; merging BH+BH systems form much more efficiently at low Z, e.g. Dominik et al., 2013; Giacobbo & Mapelli, 2018; Klencki et al., 2018; Chruslinska et al., 2019). In the NKR and NKBE models, although the fraction of merging systems may be much higher ($\sim 37\%$ for NKR model), the number of these systems is significantly smaller. It is a consequence of higher average NKs ($\sim 125$–$130$ km/s in contrast to $\sim 10$ km/s in STD model), what frequently disrupts wide binaries.

Typical center-of-mass initial velocities ($v_{\rm COM,0}$; i.e. binary velocities just after the second SN not including velocity resulting from the motion in a gravitational potential) are below $\sim 20$ km/s (Fig. 9). Only less than $15\%$ of DCOs (mainly compact; route $\mathcal{R}_{\rm DCO,2}$) possess velocities exceeding $20$ km/s. Therefore, vast majority of the DCOs remains in the vicinity of their birth places. Highest velocities ($v_{\rm COM,0}\gtrsim 200$ km/s) are obtain by BH+NS systems, whereas systems with lower velocities are dominated by BH+BH binaries, what results from typically higher NKs obtained by NSs, than BHs. Typical velocities of merging systems ($t_{\rm merge}<10\,\mathrm{Gyr}$) are below $\sim 50$ km/s for STD model, so most of the mergers are predicted to occur in the vicinity of the formation places (e.g. Perna et al., 2018). However, some systems may obtain also high velocities ($\gtrsim 200$ km/s). These are mainly BH+NS systems, which we predict to merge away from their birth environments. For NKR and NKBE models, DCOs obtain typically intermediate velocities ($\sim 50$–$130$ km/s; mainly route $\mathcal{R}_{\rm DCO,2}$; typically, $t_{\rm merge}\lesssim 10\,\mathrm{Gyr}$), although some systems (mainly BH+NS) may have velocities as high as $\sim 400$ km/s.

3.2.2 Wide non-DCO binaries

Throughout the paper we use the term wide BHB referring to systems harbouring a BH with a non-compact companion and experiencing no strong interaction during their entire evolution (the separation between the stars is large and the stars never fill their Roche lobes). In this section, we focus on wide BHBs, excluding wide DCOs, which were discussed in previous section.

Tab. 6 summarizes our results for the main models. The typical number of wide BHBs is between $5.5\times 10^{5}$ and $3.3\times 10^{6}\,M_{\odot}$. Only for models with higher average NKs (NKR and NKBE) the predicted number of such binaries is much lower ($3.1\times 10^{3}$ and $1.0\times 10^{-4}$, respectively). The probability for a binary to obtain a low NK ($\lesssim 20$ km/s) in NKR and NKBE models is much lower than in STD model (c.f. Fig. 5) and typically NKs are $\gtrsim 100$ km/s. Therefore, in these models wide binaries are more easily disrupted during the formation of a BH. In contrast, standard NK model (Hobbs et al., 2005) predicts mostly low, or negligible NKs for BHs which allows for more frequent survival of wide binaries.

The number of wide BHBs is higher in STD model than in mid-$Z$ and low-$Z$ models mainly because of the more significant orbital expansion resulting from mass loss in stellar wind, which is stronger for higher metallicities. Additionally, stellar expansion is higher in mid-$Z$ model, thus futher increasing the chance of interactions (Fig. 4). In low-$Z$ model the stellar expansion is smaller than in the mid-$Z$ model, therefore, results in higher number of wide BHBs.

BH companions in wide systems are predominantly CO WDs ($74$–$89\%$). WDs, being the final evolutionary stage of low-mass ($\lesssim 8\,M_{\odot}$) stars’ evolution, are simultaneously the longest evolutionary stage of the evolution of stars with initial masses $2\lesssim M\lesssim 8\,M_{\odot}$ during the Hubble time. For comparison, low mass stars ($M\lesssim 1\,M_{\odot}$) spend most of this time on the MS, whereas, heavier stars ($M\gtrsim 8\,M_{\odot}$) quickly ($\lesssim 50\,\mathrm{Myr}$) end their evolution and form a NS, or a BH. The typical separations of wide BHBs are $a\gtrsim 5000$.

In an exemplary evolution of a typical wide BHB, the progenitor system on ZAMS is composed of a $\sim 25\,M_{\odot}$ primary and a $6.0\,M_{\odot}$ companion (see Fig. 10). The initial separation is large ($\sim 16,000\,R_{\odot}$), which precludes any interactions. The separation further increases due to wind mass-loss from the stars and just before the SN reaches $\sim 34,000\,R_{\odot}$. The primary, being heavier, evolves faster and in $\sim 7.8\,\mathrm{Myr}$ becomes a BH. The NK and the Blaauw kick are crucial at that point. Only a small total kick will allow for a survival of the binary. The secondary is still on its MS for another $\sim 60\,\mathrm{Myr}$. At the age of about $76\,\mathrm{Myr}$ a CO WD forms. The final separation may reach $56,000\,R_{\odot}$.

Mass distribution of BHs in wide BHBs (Fig. 11) resembles this for dSBHs as most dSBHs are not interacting significantly prior to disruption, thus frequently originate from wide BHBs. Two distributions differ in the presence (dSBHs), or the lack (wide BHBs) of the high mass peak ($\sim 22\,M_{\odot}$) for STD model. Also the high-mass tail (above $\gtrsim 15\,M_{\odot}$ for STD model and $\gtrsim 40\,M_{\odot}$ for mid-$Z$ and low-$Z$ models), which originates from wind accretion by a BH, or its progenitor, in the case of non-interacting binaries is less pronounced, than in distributions for interacting binaries. The maximal BH masses are $\sim 16\,M_{\odot}$ for STD model and $\sim 46\,M_{\odot}$ in lower metallicities. The distribution of companion masses is dominated by WDs (carbon-oxygen and oxygen-neon), which form a characteristic peak between $\sim 1\,M_{\odot}$ and the Chandrasekhar mass ($M_{\rm Ch}=1.44\,M_{\odot}$). The long tail on this distribution is mainly composed of MS stars, which may have masses up to $\sim 80\,M_{\odot}$ for solar metallicity, or $\gtrsim 100\,M_{\odot}$ for lower metallicity models. The massive companions ($M\gtrsim 8\,M_{\odot}$) are going to explode in SN explosion which may disrupt the BHBs (either through a NK, or a Blaauw kick) populating dSBHs, or form a wide DCOs (discussed earlier in Sec. 3.2.1).

Although in rare cases the initial center-of-mass velocity ($v_{\rm COM,0}$) of wide BHBs may reach $\sim 80$ km/s, typically ($\gtrsim 95\%$ of cases) they are negligible ($\lesssim 10$ km/s). Due to the fact that orbital velocities are typically also small ($\lesssim 20$ km/s), the binaries rarely leave their birth places if evolved in isolation. We note, that in realistic situation a wide BHB is frequently influenced by fly-by encounters (e.g. Klencki et al., 2017).

3.2.3 Mass Transferring binaries with BH accretors

Mass transferring BHBs (MTBHB) have different evolutionary routes and properties than the typical BHBs, which are dominated by non-interacting systems (wide BHBs; Sec. 3.2.2). Although we note that BHs were detected also in high-mass XRBs where MT rates high enough to fuel an accretion disk may occur also through companion’s stellar wind (e.g. Ruhlen et al., 2011), in this study only RLOF systems were included.

For any adopted SFH the fraction of BHs accreting mass from their companions through RLOF is small. However, this small fraction is more a result of the brevity of MT phase in comparison to the evolutionary time-scales, than the rarity of BHBs that are close enough to commence RLOF. Here, for the sake of presentation, we assumed a constant star formation for the last $10\,\mathrm{Gyr}$. In such a case, only $<1\%$ (see Tab. 6) of all BHBs contain a donor which is filling its Roche lobe and transferring mass onto a BH. MTBHBs can be perceptible as an X-ray binaries (e.g. Tauris & van den Heuvel, 2006). We found out that in our models $2$ – $28\%$ of BHs have gone through a MT phase which lasts typically $\lesssim 40\,\mathrm{Myr}$. A BH accretes on average $0.7$ – $3.1\,M_{\odot}$ in that time.

Our simple model predicts $32$ – $4,200$ MTBHBs per MWEG. The number is the highest in lower metallicities ($4.4\times 10^{2}$ and $4.2\times 10^{3}$ systems for low-$Z$ and mid-$Z$ models, respectively) where it is easier to start RLOF MT (Linden et al., 2010). Also in SS0 model we see $\sim 8$ times higher predicted number of MTBHBs, than in the STD model. It is a result of assumed in this model thermal distribution of eccentricities ($P(e)\sim e$), what significantly increases average initial eccentricity, thus lowering the initial periastron distance, what makes RLOF easier.

There are two main evolutionary routes (Tab. 3) typical for all models (comprising $79\div 99\%$ of all MTBHBs). Both channels differ mainly in the mass of an accreting BH. In the $\mathcal{R}_{\rm MTBHB,1}$ route, typical BH masses are above $\sim 5\,M_{\odot}$, whereas in $\mathcal{R}_{\rm MTBHB,2}$ route, they are around $\sim 3\,M_{\odot}$. This distinctness stems from the different BH formation mechanisms. In $\mathcal{R}_{\rm MTBHB,1}$, a BH forms directly after a SN, whereas in $\mathcal{R}_{\rm MTBHB,2}$, a NS forms first and after a period of mass accretion reaches a critical mass ($M_{\rm max,NS}$) and collapses to a BH. The critical mass is in general not well constrained and depends on the applied equation of state and rotation (e.g Kalogera & Baym, 1996). In our simulations a value of $M_{\rm max,NS}=2.5\,M_{\odot}$ was adopted.

An exemplary system which becomes a MTBHB through $\mathcal{R}_{\rm MTBHB,1}$ route begins its evolution having ZAMS masses of $42\,M_{\odot}$ and $5.4\,M_{\odot}$ (Fig. 12, upper plot). The separation is moderate, $a_{\rm ZAMS}\approx 5000\,R_{\odot}$. In $\sim 4.9\,\mathrm{Myr}$, the primary evolves off the MS and commences a CE phase. Afterwards, the separation is reduced to $\sim 13\,R_{\odot}$. After additional $0.5\,\mathrm{Myr}$, the primary becomes a BH receiving a small NK. The companion needs additional $\sim 60\,\mathrm{Myr}$ to expand due to nuclear evolution and fill its Roche lobe. The MT prolongs for $16\,\mathrm{Myr}$ during which the companion loses $\sim 90\%$ of its mass and evolves off the MS. The BH accretes about $1.5\,M_{\odot}$.

As far as the $\mathcal{R}_{\rm MTBHB,2}$ route is concerned (Fig. 12, lower plot), the initial masses on ZAMS are much smaller and in a typical case equal $\sim 23\,M_{\odot}$ and $\sim 3.3\,M_{\odot}$. The initial separation is large ($a\approx 9400\,R_{\odot}$). Similarly to $\mathcal{R}_{\rm MTBHB,1}$ the primary evolves off the MS and commences a CE phase while being a $\sim 10\,M_{\odot}$ CHeB star. The outcome is a compact binary ($a\approx 40\,R_{\odot}$) composed of a $\sim 7.9\,M_{\odot}$ HeS and $\sim 3.5\,M_{\odot}$ MS star. The primary is not heavy enough to form a BH. Instead, in $100\,\mathrm{kyr}$ forms a heavy ($\sim 1.8\,M_{\odot}$) NS with a strong NK. The orbit becomes highly elongated ($a\approx 200\,R_{\odot}$, $e\approx 0.96$). The secondary, which is now the more massive star, expands due to nuclear evolution and after $\sim 200\,\mathrm{Myr}$ fills its Roche lobe during periastron passage and a MT begins. In such a situation, we assume that the orbit is immediately circularized. In $\sim 500\,\mathrm{Myr}$, the NS accretes $>0.7\,M_{\odot}$, the critical mass is reached, and it collapses to a low-mass ($2.5\,M_{\odot}$) BH. Afterwards, the MT restarts when the secondary is evolving off the MS. Finally, the BH grows to a mass of $\sim 3\,M_{\odot}$.

Fig. 13 shows distribution of masses for BHs and companions. The leftmost peak (both for BHs and companions distributions) is composed of BHs formed through $\mathcal{R}_{\rm MTBHB,2}$ route, whereas, heavier BHs and companions represent $\mathcal{R}_{\rm MTBHB,1}$ route. The maximal BH masses, naturally, exceed these for wide BHBs, or expected from single star evolution (Belczynski et al., 2010a). The heaviest companions are usually accompanied by the heaviest BHs what allows to avoid extreme mass ratios, which in our calculations are assumed to lead to dynamical instability during MT.

Models mid-$Z$ and low-$Z$ predict larger populations of MTBHBs than STD model. In lower metallicity environments it is easier to produce heavier BHs, thus companions Roche lobes are relatively smaller. Also, the nuclear expansion is slower what allows for longer phases of stable MT.

Typical initial center-of-mass velocities ($v_{\rm COM,0}$) are smaller than $20$ km/s. It is a result of the fact that both heavy BHs ($M\gtrsim 40\,M_{\odot}$) and heavy NSs ($M>1.8\,M_{\odot}$) present in MTBHBs formed through channels $\mathcal{R}_{\rm MTBHB,1}$ and $\mathcal{R}_{\rm MTBHB,2}$, respectively, and have low, or negligible NKs. The exception are models with higher average NKs (NKR and NKBE) which give much higher $v_{\rm COM,0}$ ($50$–$100$ km/s). We note, that velocities around $100$ km/s are also attainable in models with standard NK distribution (e.g. STD) in systems with lighter BHs ($M_{\rm BH}\approx 6\,M_{\odot}$).

3.3. Mergers

Mergers are a frequent outcome in population synthesis calculations of isolated binaries due to a failed CE ejection. Our results show that most typically (more than $95\%$ of cases) mergers occur when the components are on MS, HG, or are HeS (Fig. 14). Although our results are not representative for the entire binary star population, as we include only binaries with initial primary’s mass in the range $10\,\hbox{--}\,150\,M_{\odot}$, we include all potential progenitors of SBHs from merging binaries for which $M_{\rm ZAMS,a}+M_{\rm ZAMS,b}\geq 20\,M_{\odot}$. We note that even if some lower-mass binaries may produce a BH predecessor after a merger, the poorly understood merger physics does not allow to include such cases in the following analysis. Throughout the paper we use a designation mSBH to distinguish SBHs formed from merger products from these originating from disrupted binaries.

In order to include merger products in the BH population we have implemented a simplified formalism of Olejak et al. (in prep.; see Sec. 2.2). After merger, the products were evolved as single stars till the formation of a compact object, or reaching the age of $15\,\mathrm{Gyr}$. This way, we were able to roughly estimate the population of mSBHs. According to our simulations, the minimal ZAMS mass that produce a BH in single-star evolution is between $19$–$22\,M_{\odot}$ depending on the metallicity. Similarly, the minimal zero-age helium main-sequence mass that produce a BH is $\sim 9\,M_{\odot}$ (e.g. Woosley, 2019). In a typical case (MS+MS and HG+MS mergers) there are no interactions prior to merger.

The predicted numbers of BHs for all tested models are provided in Tab. 2. The model dependence is very low with $\sim 3.6\,\hbox{--}\,6.6\times 10^{7}$ BHs originating from mergers, except flat IMF and steep IMF models where the contrast is significant ($1.1\times 10^{8}$ and $1.6\times 10^{7}$, respectively) due to different (higher, or lower, respectively) number of massive stars on ZAMS.

The mass distribution (Fig. 15) shows a major peak at $\sim 7.5\,M_{\odot}$ ($\sim 15\,M_{\odot}$ for mid-$Z$ and low-$Z$ models). The distribution noticeably differs from distributions of BH masses obtained through other formation channels. Particularly interesting is the presence of BHs inside the Mass Gap ($\sim 2.5$ – $5\,M_{\odot}$) were no compact objects have been detected through observations (e.g. Bailyn et al., 1998; Özel et al., 2010; Farr et al., 2011). These BHs mainly originate from mergers of NSs with HeS thus exceeding the maximal mass for a NS ($M_{\mathrm{max,NS}}$). In a typical case, the components masses on ZAMS are $\sim 10\,M_{\odot}$ and $\sim 9.6\,M_{\odot}$. The primary, being slightly heavier, evolves faster and in $\sim 24\,\mathrm{Myr}$ fills its Roche lobe while expanding on a HG. The mass transfer leads to a mass reversal with primary becoming a $2.2\,M_{\odot}$ HeS, whereas secondary grows to $14\,M_{\odot}$ still remaining on a MS. In about $4\,\mathrm{Myr}$ primary explodes as an electron capture supernova and forms a NS with a mass of $1.26\,M_{\odot}$ and negligible NK. The separation is $770\,R_{\odot}$ at the moment. Secondary needs additional $1\,\mathrm{Myr}$ to fill its Roche lobe during a CHeB phase. The mass ratio is extreme with the donor being ten times heavier then the accretor, what leads to a CE and, in a consequence, merger of both stars. Finally, according to our prescription a $2.9\,M_{\odot}$ BH forms. Such low-mass free-floating BHs can be detectable by microlensing methods Wyrzykowski et al. (e.g. 2015). We note that the formation of low-mass BHs through this channel heavily depends on our assumptions concerning mergers (see Sec. 2.2) and, particularly, the fraction of the HeS which is actually accreted onto the NS.

Another striking feature of the distribution is the presence of high mass tails exceeded up to $\sim 50\,M_{\odot}$, $\sim 80\,M_{\odot}$, or $\sim 130\,M_{\odot}$ for STD, mid-$Z$, and low-$Z$ models, respectively (Fig. 15). Typically, in STD model, progenitors of such massive BHs are formed from a merger of a $\sim 25\,M_{\odot}$ HeS with a $46\,M_{\odot}$ MS star which leaves a $\sim 43\,M_{\odot}$ HeS as an outcome. Such a massive HeS collapses directly to a BH with $\sim 40\,M_{\odot}$ mass. Evolution for massive BHs in other metallicities are analogous. These massive stars, although rare, may produce significantly stronger microlensing events than those with typical masses ($\sim 7$–$8\,M_{\odot}$). W note that Spera et al. (2019) obtained even heavier single BHs as a result of stellar mergers using different prescription for post-merger stellar parameters. In contrast to our results, their heavies BHs (up to $\sim 145\,M_{\odot}$ for metallicity $Z=0.0001$) originate from mergers of two MS stars.

BHs originating from mergers will have typically very low initial velocities ($v_{\mathrm{BH,0}}$). According to our results more than $75\%$ of the post-merger BHs will have no additional velocity to that resulting from the movement in the galactic potential ($v_{\mathrm{BH,0}}\approx 0$). Only BHs originating from double NS mergers may obtain high $v_{\mathrm{BH,0}}$ (above $100\,\mathrm{km}\,\mathrm{s}^{-1}$), however, they constitute only $<0.1\%$ of all post-merger BHs. We note, that asymmetries in the merger process and mass ejections may be a source of significant NKs, but are not involved in this study.

4. Application: Microlensing of bulge stars by isolated BHs

Gravitational microlensing effect is an amplification of the light of a distant object by another object (luminous, or not) aligned in the line-of-sight toward the source (Paczynski, 1996). Unfortunately, lens mass measurement, which is necessary to conclude on its stellar, or BH nature, is typically prone to parameter degeneracies (e.g. Wozniak et al., 2001; Sumi et al., 2013; Wyrzykowski et al., 2015; Mróz et al., 2017) and BH lenses can be found only in particular circumstances (e.g. Bennett et al., 2002; Dong et al., 2007; Wyrzykowski et al., 2016; Rybicki et al., 2018). On the other hand, population studies allow us to estimate the expected number of microlensing events with BH lenses, what we do in the following analysis. Other similar attempts include studies by Osłowski et al. (2008) and Lu et al. (2019).

Due to the high star number density and relative proximity, the Galactic bulge is a frequent direction of observations aimed at detecting microlensing events (e.g. Gould, 2000; Wyrzykowski et al., 2015), therefore, here we also concentrate on this particular direction.

We make a series of simplifying assumptions. First of all, we assume that the source is always located in the bulge, which is represented by a flat disk of size equal $\theta_{\rm bulge}=31$ square degrees, possessing $N_{\star,\rm bulge}=1.5\times 10^{8}$ stars. The tangential velocity distribution (radial velocity can be ignored in microlensing) is on average zero and has a dispersion $\sigma_{\rm bulge,z}=\sigma_{\rm bulge,y}=80$ km/s in both directions: parallel (y-axis) and perpendicular (z-axis) to the galactic plane (Skowron et al., 2011).

Secondly, we assume that a BH lens is always located in the galactic disk, which has a total mass of $5\times 10^{10}\,M_{\odot}$ (Licquia & Newman, 2015) and its stellar number density (we assume that BHs have the same distribution as stars in general) is described by

[TABLE]

where $z$ is the distance from the galactic plane and $R$ is the distance from the galactic centre in the galactic plane, following Batista et al. (2011) and Skowron et al. (2011). We assume that components of the tangential velocity of stars in the disk are on average $v_{\rm disk,mean,z}=0$ in the direction perpendicular to the galactic plane and $v_{\rm disk,mean,y}=200$ km/s in the galactic plane. The tangential velocity dispersion is assumed to be $\sigma_{\rm disk,z}=40$ km/s perpendicular to the galactic plane and $\sigma_{\rm disk,y}=55$ km/s parallel to galactic plane.

In contrary to other studies, here we use the BHs from disrupted binaries and stellar mergers only, as BHs in binaries represent small fraction ($<20\%$) of the total population (Tab. 2). Additionally, microlensing by binary lens is more complex (caustic crossing, high amplification, etc.) and will be investigated in a separate study. Also, due to the assumed $100\%$ binary fraction for stars heavier than $10\,M_{\odot}$ on ZAMS, we do not expect any BHs formed through single-star evolution.

Stellar mergers may potentially constitute a significant part of the BH population (see Sec. 3.3), so it is necessary to include them in the calculation of the microlensing rate. In such a case, not only the number of BHs is important, but also their masses and proper motions. We calculate the final BH mass following the formulas of Olejak et al. (in prep.; see Sec. 2.2 for the post merger mass and evolutionary phase). As far as the proper motions are concerned, we assume that the $v_{\mathrm{BH,0}}$ of the post merger star is equal to the pre-merger center-of-mass velocity of the binary, i.e. merger occurs without a kick.

We use the Monte Carlo method to sample the spacial and velocity distributions of both lenses and sources. In our analysis we include only lenses which were localized between the observer, assumed to be located at $x_{\odot}=8.3$ kpc from the Galactic centre, and the Galactic bulge, ignoring less likely lenses and sources from the far Galactic disk.

An angular Einstein radius of a lens depends on its mass ($M_{\mathrm{L}}$) and distances between the observer and the lens ($D_{\mathrm{L}}$), or the source ($D_{\mathrm{S}}$):

[TABLE]

where $\kappa=\frac{4G}{c^{2}}\approx 8.1~{}\frac{\mathrm{mas}}{\,M_{\odot}}$. If the relative proper motion between the lens and the source is $\vec{\mu}_{\mathrm{rel}}=\vec{\mu}_{\mathrm{L}}-\vec{\mu}_{\mathrm{S}}$, the source crosses the Einstein radius in time called Einstein Radius crossing time, or simply, event’s time-scale, $t_{\mathrm{E}}=\frac{\theta_{\mathrm{E}}}{|\mu_{\mathrm{rel}}|}$.

The relative proper motion is computed using the velocities of the lens and the source as:

[TABLE]

and

[TABLE]

where $v_{\rm rel,y/z}$ is a component of the relative velocity along y-/z-axis. $v_{\rm BH,y/z}$, $v_{\rm source,y/z}$, and $v_{\rm Earth,y/z}$, are the velocities of the BH, the source, and the observer, along y-/z-axis, respectively. The observers velocity was assumed to be $v_{\rm Earth,y/z}=230/15.5$ km/s. $x=x_{\odot}-D_{\rm L}$ is the distance of the lens from the Galactic centre.

Finally, we calculated an estimated number of microlensing events (${\rm E}(N_{\rm ML})$) during $t=1$ yr as (assuming $\theta_{\rm E}$ is small),

[TABLE]

where, $\theta_{\rm E,i}$ is the Einstein radius generated by i-th BH, $v_{\rm rel,i}=\sqrt{v_{\rm rel,x,i}^{2}+v_{\rm rel,y,i}^{2}}$ is the relative tangential velocity of the i-th BH, and $\Omega_{\rm bulge}$ is the considered area of the bulge in steradians. The summation goes over all the included BHs, i.e localized between the observer and the bulge.

If we assume that the stellar mass of the Galactic disk is $\sim 5\times 10^{10}\,M_{\odot}$ (Licquia & Newman, 2015) and we take the number of stars in the bulge as $N_{\star,bulge}=1.5\times 10^{8}$ from the number of monitored stars by the OGLE-III survey covering 31 sq.deg of the bulge (Wyrzykowski et al., 2015), we predict that there might be as many as $\sim 14$ yr*-1* microlensing events due to BHs. In the so-far analysed OGLE-IV data (Udalski et al., 2015; Mroz et al., 2019), which covers $160$ sq.deg and monitors $400$ million sources, the number of expected events due to BHs grows to $\sim 26$ yr*-1*. Here we include only events with characteristic time $t_{\mathrm{E}}$ in range $1$–$1000$ days which includes most of the events ($>95\%$). Although BHs from disrupted binaries have higher velocities and, therefore, higher chance for microlensing, BHs from disrupted binaries (dSBH; Sec. 3.1) dominate the rate ($\sim 55\%$ of events) due to higher masses and greater abundance.

For a given source and lens distances and their velocities we computed $\theta_{\mathrm{E}}$ and $\mu_{\mathrm{rel}}$ and hence the distribution of time-scales ($t_{\mathrm{E}}$, Fig. 16). The typical values of $t_{\mathrm{E}}$ are on average $\sim 97$ days for BH originating from disrupted binaries (dSBH; Sec. 3.1) and $\sim 108$ days for these originating from stellar mergers (mSBH; Sec. 3.3.

Events with such time-scales are detectable in current microlensing surveys with efficiency of more than $75\%$ (e.g. Sumi et al., 2013; Wyrzykowski et al., 2015). This means in the OGLE-III data covering 8 years of monitoring of $150$ million of stars there should be about $84$ events due to BHs. Since, as shown above, these tend to have longer time-scales, the annual Earth motion (the parallax effect) is likely to cause the distortion to the standard microlensing light curve, therefore, it is possible that some fraction of those $84$ events are not in the Wyrzykowski et al. (2015) sample of standard events. Indeed, Wyrzykowski et al. (2016) have identified $59$ OGLE-III events with significant parallax effect, among which $13$ are likely due to NSs and BHs. The improved analysis of Wyrzykowski & Mandel (2019) has increased the number of NS and BH lenses to $18$ in that sample. The remaining undetected BH lenses might still await discovery in the OGLE-III data with moderate and small annual parallax signatures.

5. Discussion

5.1. Detectability of single BHs

Although the predicted number of SBHs in MWEG (both dSBH and mSBH) is very significant ($3.2\times 10^{7}$ – $2.5\times 10^{8}$), these objects still evade detection. Already Shapiro & Teukolsky (1983) estimated that BHs in the Galaxy should be counted in millions. Timmes et al. (1996) estimated that there are even $1.4\times 10^{9}$ BHs in the Galaxy. Samland (1998) calculated the BH population to be $1.8\times 10^{8}$. They included the changes in SFR, but excluded binary stars from their analysis. It was proposed that SBHs may be detected when interacting with interstellar matter (Shvartsman, 1971; Agol & Kamionkowski, 2002; Barkov et al., 2012; Tsuna et al., 2018; Tsuna & Kawanaka, 2019), however, the estimated X-ray emission from an accreting SBH is very small ($<10^{35}\,\mathrm{erg}\,\mathrm{s}^{-1}$; e.g. Barkov et al., 2012). Such faint sources may be potentially detected only in our Galaxy. Matsumoto et al. (2018) suggested that the accretion onto a SBH from interstellar medium may not be necessarily spherical and result in the formation of an accretion disk. In such disks the same instabilities may arise as in X-ray Novae resulting in a transient behavior with outburst luminosities reaching $\sim 10^{38}\,\mathrm{erg}\,\mathrm{s}^{-1}$.

SBHs can be also detected as a gravitational microlenses (Agol et al., 2002; Wyrzykowski et al., 2016; Rybicki et al., 2018). As shown in this paper, SBHs may have very high velocities what increases a chance for microlensing event. On the other hand, a significant part of DCOs ($\mathcal{R}_{\rm DCO,2}$ channel) have separations small enough to act as a single, very massive lens, with the mass as the sum of the masses of individual components. In combination with small proper motions of such binaries, they may generate long-lasting (t${}_{\rm{E}}\sim 1\,\mathrm{yr}$) events, however, we note that undisrupted compact DCOs make less than 15% of all BHs in the Milky Way.

Although none SBHs were detected so far, a few observational methods were proposed. The predicted number of observations in future X-ray missions is strongly dependant on the BH velocities and accretion models. Future missions may help to improve the constrains on this physics. Nevertheless, detectability of SBHs in X-rays is strongly spoiled by background AGNs, hard coronal emitters and CVs (Motch & Pakull, 2012). Another possibility, to which more attention is put in this work, is the detection of SBHs with microlensing (Paczynski, 1986). First candidates were proposed in Bennett et al. (2002) and Mao et al. (2002) and were followed-up with X-ray observations, however, no signal from ISM accretion was detected (Maeda et al., 2005; Nucita et al., 2006). Recently, several candidates were proposed (Wyrzykowski et al., 2016), but the lens’s mass estimates are heavily degenerated due to the lack of measurement of relative velocities. Consequently, WDs or NSs cannot be fully ruled out. For the currently on-going events, there is an opportunity for degeneracy breaking thanks to astrometric measurements from optical interferometry (Dong et al., 2019) or Gaia(Rybicki et al., 2018), however, this will only be possible for the brightest events. In near future with the LSST it should be possible to observe thousands of microlensing events due to BHs and, therefore, to confront the stellar evolution predictions with the observed parameters of Galactic BHs.

5.2. Detectability of BHs in binaries

All dynamically confirmed galactic BHs were detected in XRBs. We note that a BH candidate was also detected in a non-interacting binary (Thompson et al., 2018) and many more are expected to be seen by Gaia (Mashian & Loeb, 2017; Breivik et al., 2017; Yamaguchi et al., 2018). In XRBs a compact object accretes mass from its companion which results in a formation of an accretion disk and production of highly energetic radiation. If the companion is observable, radial velocities can be measured and the mass function can be estimated. Unfortunately, the majority of XRBs containing BHs are transient systems, thus are visible only during outbursts and only a few recurring systems were observed. What is more, typical transient systems are characterized by low mass donors ($M_{\mathrm{don}}\lesssim 2\,M_{\odot}$), therefore, undetectable from extra-galactic distances and outshined by the disk during an outburst. In $19$ XRBs, BHs were dynamically confirmed (e.g. Wiktorowicz et al., 2014), however, there are many more candidates (Corral-Santana et al., 2016).

Gravitational waves emission from double BH mergers may be perceived as a detection of BHs. Up-to-now, 10 binary BH merger events were detected by LIGO/Virgo (The LIGO Scientific Collaboration et al., 2018b). Using different pipeline Venumadhav et al. (2019) found 6 additional sources in publicly available data from the second observing run. Notably, in spite of one double NS merger detection (Abbott et al., 2017b), no BH+NS mergers were detected. The masses of pre-merger BHs range from $\sim 7.7\,M_{\odot}$ to $\sim 50.6\,M_{\odot}$ with the total binary mass as high as $\sim 85.1\,M_{\odot}$ (The LIGO Scientific Collaboration et al., 2018b). Many tens of detections are expected from the next observational runs (The LIGO Scientific Collaboration et al., 2018a) and several candidates have already been found (e.g. Singer, 2019).

BHs in binaries may be detected also through microlensing. They may lens differently depending on the compactness of the binary. If the separation is significantly smaller than the total mass Einstein radius, the system will act as a single lens with a mass equal to the sum of the components’ masses. On the other hand, if the system is very wide (many Einstein radii), a BH will act as a single lens unaffected by a remote companion. In-between these two options, we have systems that act as a binary lens with all the effects connected with caustics crossing. Up-to-now, no microlensing events with BH lenses were confirmed. Although there are several candidates (e.g. Miyake et al., 2012; Shvartzvald et al., 2015; Wyrzykowski et al., 2016), parameter degeneracies make it impossible to infer mass, thus distinguish between a BH and other low-luminosity objects like WDs.

5.3. Comparison with previous estimates of microlensing rates by BHs

Our prediction of $\sim 14$ microlensing events due to BHs in a year stays in contrast to previous estimates. Gould (2000) analytically calculated that only $1\%$ of microlensing events in the direction of the Bulge are due to BH lenses, what gives about 30 events in OGLE-III (3-4 per year). There are several reasons for this contradiction. First of all, Gould (2000) concentrates on BH lenses in the bulge, whereas our lenses are in the disk. $\pi_{\rm rel}$ (see Eq. 5) is typically smaller for bulge lenses, thus also the Einstein radius, which is proportional to the probability of microlensing. Secondly, Gould (2000) assumed $40\,M_{\odot}$ as a minimal ZAMS mass producing a BH remnant, whereas in our calculations we take $\sim 22\,M_{\odot}$ (for solar metallicity). Consequently, they obtained much less BHs from the same initial population. Moreover, BHs in Gould (2000) follow the velocity distribution of stars, whereas in our simulations BHs obtain additionally a kick (including NK, Blaauw kick, and the kick from binary disruption). Thirdly, the binary systems were ignored in the work of Gould (2000), while here we show that the binary systems could potentially contribute significantly to the population of lensing BHs.

A more recent work, Rybicki et al. (2018) used the same population synthesis StarTrack code to estimate the microlensing rate events due to BHs. However, they focus on astrometric microlensing events, i.e. events in which not only photometric observation is possible, but also astrometric signal can be detected with Gaia space mission. They estimated there should be a few such events per year. About an order of magnitude of difference between our results stems from the fact that only $\sim 5\%$ of bulge stars ($5\times 10^{6}$ in their work) are bright enough to potentially generate a detectable astrometric microlensing event. In contrast, we included all bulge stars as potential sources ($1.5\times 10^{8}$), thus $30$ times more, which results in significantly higher estimated rate. Moreover, Rybicki et al. (2018) used averaged, or typical values to calculate the rates, whereas in this paper we populate the Galaxy with stars, evolve them to the formation of a BH and calculate the rates directly.

We note that in our approach to estimate the microlensing rates we made a number of simplification. Especially, we used only one model from 8 analysed in this papers. Additionally, we chose a simple model of BH distribution in the Galaxy. For example, model of (Robin et al., 2003), which was used in Rybicki et al. (2018), gives microlensig rate by BH of $\sim 9$ yr*-1*, so nearly two times smaller than the rate obtained for the distribution described by Eq. 4. We also note that the rate may change by a factor of $\sim 5$ depending on the chosen evolutionary parameters, e.g. IMF, natal kicks, metallicity.

5.4. Hertzspurng Gap donors in common envelope evolution

It is still not known if binaries undergoing CE events with HG donors will merger due to the lack of clear core-envelope boundary, or they can survive, because the boundary develops during the HG phase. In order to check the importance of our assumption (model B - CE with HG donor always leads to a merger), we have tested also an opposite case (model A) where binaries are allowed to survive the CE phase even if donor is on its HG. The results are presented in Tab. 7.

The difference in the total number of BHs is small ($\lesssim 10\%$), as most of the BHs form without a CE phase during the binary evolution. More pronounced is the higher number of DCOs (particularly, the close ones), which may exceed two orders of magnitude in comparison to model B for higher NK models (NKR and NKBE). A part of close DCOs progenitors evolve through a CE phase and in model A donors are frequently HG stars.

The number of mergers leading to the formation of a BH is also affected. In model B, CE events with HG donors are treated as mergers, what is not always a case in model A. In a consequence, there are significantly fewer (up to $33\%$) BHs originating from stellar mergers in model A than in model B.

As the total number of BHs is affected only slightly by the treatment of HG donors in CE events, we predict a negligible impact on the microlensing observations. However, it may be significant if events caused by close systems are considered.

5.5. Binary fraction

In this work, we assumed that all stars heavier than $10\,M_{\odot}$ are born in binaries. Such a high binary fraction ($100\%$) for massive stars is supported by observations, which suggest their binary fraction to be higher than $90\%$ (Sana et al., 2012). For lower mass stars we assumed equal number of binaries and single stars on ZAMS (binary fraction equal $50\%$).

Nonetheless, we agree that the binary fraction may not change in such a drastic way when the mass of the primary increases. More probably, the binary fraction increases continuously. Therefore, a fraction of massive stars, which are heavy enough to produce BHs ($M_{\rm ZAMS}\gtrsim 20\,M_{\odot}$) may actually form without companions. If included in our simulations, these stars could increase the relative fraction of single BHs in the entire BH population. However, the total number of BHs is expected to be lower due to the fact that flat mass ration distribution in binaries gives higher average mass of two stars in a binary, thus more BH progenitors, than two single stars following the IMF relation. We note, that mergers of binary systems of two massive stars (both heavy enough to form a BH in single stars evolution), may decrease the total number of BHs, but mergers of NS progenitors ($M_{\rm zams}\approx 8$–$20\,M_{\odot}$) may become heavy enough to form a BH in posterior evolution, thus increasing the total number of BHs.

Many of the massive binaries may be part of triples and higher-order systems (e.g. Toonen et al., 2016). Even $\sim 50\%$ of massive (OB type) stars may exist in triples (e.g. Sana et al., 2014), but the inner binary in a hierarchical triple system may evolve effectively isolated from the third star, thus estimating the significance of higher-order systems on BH populations is complicated. Of particular interest are mergers induced by Kozai-Lidov mechanism in stellar triples (e.g. Antonini et al., 2017). However, higher-order systems are not understood yet well enough to be included in population synthesis modeling.

6. Summary

In this paper we analysed the general properties of synthetic BH populations in different stellar environments represented by models which differ in the metallicity, initial parameter distributions, and natal kick prescriptions. We note that the results are applicable for further studies like predictions for present and future survey missions, or in-depth analysis of specific BH populations (e.g. the Milky Way galaxy). We particularly focused on BHs originating from disrupted binaries and stellar mergers, which were frequently omitted in previous studies. We find that those BHs constitute a majority of the total BH population.

Particularly, we show that single BHs dominate the total population of BH ($\gtrsim 83\%$ of all BHs), even though massive stars form predominantly in binaries. Both binary disruptions and stellar mergers are important and the predicted number of BHs is only slightly affected (up to a factor of $\sim 3$) by a chosen model. Although BHs in binaries constitute only a small part ($\lesssim 17\%$) of the population, their number is heavily dependent (about two orders of magnitude) on the adopted model parameters (especially natal kick prescription and star formation history), therefore, new observations may significantly help to constrain these parameters and better understand evolutionary processes (e.g. Maccarone et al., 2019).

Using our results, we calculated the expected rate of microlensing events with BH lenses in the direction of the Galactic bulge. We expect that as many as $\sim 14$ such events per year with average crossing times around $100$ days in the OGLE-III footprint and about $26$ in the OGLE-IV data. Only some of these events may be observable from Earth due to technical limitation (e.g., low luminosity, extinction). In this estimate we have neglected BHs which remained in binaries, because, low fraction ($\lesssim 17\%$ of all BHs) and slow velocities (typically, $\lesssim 20$ km/s) allow to suppose that their influence is small.

A grid of 54 models (including 8 main models analysed in the paper) is available in a free-access database222https://universeathome.pl/universe/bhdb.php. In our future work we plan to utilize these datafiles in order to study microlensing by single BHs in the Galaxy and provide detailed predictions for surveys like Gaia, or Einstein Telescope.

We are thankful to the anonymous referee who helped to improve the paper and thousands of volunteers, who took part in the Universe@Home project333https://universeathome.pl/ and provided their computers for this research. We acknowledge the Polish NCN grant HARMONIA 2015/18/M/ST9/00544 to ŁW, OPTICON H2020 EC grant no. 730890, MNiSW grant DIR/WK/2018/12 and SONATA BIS 2 DEC-2012/07/E/ST9/01360. G.W. is partly supported by the President’s International Fellowship Initiative (PIFI) of the Chinese Academy of Sciences under grant no.2018PM0017 and by the Strategic Priority Research Program of the Chinese Academy of Science Multi-waveband Gravitational Wave Universe (Grant No. XDB23040000). KB and GW acknowledge NCN grant MAESTRO 2018/30/A/ST9/00050. MC and JK acknowledge support from the Netherlands Organisation for Scientific Research (NWO). We acknowledge the Polish NCN grant PRELUDIUM 2017/25/N/ST9/01253 to KR. This work is partly supported by the National Key Program for Science and Technology Research and Development (Grant No. 2016YFA0400704), and the National Natural Science Foundation of China (Grant No. 11690024 and 11873056).