Detecting Black Hole Binaries by Gaia

We begin by considering the distributions of the binary parameters at birth. For the primary star, we use the Kroupa (2001) IMF as the fiducial one:

$$\begin{eqnarray}{{\rm{\Psi }}}_{{\rm{K}}01}({\bar{M}}_{1})d{\bar{M}}_{1}\propto \left\{\begin{array}{ll}{\bar{M}}_{1}^{-1.3}d{\bar{M}}_{1}, & 0.08\,{M}_{\odot }\leqslant {\bar{M}}_{1}\lt 0.5\,{M}_{\odot },\\ {\bar{M}}_{1}^{-2.3}d{\bar{M}}_{1}, & 0.5\,{M}_{\odot }\leqslant {\bar{M}}_{1}\lt 100\,{M}_{\odot }.\end{array}\right.\end{eqnarray} \tag{ 1 }$$

We also take another IMF given by Kroupa et al. (1993) and Kroupa & Weidner (2003) to investigate dependence of results on IMF:

$$\begin{eqnarray}{{\rm{\Psi }}}_{{\rm{K}}93}({\bar{M}}_{1})d{\bar{M}}_{1}\propto \left\{\begin{array}{ll}{\bar{M}}_{1}^{-1.3}d{\bar{M}}_{1}, & 0.08\,{M}_{\odot }\leqslant {\bar{M}}_{1}\lt 0.5\,{M}_{\odot },\\ {\bar{M}}_{1}^{-2.2}d{\bar{M}}_{1}, & 0.5\,{M}_{\odot }\leqslant {\bar{M}}_{1}\lt 1.0\,{M}_{\odot },\\ {\bar{M}}_{1}^{-2.7}d{\bar{M}}_{1}, & 1.0\,{M}_{\odot }\leqslant {\bar{M}}_{1}\lt 100\,{M}_{\odot }.\end{array}\right.\end{eqnarray} \tag{ 2 }$$

For the secondary mass, we assume a flat mass ratio distribution as the fiducial case (Kuiper 1935; Kobulnicky & Fryer 2007):

$$\begin{eqnarray}&&{\rm{\Phi }}(q)=\displaystyle \frac{1}{(1-{\bar{M}}_{\min }/{\bar{M}}_{1})}{q}^{0}.\end{eqnarray} \tag{ 3 }$$

We also try the cases of the index of q, −1 and +1. Here, we set the lower limit of q as ${\bar{M}}_{\min }/{\bar{M}}_{1}$, where ${\bar{M}}_{\min }=0.08\,{M}_{\odot }$ is the minimal initial mass of a star.

We assume that the recent specific star formation rate in our Galaxy to be a constant (Belczynski et al. 2007),

$$\begin{eqnarray}&&\int d\bar{M}\,\bar{M}{\rm{\Psi }}(\bar{M})=3.5\,{M}_{\odot }\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 4 }$$

where Ψ = Ψ_K01 or Ψ_K93. This assumption is justified as the life times of stars that we examine are much shorter than the evolution times of our Galaxy.

The distribution of initial binary separations is assumed to be logarithmically flat (Abt 1983):

$$\begin{eqnarray}&&{\rm{\Gamma }}(\bar{A})=\displaystyle \frac{{{\rm{\Gamma }}}_{0}}{\bar{A}}.\end{eqnarray} \tag{ 5 }$$

The normalization factor Γ₀ is determined by the range of $\bar{A}$:

$$\begin{eqnarray}&&{\int }_{{\bar{A}}_{\min }}^{{\bar{A}}_{\max }}d\bar{A}{\rm{\Gamma }}(\bar{A})=1.\end{eqnarray} \tag{ 6 }$$

We set the lower limit of this integral, ${\bar{A}}_{\min }$, as the distance such that the primary fills its Roche lobe at the periastron (Eggleton 1983):

$$\begin{eqnarray}&&{\bar{A}}_{\min }=\displaystyle \frac{0.6{q}^{-2/3}+\mathrm{ln}(1+{q}^{-1/3})}{0.49{q}^{-2/3}}{R}_{1}.\end{eqnarray} \tag{ 7 }$$

With this definition of A_min and the normalization condition, Γ₀ is a function of ${\bar{M}}_{1}$ and q. We set the upper limit as ${\bar{A}}_{\max }={10}^{4}$ au. Recent studies show that the double stars with the separation up to ∼10⁴ au can stay gravitationally bound for long timescales (Andrews et al. 2017; Oelkers et al. 2017).

As we are interested in binaries that contain a BH (as a remnant of the primary star) and a secondary star that has not collapsed, the age of the system should be between t = t_L,1 and t = t_L,2, where t_L,i is the lifetime of a star with an initial mass ${\bar{M}}_{i}$. We adopt here the lifetime suggested by Eggleton (1983). For the primary star, t_L,1 is the time when it collapses into a BH.

Turning now to the current conditions of the system, we note that with no current MT from the secondary star to the BH the radius of the secondary should be smaller than its Roche lobe:

$$\begin{eqnarray}&&{R}_{L}({M}_{2}/{M}_{\mathrm{BH}},A)\gt {R}_{2}.\end{eqnarray} \tag{ 8 }$$

Using the mass–radius relation for terminal age main sequence (Demircan & Kahraman 1991) and the formula by Eggleton (1983), we can rewrite the condition (8) as

$$\begin{eqnarray}&&A\gt {A}_{\min }\\ &&\equiv \,\displaystyle \frac{0.6+{\left(\tfrac{{M}_{2}}{{M}_{\mathrm{BH}}}\right)}^{2/3}\mathrm{ln}\left[1+{\left(\tfrac{{M}_{2}}{{M}_{\mathrm{BH}}}\right)}^{-1/3}\right]}{0.49}{\left(\displaystyle \frac{{M}_{2}}{{M}_{\odot }}\right)}^{0.83}\,{R}_{\odot }.\end{eqnarray} \tag{ 9 }$$

This condition set the lower limit of the initial binary separation, $\bar{A}\gt {\bar{A}}_{\mathrm{RL}}$.

An upper limit, ${\bar{A}}_{\mathrm{period}}$, is determined by the condition that the orbital period of the binary should be shorter than P_max, which is given in Section 2.7:

$$\begin{eqnarray}&&A\lt {A}_{\max }\equiv {\left[\displaystyle \frac{G({M}_{1}+{M}_{2})}{{({P}_{\max }/2\pi )}^{2}}\right]}^{1/3}.\end{eqnarray} \tag{ 10 }$$

In the following sections, we related the current A_max to an upper limit ${\bar{A}}_{\mathrm{period}}$ on the initial separation.

We assume that the spacial distribution of BH binaries in the Galaxy traces the stellar distribution. According to Bahcall & Soneira (1980), we can write the star formation number density per unit mass bin in the Galactic disk as a function of the distance from the Galactic center r in the Galactic plane, and the distance perpendicular to the Galactic plane z as:

$$\begin{eqnarray}&&{\rho }_{d}({\boldsymbol{r}},z,\bar{M})={\rm{\Psi }}(\bar{M})\cdot {\rho }_{d,0}\exp \left[-\displaystyle \frac{z}{{h}_{z}}-\displaystyle \frac{{\boldsymbol{r}}-{r}_{0}}{{h}_{r}}\right],\end{eqnarray} \tag{ 11 }$$

where r₀ = 8.5  kpc is the distance from the Galactic center to the Sun, and (h_z, h_r) = (250 pc, 3.5 kpc) are the scale lengths for the exponential stellar distributions perpendicular and parallel to the Galactic plane, respectively. Hereafter, we consider only the disk component because the binaries in the Galactic bulge would not be observed due to the interstellar absorption. Then, we can determine the normalization factor, ρ_d,0 by

$$\begin{eqnarray}&&4\pi {\int }_{0}^{{{\boldsymbol{r}}}_{\max }}{\boldsymbol{r}}d{\boldsymbol{r}}{\int }_{0}^{{z}_{\max }}{dz}\,{\rho }_{d,0}\exp \left[-\displaystyle \frac{z}{{h}_{z}}-\displaystyle \frac{{\boldsymbol{r}}-{r}_{0}}{{h}_{r}}\right]=1.\end{eqnarray} \tag{ 12 }$$

We describe this distribution with respect to the spherical coordinate centered at the Earth, (D, b, l), where

$$\begin{eqnarray}&&{\boldsymbol{r}}={[{r}_{0}^{2}+{D}^{2}{\cos }^{2}b-2{{Dr}}_{0}\cos b\cos l]}^{1/2},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&z=D\sin b,\end{eqnarray} \tag{ 14 }$$

and D, b, and l are the distance from the Earth, Galactic latitude, and Galactic longitude.

The total number of BH binaries without mass accretion detectable by Gaia can be obtained as a multidimensional integral over the initial primary mass, the mass ratio, the initial separation, and the position:

$$\begin{eqnarray}N & = & \displaystyle \frac{2{f}_{\mathrm{bin}}}{1+{f}_{\mathrm{bin}}}{\displaystyle \int }_{{M}_{\min ,\mathrm{BH}}}^{100\,{M}_{\odot }}d{\bar{M}}_{1}{\displaystyle \int }_{{q}_{\min }}^{1}{dq}\,{{\rm{\Gamma }}}_{0}{\rm{\Phi }}(q)\\ & & \cdot ({t}_{{\rm{L}},2}-{t}_{{\rm{L}},1}){\displaystyle \int }_{\max ({\bar{A}}_{\mathrm{RL}},{\bar{A}}_{\min },{\bar{A}}_{\det })}^{\min ({\bar{A}}_{\mathrm{period}},{\bar{A}}_{\max })}d\bar{A}\,\displaystyle \frac{1}{\bar{A}}\\ & & \times {\displaystyle \int }_{0}^{2\pi }{dl}{\displaystyle \int }_{0}^{\pi /2}\cos {bdb}{\displaystyle \int }_{0}^{{D}_{\max }}{D}^{2}{dD}{\rho }_{d}(D,b,l,{\bar{M}}_{1}),\end{eqnarray} \tag{ 15 }$$

where we set the lower limit of the integration with respect to ${\bar{M}}_{1}$ as M_min,BH = 20 M_⊙, above which a primary star would form a BH after its collapse, and f_bin is the binary fraction. Hereafter, we assume f_bin = 0.5, which means that we have 50 binaries and 50 single stars out of 150 stars. In addition, q_min represents the minimum mass ratio defined in Section 2.6, and ${\bar{A}}_{\det }$ is defined by considering the condition for the BH identification with Gaia (Section 2.7). D_max is set as 10 kpc in Section 3. Note that while the integration looks simple, it involves numerous implicit dependences, e.g., t_L,i (i = 1, 2) depend on ${\bar{M}}_{i}.$

As shown below, the final binary separation, A, can be described as

$$\begin{eqnarray}&&A=\bar{A}\cdot a(q),\end{eqnarray} \tag{ 16 }$$

where a(q) is a function of the initial mass ratio q; then, we can describe the differential number distribution of BH binaries of interest as

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dM}}_{\mathrm{BH}}{{dM}}_{2}{dAdR}}=\displaystyle \frac{{dN}}{d{\bar{M}}_{1}d{\bar{M}}_{2}d\bar{A}{dR}}\cdot \displaystyle \frac{\partial ({\bar{M}}_{1},{\bar{M}}_{2},\bar{A})}{\partial ({M}_{\mathrm{BH}},{M}_{2},A)}\\ &&\quad =\,\displaystyle \frac{2{f}_{\mathrm{bin}}{{\rm{\Gamma }}}_{0}}{1+{f}_{\mathrm{bin}}}\cdot 2{\displaystyle \int }_{0}^{2\pi }\cos {bdb}{\displaystyle \int }_{0}^{\pi /2}{{dlR}}^{2}{\rho }_{d}(R,b,l,{\bar{M}}_{1})\\ &&\qquad \cdot \,\displaystyle \frac{1}{{\bar{M}}_{1}}({t}_{{\rm{L}},2}-{t}_{{\rm{L}},1})\displaystyle \frac{1}{\bar{A}}\cdot \displaystyle \frac{\partial ({\bar{M}}_{1},{\bar{M}}_{2},\bar{A})}{\partial ({M}_{\mathrm{BH}},{M}_{2},A)},\end{eqnarray} \tag{ 17 }$$

where $\partial ({\bar{M}}_{1},{\bar{M}}_{2},\bar{A})/\partial ({M}_{\mathrm{BH}},{M}_{2},{A}_{f})$ is the Jacobian of the variable transformation from the initial to final parameters.

We assume for simplicity that the BH mass satisfies the equation

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}=k\bar{{M}_{1}},\end{eqnarray} \tag{ 18 }$$

where we adopt k = 0.2 as the fiducial case. This is a conservative assumption. A good approximation is given by Equation (2.42) in Eggleton (2011). In short, k increases with mass from about 0.1 at M₁ = 0.8 ${M}_{\odot }$ to about 0.4 at M₁ = 40 ${M}_{\odot }$, and then it flattens toward 0.5 at very large masses. In reality, the WR stars will loose a lot of matter in the wind later, and the BH will have a smaller mass than the mass of kM. Our choice of k = 0.2 takes into account these two processes. We also try k = 0.1 and 0.5.

However, the real relation between the ZAMS mass and BH mass should be more complicated. Thus, we assume the mass relation as follows:

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}=\displaystyle \frac{2}{\mathrm{ln}3}\mathrm{ln}({\bar{M}}_{1}-19)+2,\end{eqnarray} \tag{ 19 }$$

where we use the mass relation in Belczynski et al. (2008) as a reference, and this function satisfies ${M}_{\mathrm{BH}}=2\ \mathrm{for}\ \bar{{M}_{1}}=20$ and ${M}_{\mathrm{BH}}=10\ \mathrm{for}\ \bar{{M}_{1}}=100$. The model with this mass relation is named "curved."

A general scenario for formation of binaries with a BH involves a binary that initially contains a massive star and a companion. We assume that the initial orbit is circular. The more massive star—the BH progenitor—will evolve faster and will initiate the first MT when it sufficiently expands and fills the Roche lobe. Here, we adopt $3\times {10}^{3}\,{R}_{\odot }$ as the maximum radius of the stars with ${\bar{M}}_{1}\gt 20\,{M}_{\odot }$. This value is given by the formula of the radius of the asymptotic giant branch star with mass of 20 ${M}_{\odot }$ and the luminosity of ${10}^{5.5}\,{L}_{\odot }$ (Hurley et al. 2000). Although the radius depends on the mass and luminosity of a star, the uncertainty is no more than a factor of a few, which has little influence on the result. We discuss this issue in detail in Section 4. We discuss the outcome of the MT in detail for two cases of the mass ratio in Sections 2.4 and 2.5.

In our calculation, we make several simplifying assumptions. We neglect the wind mass loss from the stars. The wind mass loss in the pre-MT phase will lead to tightening of the orbit, and therefore it is degenerate with the initial orbital separation. The mass loss from the primary will also decrease the mass that can be transferred to the companion in the large mass ratio case. However, the contribution to the observed number of binaries in this case is small, which is discussed in detail in Section 4. We assume that a BH forms through direct collapse of the compact core with no mass loss during the process. Additionally, we assume that the BH receive no natal kicks during formation. The current understanding of BH formation is that such kicks are small (Section 4).

If the mass ratio is larger than ≳0.3, the MT from a giant primary star will initially be unstable, but then it will stabilize because the mass ratio will be reversed. In such a MT, the orbit initially tightens (i.e., the orbital separation decreases) and the system loses its mass rapidly. Once the mass ratio of the binary reaches unity, the separation starts to increase and the MT finally ceases. Using the formulae presented by Podsiadlowski et al. (1992), the ratio of the orbital separation to the initial one is given by

$$\begin{eqnarray}&&a(q)=\displaystyle \frac{{M}_{\mathrm{BH}}+{M}_{2}}{{\bar{M}}_{1}+{\bar{M}}_{2}}{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{\bar{M}}_{1}}\right)}^{{c}_{1}}{\left(\displaystyle \frac{{M}_{2}}{{\bar{M}}_{2}}\right)}^{{c}_{2}},\end{eqnarray} \tag{ 20 }$$

where c₁ = α(1 − β) − 2 and c₂ = − α(1 − β)/β − 2. Here, α is the specific angular momentum per unit mass lost from the system, and β is the fraction of mass that goes to the acceptor from the donor. We assume that the MT would continue until the primary mass is equal to the mass of its remnant (i.e., a BH). After the MT, the mass of the secondary becomes ${M}_{2}={\bar{M}}_{2}+\beta (1-k){\bar{M}}_{1}$, and then we obtain

$$\begin{eqnarray}&&a(q)=\displaystyle \frac{(k+\beta (1-k))+q}{1+q}{k}^{{c}_{1}}{\left(\displaystyle \frac{\beta (1-k)}{q}+1\right)}^{{c}_{2}}.\end{eqnarray} \tag{ 21 }$$

Hereafter, we use α = 1.0 and β = 0.5, corresponding to the standard evolution model (e.g., Belczynski et al. 2002). In this case, the Jacobian (Equation (17)) is

$$\begin{eqnarray}&&\displaystyle \frac{\partial ({\bar{M}}_{1},{\bar{M}}_{2},\bar{A})}{\partial ({M}_{\mathrm{BH}},{M}_{2},{A}_{f})}=\displaystyle \frac{1}{k\cdot a(q)}\\ &&\quad =\displaystyle \frac{1+q}{(k+\beta (1-k))+q}{k}^{-{c}_{1}-1}{\left(\displaystyle \frac{\beta (1-k)}{q}+1\right)}^{-{c}_{2}}.\end{eqnarray} \tag{ 22 }$$

If the mass ratio is smaller than ≲0.5, the system will undergo a violent MT from the primary to the secondary. This will leads to a CE phase. During the CE phase, the binary orbital energy is used to expel the envelope. The evolution of the binary separation can be modeled following Webbink (1984). Let us assume that the primary star has a core mass M_c,1 and an envelope mass M_env,1, and that the initial and final orbital separations are A_i and A_f, respectively. From energy conservation, we have

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CE}}\left(\displaystyle \frac{{{GM}}_{c,1}{M}_{2}}{2{A}_{f}}-\displaystyle \frac{{{GM}}_{1}{M}_{2}}{2{A}_{i}}\right)=\displaystyle \frac{{{GM}}_{1}{M}_{\mathrm{env},1}}{\lambda {R}_{L}(1/q,{A}_{i})},\end{eqnarray} \tag{ 23 }$$

where α_CE is the efficiency of converting the orbital energy into the kinetic energy of an envelope during the CE phase, and λ is a parameter which is determined by the structure of the primary star. Following Belczynski et al. (2002) and Belczynski et al. (2008), we assume that α_CEλ = 1 as the fiducial case. For the massive star with 20 ${M}_{\odot }$ or larger, λ can be ∼0.1 if the stellar radius is larger than ∼1 au (Dominik et al. 2012). Thus, the cases that α_CEλ = 0.1 are also investigated in Section 3. We also examine the case α_CEλ = 2.0. Then the binary separation shrinks by a factor of

$$\begin{eqnarray}&&a(q)={\left[\displaystyle \frac{2(1-k)}{{\alpha }_{\mathrm{CE}}\lambda {r}_{L}{kq}}+\displaystyle \frac{1}{k}\right]}^{-1},\end{eqnarray} \tag{ 24 }$$

where r_L = R_L(q₁, A_i)/A_i. For the binary not to merge after the CE phase, the final separation should be larger than the sum of the stellar radii ${A}_{f}\gt {R}_{1}^{{\prime} }+{R}_{2}^{{\prime} }$, where ${R}_{i}^{{\prime} }$ (i = 1, 2) are the stellar radii right after the CE phase. In this case, the Jacobian (Equation (17)) is

$$\begin{eqnarray}&&\displaystyle \frac{\partial ({\bar{M}}_{1},{\bar{M}}_{2},\bar{A})}{\partial ({M}_{\mathrm{BH}},{M}_{2},{A}_{f})}=\displaystyle \frac{1}{k\cdot a(q)}=\displaystyle \frac{1}{k}\left[\displaystyle \frac{2(1-k)}{{r}_{L}{kq}}+\displaystyle \frac{1}{k}\right].\end{eqnarray} \tag{ 25 }$$

Gaia is observing in optical wavelengths, therefore interstellar extinction reduces the total number of BH binaries detectable by Gaia. Gaia's limiting magnitude is 20 mag, and the average extinction of the Milky Way disk is ∼1 mag per 1 kpc in the V-band (see Spitzer 1978, and, for example, Shafter 2017). Thus, the fraction of stars detectable by Gaia in all stars may be drastically reduced for the distance farther from 1 kpc. As the lower limit of the luminosity of the observable star gets higher, the lower limit of the corresponding stellar mass gets higher. In this paper, we equate the Gaia band with the V-band, which is valid when the star is bluer than G-type stars whose color V – I ≲ 1 (Jordi et al. 2010).

The interstellar extinction in the V-band A_V affects the relation between the absolute magnitude M_V, the apparent magnitude m_V, and the distance D:

$$\begin{eqnarray}&&{M}_{V}={m}_{V}-5(2+{\mathrm{log}}_{10}{D}_{\mathrm{kpc}})-{A}_{V}({D}_{\mathrm{kpc}}),\end{eqnarray} \tag{ 26 }$$

where the distance is normalized by 1 kpc and we adopt A_V(D_kpc) = D_kpc. In addition, the absolute magnitude is related to the companion mass, M₂, as

$$\begin{eqnarray}{M}_{2}=\left\{\begin{array}{ll}{10}^{-0.1({M}_{V}-4.8)} & {M}_{V}\ \lt 8.5\\ 1.9\times {10}^{-0.17({M}_{V}-4.8)} & {M}_{V}\ \gt 8.5,\end{array}\right.\end{eqnarray} \tag{ 27 }$$

where we adopt the empirical mass–luminosity relation in Smith (1983). The absolute magnitude ${M}_{V}$ = 8.5 corresponds to a companion mass M₂ = 0.4 ${M}_{\odot }$. By combining Equations (26) and (27), we find the companion mass whose apparent magnitude and distance are m_V and D_kpc, respectively. Thus, given the limiting magnitude of Gaia and distance of the system, we obtain M_2,min, which is defined to be the minimum mass of the companion observable by Gaia. Therefore, q_min in Equation (15) is now defined as $\max (\bar{{M}_{2}}({M}_{2,\min })/\bar{{M}_{1}},{\bar{M}}_{\min }/\bar{{M}_{1}})$.

We need impose constraints on various parameters of the binaries to identify the primary of the binary system as a BH. The robust way to do so is to measure its mass. The astrometric observations of the companion star enables us to estimate the mass of the other unseen object through the measurements of the semimajor axis and the orbital period. The mass M_BH can be expressed by M₂, the orbital period P_orb, and the angular semimajor axis, a_*:

$$\begin{eqnarray}&&\displaystyle \frac{{({M}_{\mathrm{BH}}+{M}_{2})}^{2}}{{M}_{\mathrm{BH}}^{3}}=\displaystyle \frac{G}{4{\pi }^{2}}\displaystyle \frac{{P}_{\mathrm{orb}}^{2}}{{({a}_{* }D)}^{3}},\end{eqnarray} \tag{ 28 }$$

where G is the gravitational constant. This equation means that the identification of BH requires measurements of M₂, P_orb, a_*, and D with a sufficient accuracy. In what follows, we estimate the required standard errors of these parameters for the BH identification and constraints on these parameters.

If the mass of the hidden companion is larger than 3 ${M}_{\odot }$ with a n–σ confidence level, the object can be identified as a BH. Note that 3 ${M}_{\odot }$ is the fiducial minimum mass of BH expected from the maximum mass of neutron stars (Kalogera & Baym 1996). This condition is expressed as

$$\begin{eqnarray}&&{M}_{\mathrm{BH}}-n{\sigma }_{\mathrm{MBH}}\gt 3\,{M}_{\odot },\end{eqnarray} \tag{ 29 }$$

where σ_MBH is a standard error of the mass estimate of the unseen primary, and we adopt n = 1. Using Equation (28) σ_MB/M_BH is related to σ_M2, σ_P, σ_a, and σ_D (the standard errors of the companion mass, orbital period, semimajor axis, and distance, respectively) as

$$\begin{eqnarray}{\left(\displaystyle \frac{{\sigma }_{\mathrm{MB}}}{{M}_{\mathrm{BH}}}\right)}^{2} & = & {\left(\displaystyle \frac{3}{2}-\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\mathrm{BH}}+{M}_{2}}\right)}^{-2}\left[{\left(\displaystyle \frac{{M}_{2}}{{M}_{\mathrm{BH}}+{M}_{2}}\right)}^{2}\displaystyle \frac{{\sigma }_{{\rm{M}}2}^{2}}{{M}_{2}^{2}}\right.\\ & & \left.+\displaystyle \frac{{\sigma }_{P}^{2}}{{P}_{\mathrm{orb}}^{2}}+\displaystyle \frac{9}{4}\left(\displaystyle \frac{{\sigma }_{a}^{2}}{{a}_{* }^{2}}+\displaystyle \frac{{\sigma }_{D}^{2}}{{D}^{2}}\right)\right],\end{eqnarray} \tag{ 30 }$$

where we assume that for all parameters the ratios of the standard errors to the parameters themselves are smaller than 1, and the correlation between errors of these parameters can be neglected. From Equations (29) and (30), we can constrain the ratios of the standard errors to the parameters, σ_M2/M₂, σ_P/P_orb, σ_a/a_*, and σ_D/D. If the following constraints are satisfied, Equation (29) is also satisfied, when M_BH > 5 ${M}_{\odot }$:

$$\begin{eqnarray}&&{M}_{2}\lt 0.1{\sigma }_{{\rm{M}}2},{P}_{\mathrm{orb}}\lt 0.1{\sigma }_{P},{a}_{* }\lt 0.1{\sigma }_{a},\mathrm{and}\,\,D\lt 0.1{\sigma }_{D}.\end{eqnarray} \tag{ 31 }$$

When M_BH < 5 ${M}_{\odot }$, conditions more stringent than in Equation (31) are required for satisfying Equation (29). However, empirically few BHs weigh less than 5 ${M}_{\odot }$ (Özel et al. 2010), so that we expect that the conditions in Equation (31) is useful for identifying most BHs. We note that these condition equations (Equation (31)) are given to simplify the following reduction, and therefore just necessary conditions.

The constraints on the standard errors of orbital period, semimajor axis, and distance in Equation (31) are reduced to conditions on BH binaries detectable by Gaia. The constraint on the standard error of the distance imposes a condition on the magnitude of the companion at a distance D. Because the distance is inversely proportional to the parallax π, σ_D/D can be expressed as σ_π/π, where σ_π is the standard error of the parallax for the Gaia astrometry, which is related to the apparent magnitude of the Gaia band (de Bruijne 2012). Assuming that the apparent magnitude of the Gaia band is equal to V-band magnitude as done in the previous section, the condition σ_π/π is represented as

$$\begin{eqnarray}&&{(-1.631+680.8\cdot z({m}_{V})+32.73\cdot {z}^{2}({m}_{V}))}^{1/2}\lt \displaystyle \frac{{10}^{2}}{{D}_{\mathrm{kpc}}},\end{eqnarray} \tag{ 32 }$$

where the expression of σ_π appears in Gaia Collaboration et al. (2016), and the function z(m_V) is

$$\begin{eqnarray}&&z({m}_{V})={10}^{0.4(\max [12.09,{m}_{V}]-15)}.\end{eqnarray} \tag{ 33 }$$

In addition, we neglect the factor including V − I of the original expression of σ_π(G) because this factor changes σ_π only by a few percent. We note that Equation (32) gives the maximum detectable apparent V magnitude m_V,max for a fixed distance as

$$\begin{eqnarray}&&{m}_{V,{\rm{\max }}}\sim 17.5+2.5{\mathrm{log}}_{10}[\sqrt{1+{(0.6{D}_{\mathrm{kpc}})}^{-2}}-1],\end{eqnarray} \tag{ 34 }$$

where we assume that the maximum function in Equation (33) is simply equal to m_V, which is valid for D_kpc ≲ 10. This maximum magnitude enables us to obtain the minimum companion mass M_2,min by using Equations (26) and (27), that is, M_2,min = M₂(m_V,max).

The conditions for the semimajor axis and orbital period in Equation (31) constrain the range of the semimajor axis. The standard error of the semimajor axis of the stellar orbit σ_a is expected to be similar to σ_π, because the semimajor axis of the stellar orbit a_* is roughly the size of the orbit on the celestial sphere. Thus, we can assume σ_a ∼ σ_π, and therefore, the constraint on the semimajor axis of the binary system A can be written as

$$\begin{eqnarray}&&A\gt 10\displaystyle \frac{{M}_{\mathrm{BH}}+{M}_{2}}{{M}_{\mathrm{BH}}}{\sigma }_{\pi }({m}_{V})D\equiv {A}_{\mathrm{ast}}.\end{eqnarray} \tag{ 35 }$$

According to orbital solutions in the Hipparcos and Tycho catalog (ESA 1997), all binaries with the orbital period less than 2/3 of the mission period of Hipparcos show the standard error of the orbital period less than 1/10 of the orbital period. Thus, we expect that for BH binaries with the orbital period less than ∼3 years, the orbital period will be measured with Gaia at the standard error less than 10%. Therefore, we adopt 3 years as the maximum orbital period. The minimum orbital period might be determined by Gaia's cadence for each object, which is roughly 50 days, so that we adopt 50 days as the minimum orbital period. These conditions give a range of the semimajor axis:

$$\begin{eqnarray}&&A({P}_{\mathrm{orb}}=50\ \mathrm{days})\lt A\lt A({P}_{\mathrm{orb}}=3\ \mathrm{years}).\end{eqnarray} \tag{ 36 }$$

The typical standard error of the stellar mass is ∼10% (Tetzlaff et al. 2011), where stellar masses are measured by using their luminosities and temperature, so that we expect that the standard error of the stellar mass measured by Gaia is also ∼10%. Therefore, ${\bar{A}}_{\det }$ is defined to be $\max (\bar{A}({A}_{\mathrm{ast}}),\bar{A}(A({P}_{\mathrm{orb}}=50\ \mathrm{days})))$.

In the next section, we obtain the numbers of BH binaries detectable with Gaia by integrating Equation (17) for various models in which the parameters are different from each other. The parameters in models are shown in Table 1. We show just distributions of the BH mass for models other than the fiducial one. We have four quantities as integration variables. The integral range of M_BH are assumed to be [4, 30 M_⊙] in the case of the fiducial model. The minimal companion mass M₂, M_2,min, is given in Section 2.7. The maximal companion mass is $\bar{{M}_{1}}[1+\beta (1-k)]$ for a mass ratio larger than 0.5 or $\bar{{M}_{1}}$ for a mass ratio smaller than 0.5. The integral interval of the semimajor axis is the range such that Equations (35) and (36) are satisfied, where we note that the semimajor axis of all binaries should be between ${\bar{A}}_{\min }a(q)$ and ${\bar{A}}_{\max }a(q)$.

The total numbers of detectable BH binaries for the models, lin01, lin05, and curved, are shown in Table 2. We see that the total number correlates with the coefficient k. This can be interpreted to be an effect of astrometric observation. If BH mass is larger, the orbit of companion is also larger, which increases the detectability of BH by Gaia. This can be also understood by Equation (35). Thus, as the coefficient k increases, the total number increases.

Figure 2 shows distributions of BH mass when the relation between the ZAMS mass and BH mass is changed. In the model of the linear mass relation (fiducial, lin01, and lin05), the BH masses show the single power-law distributions, and their power-law indexes are similar to each other. On the other hand, the distribution of the model "curved" shows a peak at ∼7 ${M}_{\odot }$. This feature arises clearly due to the projection of IMF onto the BH mass with the relation of Equation (19). The small number below 7 ${M}_{\odot }$ reflects the steep slope in the function of Equation (19), where we note that the minimal BH mass is 2 ${M}_{\odot }$, which is within the mass range of a neutron star. The drastic decrease over 7 ${M}_{\odot }$ is caused by the projection of IMF, whose mass range is [40 ${M}_{\odot }$: 150 ${M}_{\odot }$] onto the narrow mass range, i.e., [7 ${M}_{\odot }$: 10 ${M}_{\odot }$].

Zoom In Zoom Out Reset image size

We also examine the IMF shown in Kroupa & Weidner (2003). The total number of detectable BH binaries is estimated to be ∼200 (Table 2). The power-law index of this IMF for stars more massive than 1.0 ${M}_{\odot }$ is smaller than that of the fiducial one, so that the total number of primary star whose main sequence mass is >20 ${M}_{\odot }$ is less than that of the fiducial case. This leads to a smaller number of the detectable BH binaries than in the fiducial case. This smaller power-law index also affects the distribution of the BH mass, which is shown in the left panel of Figure 3. The slope of the distribution in the case of K03 is steeper than that of the fiducial case. Thus, the power-law index of IMF directly affects the mass distribution of BH detectable with Gaia.

Zoom In Zoom Out Reset image size

The parameter αλ does not affect the total number of detectable BH binaries. This is because the detectable binaries do not undergo the CE phase. Therefore, the distributions of BH mass (Figure 3) show no difference between fiducial, al01, and al20 cases.

Table 2 also shows the total numbers of the cases in which the distribution of mass ratio q is different. This indicates that the smaller the power-law index of the q distribution, the smaller the total number. As the index becomes smaller, the fraction of massive stars becomes smaller. The massive stars contribute the total number because they undergo the MT phase, so that the total number decreases as a result. Figure 4 shows the corresponding distributions of BH mass, whose slopes are not so different from that of the fiducial case. This means that the q distribution does not affect the shape of the distribution of BH masses.

Zoom In Zoom Out Reset image size

The number estimate of detectable BHs depends on various parameters other than those changed in this paper. Among those are the star formation rate (Equation (4)), the distribution of semimajor axis (Equation (5)), the binaries spatial distribution (Equation (11)), the range of the orbital period required for the detection of BHs (Equation (36)), and Gaia's astrometric precision σ_π. Thus, the total number shown in Table 2 can easily vary if those value are different from the one that we used. However, we do not expect that reasonable variation of these parameters would change significantly the number of detected BHs.

Our calculations rely on several simplifying assumptions: we ignored the effects of the stellar wind mass loss, natal kicks, and the initial eccentricities. Before concluding, we review these assumptions here. First, we note that the orbital expansion due to the stellar wind mass loss does not affect the number count of detectable BH binaries significantly. The factor by which the orbital separation expands depends only on the ratio of the binary mass before mass loss to that after the mass loss, and it does not depend on the separation. Because we assume a logarithmically flat distribution of initial separation, such a modification of the orbital separation does not change the final count of detectable BH binaries. Second, the stellar wind mass loss may affect whether the binaries would experience a CE phase or a MT phase because the mass ratio would be modified from that in the ZAMS phase. However, when the mass is maximally lost from the primary star by the time of the CE or MT phase, the mass would be $k{\bar{M}}_{1}$, and then the initial mass ratio that divides two cases (CE or MT) would be no less than ≳0.3k. In the fiducial model, the distribution of the initial mass ratio is flat, and then such a modification would affect the final result only by a factor. Third, as for the effect of natal kicks, Figure 4 of Breivik et al. (2017) shows that kicks do not affect the final results such as the number of detectable BH binaries or its mass distribution. Finally, as for the initial eccentricity, Figure 4 of Breivik et al. (2017) shows that most of the detectable BH binaries have eccentricities that are nearly equal to zero. This means that our estimates obtained by assuming circular orbits for all binaries do not deviate so much from the estimates obtained taking into account the non-zero eccentricity.

While we assume the uniform binary fraction (50%) for entire stellar mass for simplicity, recent studies show that massive stars are expected to be found in binaries at higher rates (e.g., Sana et al. 2012; Moe & Di Stefano 2017). A larger binary fraction increases linearly the number of BH binaries. If we assume the binary fraction ∼0.7 (Sana et al. 2012), we can easily obtain the total numbers of BH binaries by multiplying those in Section 3 by ∼1.2.

We adopt a uniform value of $3\times {10}^{3}\,{R}_{\odot }$ as the maximal radius of massive stars. Although this value depends on the stellar mass (and the luminosity) and can be larger for more massive stars, a variation of this value does not affect the results. Our results show that the detectable binaries should have passed through a MT phase, which implies that the separation does not change so much after that phase for a binary with a massive companion. In addition, the orbital period of the binary whose Roche radius is $\sim 3\times {10}^{3}\,{R}_{\odot }$ is ∼10 yrs. This is larger than the period of the detectable binaries, which is given in Section 2.7. Thus, even if the maximal radius of massive stars is larger than $3\times {10}^{3}\,{R}_{\odot }$, Gaia cannot identify them as BH binaries, therefore the results are not influenced by this uncertainty.

We assume that the G-band magnitude is equal to the V-band magnitude. This is basically valid for nearby blue stars. Our results show that the companions of all detectable binaries are blue stars, but most of them are located at ∼7 kpc, which means that these stars suffer from significant interstellar extinction. The extinction can be estimated to be A_V ∼ 7 using "1 magnitude per kpc." The ratio of extinction A_G/A_V given in Jordi et al. (2010) shows that A_G/A_V ∼ 0.9 for the bluest stars (V − I ∼ 0.4; T_eff = 50,000 K in Jordi et al. (2010)) if A_V = 5. Correspondingly, we expect that, due to the linearity of A_G/A_V with respect to A_V, A_G/A_V ∼ 0.8 for A_V ∼ 7. Thus, we expect A_G ∼ 6. This difference in the extinction does not affect the results because the apparent V magnitude of these stars is ∼13 (${L}_{\mathrm{star}}\sim {10}^{5}$ and distance ∼7 kpc), which is much brighter than Gaia's limiting magnitude.

The semimajor axis and the parallax are not generally degenerate in the Gaia astrometric measurements. There are two reasons. One is that the periods of projected motion are different. The period of parallax is just 1 year and the orbital period is generally different from that. Thus, we can separate these two motions by the Fourier analysis. The other is that the phase and shape of these two motions are different. The phase and shape (ellipticity and position angle) of motion due to parallax is determined by the ecliptic longitude and latitude, respectively. On the other hand, those of the orbital motion are generally independent of the ecliptic coordinate. In this paper, we assume that the astrometric signature due to the orbital motion is large enough to be detected, so that these two motions can be separated. Of course, we note that if the orbital period is just the same as one year, and if the shape of orbital motion is the same as that of parallax, which occurs when the orbit is circular (the same as Earth orbit) and the inclination angle is equal to the ecliptic latitude, we expect that these two motions are degenerate. However, such case must be rare.