© Y. Wang et al. 2022

1 Introduction

The first space astrometric mission, Hipparcos, observed about 120 thousand sources with an effective distance up to 1 kpc in 1989 (Perryman et al. 1997). Of these sources, 18 thousand are recorded as non-single stars, 235 of which come with a binary orbit solution (Lindegren 1997). After more than 20 years, its successor Gaia has already enlarged the astrometry sample at least 12 thousand times, providing over 1.46 billion sources with full astrometric data and 344 million sources with only mean positions in the Gaia Early Data Release 3 (Lindegren et al. 2021).

There are two steps in measuring astrometric binaries, detecting the potential astrometric binary systems and then solving their orbits. In the first step the basic astrometric information of current releases of Gaia data are obtained by treating all of them as single stars (Lindegren et al. 2012). It is still possible to search for astrometric binaries from such astrometric solution by using renormalised unit weight error (RUWE) and excess_noise (Penoyre et al. 2020; Belokurov et al. 2020; Stassun & Torres 2021). Furthermore, using the proper motion anomaly, it is also possible to characterise the presence of companions with stellar or substellar mass (Kervella et al. 2019).

In the second step there have been various efforts to solve the orbital parameters of the binary systems using the Gaia astro-metric epoch data. Pourbaix (2011) suggests that Gaia could obtain the orbit solution for 8.8 million unresolved binaries. Gaia is also able to obtain the binary orbit solution for BHs, NSs, or white dwarfs with a luminous companion (BH/NS/WD-LCs) (Barstow et al. 2014), or even to grab the signal of exoplanets in some nearby sources (Casertano et al. 2008). Andrews et al. (2019; hereafter A19) made a realistic Markov chain Monte Carlo (MCMC) simulation of such BH/NS/WD-LC systems with mock Gaia epoch data. Their result shows that the orbit of BH-LCs can be solved up to about 1 kpc, while the uncertainty of NS-LCs can reach 0.1 M_⊙ at this distance.

To predict the number of Gaia solvable BH/NS-LC systems, people need not only the Gaia observation constraints, but also the model of binary evolution and Milky Way stellar population synthesis. Binary evolution codes like Massive Objects in Binary Stellar Evolution (MOBSE), StarTrack, and Stellar EVolution N-body (SEVN) could provide the evolution track of a binary system from birth to death (Giacobbo & Mapelli 2018, 2020; Belczynski et al. 2008; Spera et al. 2016). For Milky Way stellar populations, dozens of papers are listed in Wiktorowicz et al. (2020), Breivik et al. (2019), and Olejak et al. (2020; hereafter O20). Combining the Gaia observation constraints with the simulated Milky Way BH-LC population, Breivik et al. (2017) suggests that 3800–12,000 BH-LCs could be discovered by Gaia at the end of the five-year mission. Wiktorowicz et al. (2020) proposed that the number of Gaia solvable BH-LCs is in the range 41–340, after further studying Gaia’s observation capability. Other works not only give the prediction of dozens to hundreds of Gaia detectable BH/NS-LCs, but also compare the results between different evolution models (Gould & Salim 2002; Yalinewich et al. 2018; Yamaguchi et al. 2018; Breivik et al. 2017, 2019; Chawla et al. 2022; Shikauchi et al. 2020, 2022). With the MCMC simulation of Gaia observations, A19 applied their constraints to a simulated BH/NS-Giant population, and estimated that there are 74 BH-Gs and 190 NSGs with a relative precision better than 0.3, which are somewhat solvable in Gaia’s vision. This relative precision is defined by Eq. (14) in A19, which is parametrised by the angular size of the orbital separation and the number of Gaia observations for each source in a five-year mission.

In this paper, to study whether a BH/NS-LC source is solvable with Gaia observation, we focus on the second step. We used the MCMC simulation method to explore Gaia’s ability in a wider parameter range and studied the Gaia-solvable BH/NS-LCs with a simulated Galactic BH/NS-LC population. In Sect. 2 we try to construct a realistic observation data model for astrometric binary systems. We give a useful expression for the relative error of the dark companion, which can be used to test Gaia’s ability directly. In Sect. 3 we describe our simulation in detail, and provide a full error-prediction model. We apply it to our binary synthesis population in Sect. 4.3, which is constructed following the instructions in O20. In Sect. 5 we present our conclusions.

2 Mock astrometric data

First, we introduce the scanning observation mode of Gaia in Sect. 2.1 and the kinematics of binary systems in Sect. 2.2. Then, in Sect. 2.3 we derive the theoretical relation between the uncertainty of mass measurement and the observable variables. Finally, we describe the method of generating the mock data in Sect. 2.4.

2.1 Gaia observations

Gaia is a scanning satellite. When a star passes across these CCDs, the information of astrometry, photometry, and spectroscopy is recorded (Gaia Collaboration 2016). During a 5 yr mission, each source could have an average of 75 astrometry observations. The satellite scans the sky at a fixed speed following its own scanning law, which makes the observation uneven and different at each place in the sky. There are two directions of the CCDs, the along-scan direction (AL) and the across-scan direction (AC)¹. In Fig. 1 we provide a schematic diagram of the mock Gaia observation of the luminous star in a fictitious NS-LC system².

We followed the instructions in Appendix A of Perryman et al. (2014) that describes the stellar motion in the AL direction³ rather than the right ascension-declination (RA-Dec) coordinates used in A19. The method of Perryman et al. (2014) allows us to use Gaia observation auxiliary data directly. The displacement η of a source in the AL direction on the tangent plane can be described by the following equations: $$\begin{array}{lll}\eta \left(t\right)\hfill & =\hfill & \left[\text{Δ}{\alpha }_0\cos \delta +{\mu }_{\alpha }\cos \delta \cdot \left(t-t_0\right)+B\cdot X+G\cdot \gamma \right]\sin \theta \hfill \\ \hfill & \hfill & +\left[\text{Δ}{\delta }_0+{\mu }_{\delta }\cdot \left(t-t_0\right)+A\cdot X+F\cdot Y\right]\cos \theta \hfill \\ \hfill & \hfill & +{\displaystyle \prod _{\eta }\varpi .}\hfill \end{array}$$(1)

Here θ is the position angle of the scanning direction and Π_η is the along-scan parallax factor. The Gaia Observation forecast Tool (GOST) can give us the mock Gaia data of each point in the sky, containing θ, Π_η, and the time t. In Eq. (1) we used five parameters, [Δα₀ cos δ, Δδ₀,µ_α cos δ, ß_δ,ϖ], to describe the system motion. Except for the parallax ϖ, [Δα₀ cos δ, Δδ₀] stands for the start point along the RA-Dec direction, while [µ_α cos δ, µ_δ] gives the proper motion. The orbital motion of the luminous component in a BH/NS-LC system are calculated by six intermediate parameters [A, B, F, G, X, Y], which are introduced in Sect. 2.2.

It is possible to obtain the radial velocity (RV) data from the archive or follow-up observations; these data are very helpful to restrict the binary orbital parameters. Gaia itself is also able to provide enough RV data since it can obtain the astrometric, photometric, and spectropic data at the same time. At the end of the 4–7 CCD rows on the Gaia satellite is an area for collecting spectra (Gaia Collaboration 2016). Although a source (V < 15 mag) is supposed to get 40 transits on this area during the 5 yr observation period, the actual number might be smaller, due to many filters which decrease the number of spectra (Katz et al. 2019). So we analysed the sources that have RV data in Gaia DR2. The average brightness of these sources is V = 12 mag, whose precision could reach 1 km s⁻¹. The average number of astrometry observations of these sources within 22 months is 30, which corresponds to 8.3 in spectroscopic data. This leads to 20 spectra for sources that have 75 astrometry observations. In this paper the position we used only has 62 astrometry observations in the 5 yr mission, thus it would only get 17 spectra.

Fig. 1 Schematic diagram of the sky-projected orbit of the luminous component and the Gaia epoch data. The central blue cross is the barycentre. The blue line shows the binary orbit of the luminous star. The black points are the mock Gaia epoch data, which show offsets from the orbit. These offsets are caused by the uncertainty in AL direction, indicated by the orientation of the error bar. All the points omit the proper motion and the annual motion.

2.2 Binary kinematics

We used Thiele-Innes elements [A, B, F, G, X, Y] in Eq. (1) to describe the luminous secondary motion in an ellipse orbit around the barycentre. In Heintz (1978), the Thiele-Innes elements are defined by the following equations: $$\begin{array}{lll}A\hfill & =\hfill & a_s\cdot \left(\cos \omega \cdot \cos \text{Ω}-\sin \omega \cdot \sin \text{Ω}\cdot \cos i\right),\hfill \\ B\hfill & =\hfill & a_s\cdot \left(\cos \omega \cdot \sin \text{Ω}+\sin \omega \cdot \cos \text{Ω}\cdot \cos i\right),\hfill \\ F\hfill & =\hfill & a_s\cdot \left(-\sin \omega \cdot \cos \text{Ω}-\cos \omega \cdot \sin \text{Ω}\cdot \cos i\right),\hfill \\ G\hfill & =\hfill & a_s\cdot \left(-\sin \omega \cdot sin\text{Ω}+\cos \omega \cdot \cos \text{Ω}\cdot \cos i\right),\hfill \\ X\hfill & =\hfill & \cos E-e,\hfill \\ Y\hfill & =\hfill & \sqrt{1-e^2}\sin E.\hfill \end{array}$$(2)

Here a_s is the semi-major axis of the secondary expressed in radians, e is the eccentricity, i is the inclination, ω and Ω are used to describe the rotation between the true orbit and the projection orbit, P is the period, and T_p is the epoch of periastron passage. The parameter E is the eccentric anomaly used to describe the position of the secondary on the orbit, which is defined by $$E=e\sin E+\frac{2\pi }{P}\left(t-T_p\right).$$

The RV is mainly determined with [K, e, P, ω, γ], where K and γ are the Keplerian speed and the system radial velocity. The Keplerian speed K of the secondary is $$K=\frac{2\pi a_s\sin i}{P\varpi \sqrt{\left(1-e^2\right)}}.$$(3)

Then the RV can be calculated as $$RV=\gamma +K\left[\cos \left(f+\omega \right)+e\cdot \cos \omega \right],$$(4)

where f is the true anomaly: $$f=2\text{\hspace{0.17em}arctan}\left[\sqrt{\frac{1+e}{1-e}}\tan (\frac{E}{2})\right].$$

2.3 Mass measurement

If we solve the binary orbit from the mock data, then we can get the following relation for the total mass M_tot from the third Keplerian law: $$M_{\text{tot}}=m_p+m_s=\frac{4{\pi }^2{\left(R_p-R_s\right)}^2}{G{P}^2}.$$(5)

The subscripts p and s stand for primary and secondary, respectively, and R is the length of the semi-major axis of the star orbit (in AU). Combining Eq. (5) with the relation between R_s and R_p + R_s yields $$R_s=\frac{a_s}{\varpi }=\frac{m_p}{m_p+m_s}\cdot \left(R_P+R_s\right).$$

We then define a new parameter KF: $$KF\equiv \frac{m_p^3}{M_{\text{tot}}^3}=\frac{4{\pi }^2}{G{\varpi }^3}\frac{a_s^3}{p^2}.$$(6)

In Eq. (6) the expression $m_p^3{M}_{\text{tot}}^{-2}$ is proportional to $a_s^3{P}^{-2}$, which is similar to the right side of Eq. (5), the third Keplerian law, hence we call it the Keplerian factor KF. A19 studied the relative error of KF and used it as a criterion to show whether a BH/NS-LC system is solvable. If the error is smaller than 0.3, they think the system is solvable. However, KF contains both M_tot and mp, so its relative error is influenced strongly by the mass ratio $\frac{m_s}{m_p}$. Thus, using the error propagation method and Eq. (6), we derive the relative error of m_p, $$\frac{\sigma \left(m_p\right)}{m_p}=\frac{1}{m_p+3m_s}\cdot \left[\left(m_p+m_s\right)\cdot \Xi +2\sigma \left(m_s\right)\right],$$(7)

which is a more direct variable and easier to interpret. In Eq. (7) we use the notation σ(x) to present the error, which is half of the difference between the 16th and 84th percentiles, rather than the standard deviation σ_x. We also define a new variable $\Xi \equiv \frac{\left(KF\right)}{KF}$ for the relative error of KF for simplicity in the rest of this article.

In the simulation m_p, m_s, and σ(m_s) are fixed input parameters, while Ξ is obtained from directly observable binary parameters. Equation (7) helps us get the relative error of m_p.

2.4 Generating the mock data

To generate the mock data we needed to set the input parameters, generate the AL direction data using Eq. (1), and add the mock error to the data. We used Eq. (4) to generate the RV data.

Some of the input parameters always stay the same for simplicity, for example [Δα₀ cos δ, Δδ₀, ω, Ω], which are fixed at [0.1 mas,0.2mas, 11.5°, 11.5°]. Similar to the values of [0 mas, 0 mas, 30°, 30°] in A19, the values of these parameters are selected manually and kept the same for consistency in the whole study. Other parameters, such as [P, a_s, ϖ, e, i,ß], are varied in Sect. 3 to explore the parameter space. Here, β is the barycentric ecliptic latitude. A19 shows that different proper motion values might have different influence on a binary system. So we analysed a Gaia DR2 subsample containing all sources that satisfy 8 < m_G < 14, and obtained the mean proper motion when the parallax is [10, 1, 0.5, 0.25] mas. Here, m_G is the G band apparent magnitude of a source. At these parallax values we set the proper motion β as [56.8, 7.1, 4.9, 4.7] mas yr⁻¹, respectively, rather than the uniform 10 mas yr⁻¹ used in A19. The different components of proper motion [µ_α cos δ, µ_δ] are [µ · sin $\frac{\pi }{6}$, µ · cos $\frac{\pi }{6}$]. In order to use Eq. (1) to generate the AL direction data, we use GOST to get the observation information at [α₀, δ₀] = [6:11:49.07, 22:49:32.68]. It is the position of LB-1 (Liu et al. 2019), a wide star-black hole binary system candidate. Our initial purpose is simulating the astrometric motion of LB-1, and then extending it to a full simulation.

We added the error ρ_AL to η, following a normal distribution ρ_AL =N(0, σ_AL), where σ_AL is the standard deviation of η, mostly fixed at 0.1 mas for scalability. We only changed σ_AL in Sect. 3.6, where we studied the influence of σ_AL itself. Compared to the 0.1–10mas of σ_AL (Lindegren et al. 2021), A19 gives a mean error of 612 mas in the AC direction by analysing the epoch position error ellipses of the solar objects, which is hundreds of times the error in AL direction. Thus, we generated the mock data η and performed the simulation only in the AL direction.

We used Eq. (4) to generate the RV data, and added the RV error from a normal distribution, $$\text{RV}{\left(t\prime \right)}_{\text{new}}=\text{RV}\left(t\prime \right)+\mathcal{N}\left(0,{\sigma }_{\text{RV}}\right)\left(\text{km\hspace{0.17em}}{\text{s}}^{-1}\right),$$(8)

where t′ is the observation time of the spectrum. Gaia can also provide a large amount of RV data, and the RV error is varied with the apparent magnitude, with a median σ_RV around 1 km s⁻¹.

3 MCMC simulation

In this section we try to obtain a relation between Ξ and the other observable parameters with a series of MCMC simulations. First, we introduce the Bayesian model in our simulation in Sect. 3.1. Then we give the feasible period range in Sect. 3.2. In the following subsections we study the variation of Ξ with other parameters and construct a function, which is called the Ξ-relation for simplicity. A19 provide an analytic estimate for Ξ, $$\Xi \left(a_s,{\sigma }_{\text{AL}},N,N_0\right)=0.9\cdot \left(\frac{{\sigma }_{\text{AL}}}{a_s}\right){\left(\frac{N_0}{N}\right)}^{\frac{1}{2}},P<5\text{yr,}$$(9)

where P ≤ 5 yr and N₀ = 75. Here N is the number of Gaia observations, while N₀ is the number of observations used in the MCMC simulation. We not only re-examined the relation between Ξ and [a_s, σ_AL, N], but also explored the relation between Ξ and other parameters, such as [P, e, i, ß]: $$\begin{array}{lll}\Xi \hfill & =\hfill & \Xi \left(a_s,{\sigma }_{\text{AL}},N,P,e,i,\beta \right)\hfill \\ \hfill & =\hfill & {\text{Φ}}_0\left(a_s,P,\right)\cdot {\text{Φ}}_1\left(P\right)\cdot {\text{Φ}}_2\left(N,N_0\right)\hfill \\ \hfill & \hfill & {\text{.Φ}}_3\left({\sigma }_{\text{AL}}\right)\cdot {\text{Φ}}_4\left(P,e,i\right)\cdot {\text{Φ}}_5\left(\beta \right).\hfill \end{array}$$(10)

In Sects. 3.3, 3.5, and 3.6 we study the parameters [a_s, N, σ_AL] and obtain the functions Φ₀, Φ₂, and Φ₃. In Sects. 3.7 and 3.8, we derive the functions Φ₄ and Φ₅ for parameters [e, i, ß]. We show the 1yr degeneracy in Sect. 3.4 and use the function Φ₁ to avoid this problem.

3.1 Bayesian model

For N mock Gaia observations, the chi-square χ² is defined by the equation $${\chi }^2={\displaystyle \sum _i^N{\left(\frac{\eta \left(t_i\right)-\eta {\left(t_i\right)}_{\text{true}}}{{\sigma }_{\text{AL}}}\right)}^2.}$$(11)

where η(t_i) is the predicted displacement, while η(t_i)_true is the input mock displacement and t_i is the time of the i-th observation. Then the Bayesian posterior probability is the same as Eq. (12) in A19, ln(posterior) ~ ln(prior) + ln(likelihood) $$\begin{array}{lll}\ln \left(\text{posterior}\right)\hfill & ~\hfill & \ln \left(\text{prior}\right)+\ln \left(\text{likelihood}\right)\hfill \\ \hfill & =\hfill & \ln \left(\text{prior}\right)-\frac{1}{2}\cdot {\chi }^2,\hfill \end{array}$$(12)

where ln(prior) is the prior of the parameters. A non-informative flat prior is used for [Δα₀ cos δ, Δδ, µ_α cos δ, µ_δ], while uniform distributions with upper and lower limits are used as the prior for the other parameters.

If we have additional N_RV mock RVs, then χ² will be $${\chi }^2={\displaystyle \sum _i^N{\left(\frac{\eta \left(t_i\right)-\eta {\left(t_i\right)}_{\text{true}}}{{\sigma }_{\text{AL}}}\right)}^2+{{\displaystyle \sum _j^{N_{\text{RV}}}\left(\frac{\text{RV}\left(t_j\right)-\text{RV}{\left(t_j\right)}_{\text{true}}}{{\sigma }_{\text{RV}}}\right)}}^2.}$$(13)

where t_j is the epoch of the RV. We focus on the astrometry data in this work and give a short discussion on the influence brought by additional RV data.

Fig. 2 Relative bias Δ(mp) at different P for the parameter combinations of [mp, ms, P, ϖ]. The colour of each point stands for a different distance D (see colour bar), which is 1 /ϖ. The size of the point gives the type of the primary, a 1.4 M⊙ NS or a 8 M⊙ BH. The grey horizontal line indicates the 0.1 limit, while the red dotted line is P = 50 d.

3.2 Model feasibility

In order to test the feasibility of the astrometry mass measurement method, we performed simulations for different parameter combinations without adding any observational errors in the position data. From O20 we collected some typical combinations for a BH/NS-LC binary systems. For the case of BH as the primary, m_p is set at 8 M_⊙, while m_s = [1, 8, 50] M_⊙. For the case of NS as the primary, m_p is set at 1.4 M_⊙ and m_s is changed to [1, 5, 10] M_⊙. The orbit period P is between 2 days and 4500 days. We used short chains here since we started from the true solution to test whether the solution can remain stable around it. For each input source we ran 30 chains for 6500 steps, with 1500 steps as the burn-in.

In Fig. 2 we plot Δ(m_p) = ∣m_p,pre – m_p,in∣/m_p,in, the relative bias of m_p. Here m_p,pre is the m_p predicted by simulation result, while mp,in is the input parameter m_p. The bias Δ(m_p) decreases with the input period. When the period is longer than 50 days, Δ(m_p) is less than 10% of the input m_p. For sources with shorter periods Δ(m_p) can exceed 30%, which means the solution might be useless. As a result, we set 50 days as the lower limit of the period of Gaia solvable sources, which is also consistent with the cadence of Gaia (Yamaguchi et al. 2018; Shikauchi et al. 2022). For the sources below this period limit, the observation data may contain more than one unique solution or lead to other incorrect solutions.

We derived the upper limit using the result from Lucy (2014), who studied the mass measurement from incomplete astrometric binary orbits. The author points out that the period of a solvable astrometry binary should not be longer than 2.5 times of the observation time in order to keep a small bias in mass measurement. We applied this criterion to a 5 yr Gaia mission, which indicates a 12.5 yr period upper limit.

Fig. 3 Main relation Φ0 constructed by a set of fitted lines. The simulation results are shown as points of different colours for different periods, from 50 days to 4500 days (see legend at right), while a series of fitted lines in the log – log space are also plotted in the same colour. The thick blue line is the fitting result of P = 1 yr. For comparison, the dotted line shows the result of A19, calculated by Eq. (9).

3.3 Φ₀; The main relation of a_s and P

We studied the mass combinations in Sect. 3.2, with P from 50 days to 4500 days, and thus the semi-major axis a_s varies from 0.009 AU to 9.9 AU. Except for [m_p, m_s, P, a_s], the other parameters [Δα₀ cos δ, Δ δ₀,µ_δ cos δ,µ_δ,ϖ, e, i, ω, Ω] are described in Sect. 2.4. The epoch of periastron passage T_p is fixed at 0.3P for different P. In the simulation of the main relation Φ₀, all of the parameters and their ranges are listed in Table 1. These upper and lower limits are able to cover almost the full range of all the MCMC points.

We performed the MCMC simulation six times for each parameter combination. Every time we generated a new set of mock data and ran 30 walkers starting from the neighbourhood around the true parameters. As was done in A19, each sampler walks 40 000 steps, which contains 10 000 steps for burn-in. The walker number and the chain length are also used in the following MCMC simulation of this work. We generated new mock data in each simulation in order to determine Gaia’s ability to observe a system rather than a specific data set, which helps us avoid the bias of extreme mock data. We only repeated six times at every point for the purpose of covering the target parameter space efficiently. This method is proved feasible in Sect. 3.10.

In Fig. 3 we used a series of straight lines to fit the Ξ–a_s relation with different periods. They have the following forms:

$${\log }_{10}\Xi =\text{k}\cdot {\log }_{10}{a}_s+b.$$(14)

Therefore, we obtained Φ₀(a_s, P) = Ξ, which is the main relation of the astrometric mass measurement ability of Gaia. Here, a_s is related to [P, m_p, m_s, ϖ] in the form of $$a_s=\varpi \cdot {\left[\frac{G{P}^2{m}_p^3}{4{\pi }^2{\left(m_p+m_s\right)}^2}\right]}^{\frac{1}{2}}.$$(15)

Combining the results in Fig. 4 and ignoring the special case around 1yr, we find that Ξ is indeed proportional to $a_s^{\text{k}}$ for P < 5 yr with an index k around −1. The intercept b is [−0.8, −0.6] for P < 1500 d, which is why A19 can use a single line to fit the data within 5 yr, while it changes greatly at longer period. In addition, we also find that Ξ around 1 yr shows a specific dispersion due to the degeneracy of the binary orbit and the annual motion. We discuss this phenomenon in Sect. 3.4.

There are several possible reasons why our main relation within 5 yr is different from the result in A19 (see dotted line in Fig. 3). First, A19 simulated the sources at [10, 100, 1000] pc while we simulated sources at [10, 1, 0.5, 0.25] mas corresponding to [100, 1000, 2000, 4000] pc. Second, we avoided using an exponential distance prior (Andrews et al. 2019; Bailer-Jones et al. 2018) by using the astrometry observable parallax directly instead of the distance. In order to assign the value of σ_AL, A19 used a σ_AL–m_G relation similar to Figure 9 in Lindegren et al. (2018) which varied with apparent magnitude. This σ_AL reaches 0.04–0.06 mas within 8–12.5 mag in the G band, while we used a constant 0.1 mas for all the σ_AL. We only used different σ_AL in Sect. 3.6, which studies the σ_AL itself. Therefore, our result reflects the true relation between Ξ and [a_s, P], and is easier to apply to different σ_AL determined by the luminosity of the secondary and the distance.

Fig. 4 Coefficients k and b of each line fitted in Fig. 3. The dashed grey lines show the periods of 1 yr and 5 yr. The outlier is the result at P = 1 yr, which was excluded from the study.

Fig. 5 Value of Ξ affected by the 1 yr degeneracy. The blue squares are the median values of the 50 simulations of each period. The blue shadow is the area between the 16th percentile and the 84th percentile. The orange line is calculated with Eq. (14), with k and b obtained from the red line in Fig. 4, ignoring the outlier of P = 1 yr.

3.4 Φ₁: Degeneracy around one year

In Fig. 3 the line at 1 yr shows an outlier of these fitting lines. The Ξ becomes quite large due to the coupling of the annual motion of the Earth and the 1 yr period orbit of the source. A similar phenomenon, shown by Holl (2011), is that the detectability of the astrometric planet with Gaia decreased at around 1 yr. In Fig. 5 we chose a system with [m_p, m_s, m] as [8 M_⊙, 8 M_⊙, 1 mas], and varied the period around 1 yr. We repeated each simulation ten times in order to get the uncertainty of Ξ at the same time.

In Fig. 5 there is a clear peak around 1 yr. The simulation result is close to the interpolation model at P = 320d and P = 420 d, which is calculated by Eq. (14) with k and b obtained from the red line in Fig. 4. To avoid the influence from the 1yr degeneracy, we just omitted the sources with periods within 320 d and 420 d by $${\text{Φ}}_1\left(P\right)=0,320d<P<420\text{\hspace{0.17em}d.}$$(16)

We did this for two reasons. First, Ξ around 1 yr is enlarged by several times due to the 1 yr degeneracy; only a few sources in this period range would satisfy our solvable criteria in Sect. 4.3. Second, we find that only 3% of the BH/NS-LC systems are P ∈ [320 d, 420 d] in our simulated sample in Sect. 4.3. Since the 1yr degeneracy has limited influence on our final BH/NS-LC population, we just omitted this period range.

3.5 Φ₂: Observation number N

Gaia shows a strong scanning pattern in observation numbers, astrometric errors, among others; therefore, the sources at different places usually have different results in observation numbers and scanning angle distributions (e.g. Lindegren et al. 2018, 2021). For example, we have 62 observations at the location of LB-1 from GOST, which is different from the observation information used in A19. Here we ran the simulations for all the parameter combinations in Sect. 3.3 with 124 observation epochs, which is twice the number of epochs at this position. Instead of the original epochs and angles in the GOST file, we added an additional epoch for each epoch with the same scanning angle. The additional epoch is sampled from a Gaussian distribution, N(0d, 20 d). We calculate the value of Ξ_{N = 124}/Ξ_{N = 62}, which has a median ratio ${0.70}_{-0.18}^{+0.29}$.

If the mission extended to ten years rather than five years, the longer observation time would also accumulate more observations. Thus, Lindegren et al. (2021) gives an expected improvement scale as T^−1/2 for the parallaxes, positions, and their uncertainty. Here we extrapolated this scale to our mass measurement $$\frac{{\text{Ξ}}_T}{{\text{Ξ}}_{T=T_0}}={\left(\frac{T}{T_0}\right)}^{-\frac{1}{2}},\left(T_0=5\text{yr}\right).$$(17)

If extending the observation mission were the same as increasing the observation number, Eq. (17) would be the same as the following equation: $${\text{Φ}}_2\left(N,N_0\right)={\left(\frac{N}{N_0}\right)}^{-\frac{1}{2}},\left(N_0=62\right).$$(18)

Our scale is ${0.70}_{-0.18}^{+0.29}$, very close to 1/ $\sqrt{2}$ or 0.707. Thus, we confirm the scale of Eq. (18). For simplicity, we only considered the variation caused by increasing the number of epochs for a 10 yr Gaia mission in this work. If we do a simulation up to 25 yr according to the analysis in Sect. 3.2, we would get a more realistic results for a 10 yr Gaia mission.

3.6 Φ₃: Gaia epoch data uncertainty σ_AL

The uncertainty of Gaia epoch data σ_AL is a direct variable that affects the mass measurement, which appears in Eq. (9): $\Xi ~\frac{a_s}{{\sigma }_{\text{AL}}}$. In the study of astrometric planets, Perryman et al. (2014) defined the astrometric signal-to-noise ratio, $S\text{/}N\equiv \frac{\alpha }{{\sigma }_{\text{FOV}}}=\frac{{\alpha }_s}{{\sigma }_{\text{AL}}}$ (the form in this work). We used a constant σ_AL in most of our simulations; therefore, we decided to study this parameter in an additional simulation test, which is equivalent to studying the effect of the S NR.

We studied a BH-LC system with [m_p, m_s] as [8 M_⊙, 8 M_⊙], where a_s is proportional to ϖ · $P^{\frac{2}{3}}$. So we chose different periods for sources at [0.5, 1, 10] mas, of which the SNR falls in the range [1,100]. In Fig. 6 we normalised the Ξ to Ξ(σ_AL = 0.1), and we got a quite good relation of Φ₃(σ_AL): $${\text{Φ}}_3\left({\sigma }_{\text{AL}}\right)=\frac{0.1}{{\sigma }_{\text{AL}}}.$$(19)

It should be noted that only the sources with ϖ > σ_AL and a_s > 1.5 σ_AL follow this relation. The sources with smaller ϖ or a_s would get higher Ξ.

Fig. 6 Values of Ξ normalised to ${\Xi }_{{\sigma }_{\text{AL}}=0.1}$. The blue points ares the median of all the values of $\Xi \text{/}{\Xi }_{{\sigma }_{\text{AL}}=0.1}$ with the same σAL. The red line is the Φ3 that equals $\frac{0.1}{{\sigma }_{\text{AL}}}$.

3.7 Φ₄: Inclination i and eccentricity e

In this section we describe the simulations we performed for different combinations of [i, e] when P is [500, 1825, 3650, 4500] d. For each period P, we fixed the [m_p, m_s, ϖ] at [8 M_⊙, 8 M_⊙, 1 mas]. The other parameters are the same as in Sect. 3.3, except for the inclination angle i and the eccentricity e. We simulated six times for every combination of [i, e] when i = [0°, 30°, 45°, 70°, 90^°] and e = [0.01, 0.5, 0.6, 0.75, 0.90].

We normalised the Ξ of each point to Ξ = 0.01, i = 30°) and find that it is hard to give a simple equation for Φ₄ (P, e, i). Therefore, we divided the period into four parts and use interpolation functions to describe Φ₄ (see Fig. 7).

There is a total trend in the four panels of Fig. 7 that Φ₄ increases when i gets larger. At high inclination, Φ₄ increases to at high eccentricity in panel a, while it has an opposite relation in panels c and d. In panel b, Φ₄ increases at both the higher and lower end of e, and decreases in the middle. The ratio of the maximum to the minimum of Φ₄ in the four panels could reach 4.5–8. Thus, Φ₄ would influence the final distribution of eccentricity, which should not be omitted in statistical studies.

3.8 Φ₅: Ecliptic latitude β

Considering in addition a [8 M_⊙, 8 M_⊙] system at 1 mas, we select nine positions for the simulations such that the absolute value of their ecliptic latitude ∣β∣ covers the range [0°, 85°]. The period falls in the range of[100d,4500d], while the other parameters are the same as in Sect. 3.3. In order to perform simulations at different positions, we used new observation data from the GOST program of each position.

In Fig. 8 we show the ratio Ξ/Ξ_{∣β∣ = 1.4°} versus ∣β∣. All of the ratio values have been corrected with observation numbers according to Sect. 3.5. The Ξ is nearly the same as Ξ_{∣β∣ = 1.4°} from ∣β∣ = 1.4° to ∣β∣ = 43°, while it drops to 0.7Ξ_{∣β∣=1.4°} over ∣β∣ = 62°. Here the point with β = 1.4° corresponds to the position [α₀, δ0]=[92.95°, 22.82°]. Thus, Φ₅(β) is apiecewise function: $${\text{Φ}}_5\left(\beta \right)=\{\begin{array}{l}\begin{array}{ll}1,\hfill & \left|\beta \right|\le 43°\hfill \end{array}\hfill \\ -0.0158\left|\beta \right|+1.679,43°<\left|\beta \right|\le 62°\hfill \\ \begin{array}{ll}0.7,\hfill & \left|\beta \right|\le 43°\hfill \end{array}\hfill \end{array}.$$(20)

3.9 Radial velocity

In this section we focus on a NS-LC system located at [0.5, 1, 10] mas which contains a 1.4 M_⊙ NS and a 1 M_⊙ companion. We used the same parameters in Sect. 3.3 with 17 additional RVs, and repeated the simulation six times for each P. The epoch of RV observation is selected from the GOST file randomly in each simulation for simplicity and universality. The RV uncertainty is a Gaussian error with σ_RV = 1 km s⁻¹.

In Fig. 9 it is clear that the additional RV data could improve the solution for sources at [0.5, 1] mas, while it helps little for sources at 10 mas. Considering Eq. (6), Ξ is proportional to $\frac{a_s^3}{{\varpi }^3{P}^2}$. Therefore, the RV data cannot help constrain the parameters better, such as [a_s, ϖ, P], if they have already reached the astrometry limit. Furthermore, the effect of RV data are not as significant at larger a_s as at smaller a_s. The reason is that the smaller a_s corresponds to a shorter period, which leads to a higher $\frac{K}{{\sigma }_{\text{RV}}}$. The greatest improvement happens at P = 1 yr, showing that the RV decouples the 1yr degeneracy of a_s and ϖ. After checking the uncertainty of many parameters or parameter combinations, we find that at P = 1 yr the RV helps constrain the uncertainty of orbit radius $R_s\left(\equiv \frac{a_s}{\varpi }\right)$ down to 5% of the same value calculated from the simulation without RV data. These phenomena suggest that RV is very useful to get a better binary solution.

3.10 Summary of the Ξ-relation

The final result of our simulation, the Ξ-relation, is a function containing six subfunctions. In summary, Φ₀(a_s, P) and Φ₁ (P) describe the overall profile of the full relation; Φ₂(N) and Φ₃ (σ_AL) give us the result of the variation in the Gaia observations; and Φ₄ (P e, i) and Φ₅(β) tell the influence of the orbit shape and the system location. The epoch observation error σ_AL is determined by the apparent magnitude m_G from Fig. A.1 in Lindegren et al. (2021). Therefore, we limit the mG within [6^m, 21^m]. In fact, most sources in Gaia would be fainter than 6^m. The sources brighter than 6^m might have worse epoch precision according to the same figure. The domain of our function is summarised in Table 2.

To show the ability of the Ξ-relation, we give an example with fixed [m_p, m_s, m_s, ϖ] at [8 M_⊙, 8 M_⊙, 1 mas] and varied period in the range of [100 d, 4500 d]. We repeated the simulation 25 times for each parameter combination to estimate the uncertainty of Ξ. Figure 10a not only gives the value of Ξ, but it also tells us the uncertainty of this error. Since we omitted the period around 1 yr, most uncertainty of Ξ is below 1%, which means that it is feasible for us to get the Ξ-relation with only six repeats of each simulated data point mentioned in Sect. 3.3.

In Fig. 10a we plot the predicted Ξ from our final relation. The Ξ-relation is almost consistent with the simulation except for a tiny difference at P < 100 d. Figure 10b shows the $\frac{\sigma \left(m_p\right)}{m_p}$, calculated from Ξ using Eq. (7), which is also consistent with the point obtained directly from the simulation. This suggests the correctness of Eq. (7), which is very useful in studying other astrophysical problems. Here we simply set σ (m_s) equal to 0. Compared to the uniform model in A19, our Ξ-relation extends to longer periods (up to 12.5 yr), and it is also more adaptable when the period is shorter than 5 yr.

Fig. 7 Values of Φ4 varied with i and e in different period ranges. For 50 d < P ≤ 1150 d a 2D nearest-neighbour interpolation function is used to fit Φ4. When P falls in the range of [1150 d,2700 d], [2700 d,4050 d] and [4050 d,12.5 yr], linear interpolation functions are used to fit Φ4.

Fig. 8 Φ5 fitted to the ratio of Ξ∣β∣ and ${\Xi }_{\left|{\beta }_0\right|}$. Each point is the median of the ratio for all the periods in [100 d,4500 d] with the same ∣β∣. The Ξ drops to 0.7 ${\Xi }_{\left|{\beta }_0\right|}$ when ∣β∣ is larger than 62°.

Fig. 9 Values of Ξ of a NS-dwarf system (1.4 M⊙ – 1 M⊙) with inclination angle i = 30°. Sources at [10, 1, 0.5] mas are plotted in red, yellow, and blue, respectively. As a comparison, the dashed line is plotted in the same colour as the Ξ of the same source with no RV data.

Fig. 10 One example of a simulated source and the Ξ-relation. (a) Median of Ξ (black points) calculated from the simulation result in Sect. 3.10 and the prediction of our model (red line) and the A19 model (blue dashed line). (b) Median of $\frac{\sigma \left(m_p\right)}{m_p}$ (black points) obtained from the simulation result directly, and the predicted $\frac{\sigma \left(m_p\right)}{m_p}$ by applying Eq. (7) to our predicted Ξ (red line) and that of A19 (blue dashed line). In both panels the grey vertical line indicates the 5 yr limit in the A19 model, and the grey horizontal dashed line shows the difference in Ξ and $\frac{\sigma \left(m_p\right)}{m_p}$. The interruption of our model around 1yr is determined by Φ1(P).

4 Application of the Ξ-relation

In this section we give two cases of application of the Ξ-relation. We introduce the method to obtain the apparent magnitude m_G at the beginning. Then in the first case we show the solvable area of some specific systems. In the second case we applied the Ξ-relation to a mock Galactic BH/NS-LC population to predict the number of Gaia solvable BH/NS-LCs. In both case we use Eq. (7) to calculate $\frac{\sigma \left(m_p\right)}{m_p}$ from Ξ, in which we assume σ(m_s), the error of m_s, is 0.1 m_s obtained by other method.

4.1 Apparent magnitude

To obtain σ_AL we used the apparent magnitude m_G to interpolate from Fig. A.1 in Lindegren et al. (2021). Given that the absolute magnitude of the secondary is M_G, the apparent magnitude m_G is determined by $$m_G=m_G+5\cdot {\log }_{10}\frac{D}{10\text{pc}}+A_G,$$(21)

where D is the distance to our Solar System and A_G is the interstellar extinction. We use a dust plane model to calculate this extinction (e.g. Chen et al. 2018; Nataf et al. 2013), $$\begin{array}{lll}{A}_v\hfill & =\hfill & {\rho }_0{\displaystyle {\int }_0^D\exp (-\frac{r\left|\sin \left(b\right)\right|}{H})\text{d}r}\hfill \\ \hfill & =\hfill & \{\begin{array}{ll}\rho \times D,\hfill & b=0°\hfill \\ \frac{{\rho }_{0}H}{\left|\sin \left(b\right)\right|}(1-exp(-\frac{D.\left|\sin \left(b\right)\right|}{H})),\hfill & b\ne 0°\hfill \end{array}\hfill \end{array}$$(22)

where b is the source’s Galactic latitude. The scale height of the Galactic dust plane model in Nataf et al. (2013) is H = 164 pc. This small scale height means that the dust concentrate around the Galactic plane, and the source with higher Galactic latitude would be less influenced by the extinction. The mean density of the dust plane, estimated by Chen et al. (2016) using 2000 nearby open clusters, is ρ₀ = 0.54 mag kpc⁻¹. We use the relative extinction value A_G/A_V = 0.789 (Wang & Chen 2019) to convert the V-band extinction A_V to A_G in Gaia’s G band. Finally, the σ_AL could be obtained from Fig. A.1 in Lindegren et al. (2021) with mG.

Fig. 11 Process to explore the Gaia solvable area for the dark-luminous systems in Sect. 4.2.

4.2 Gaia’s ability to solve systems with a dark companian

We test the ability of Gaia to measure the mass of the dark component in a binary as an example of using the Ξ-relation. This dark component could be a BH or NS, or even a WD or a planet. Here, we set $\frac{\sigma \left(m_p\right)}{m_p}$ <30% as the criterion for being astrometrically solvable.

Figure 11 gives the process to explore the solvable area for different systems. In Step 1 we set the system parameters; in particular, we varied the system period in 50 d ~ 12.5 yr and the parallax in 0.05 ~ 100 mas. In Step 2 the typical absolute magnitude in G band of the secondary are derived from the PARSEC model (Marigo et al. 2017). Then, the apparent magnitude m_G and σ_AL are obtained according to Sect. 4.1 in Step 3. After checking the solvable criteria in Step 4, it is able to give the solvable area for such a dark-luminous system. In this section the location is fixed at [α₀, δ₀], which is the location of LB-1. Although the Ξ-relation works well for nearby sources, we still set 10pc as the minimum distance for it. Since our model describes the motion of the source in a tangent plane, it might not hold for sources too close to the Solar System.

We calculated the value $\frac{\sigma \left(m_p\right)}{m_p}$, with Eq. (7), for different [m_p, m_s] combinations in Table 3. For a BH/NS/WD-LC system, if the luminous component is a 1 M_⊙ dwarf, the solvable area would be from 20pc to 2.06, 1.51, and 1.24 kpc for BHs, NSs, and WDs, respectively. On the other hand, if the luminous component is a 1 M_⊙ giant, this area would be 0.21 kpc to 6.44, 4.54, and 3.40kpc, respectively. The farthestlimitofthis solvable area is determined by the source with P ~ 1200 d, while the nearest limit is determined by the bright end of the σ_AL – m_G relation. For exoplanets we study a 0.003 M_⊙ planet, about 3.2 M_Jup, which orbits a 1 M_⊙ dwarf or a 0.6 M_⊙ dwarf. For consistency, we still use m_s for the luminous star and m_p for the planet. We find that Gaia is able to do the mass measurement of such a planet around the 0.6 M_⊙ dwarf in 68 pc with a precision better than 30%. The solvable area of a planet around 1 M_⊙ dwarf shrinks to only 31 pc because the astrometric wobble of the dwarf caused by the planet is smaller than the wobble of a 0.6 M_⊙ dwarf. We didn’t run the simulation for a giant with a planet, since the nearest limit of a giant with BH, NS, and WD is 220 pc, which is much farther than the upper limit of a planet-LC system.

In Fig. 12 we show the astrometrically solvable domain for different systems. These plots have the following features in common: (1) the mountain-like profile, determined by our Ξ-relation; (2) a gap around 1yr caused by the 1yr degeneracy; and (3) a clear lower distance limit for a bright secondary. In addition, taking panel 6 in Fig. 12 as an example, when the period is shorter than 1 yr, the precision is worse in the middle. Then it gets better at some distance, and finally goes to 30%. This feature is very clear in panels 4–7, which is caused by the non-monotonic σ_AL – m_G relation.

Fig. 12 Solvable domain for eight combinations of primary and secondary. The parameters and results are listed in Table 3. The luminous companion in panels [1,2,3,4,8] is a dwarf, and is a giant in panels [5,6,7]. The colour stands for the value of σ(mp)/mp, from 0 to 0.3. Different y limits are used in each panel, as indicated by ticks on the y-axis.

4.3 The Gaia solvable BH/NS-LC population

In order to examine the real ability of Gaia, we applied our Ξ-relation to a mock sample of a Galactic BH/NS-LC population using the binary evolution code MOBSE⁴, which is based on the code Binary Star Evolution (BSE) (Hurley et al. 2002), with some important upgrades to the models, including natal kicks, stellar winds, and supernovae (e.g. Giacobbo et al. 2018; Giacobbo & Mapelli 2018, 2019, 2020).

In MOBSE we adopted the following models. For the stellar winds we applied a metalicity-dependent model (see Belczynski et al. 2010; Chen et al. 2015). During the common envelope (CE) phase, the CE efficiency parameter a, which usually has a great uncertainty, is set to 3 for simplicity (e.g. Giacobbo et al. 2018; Giacobbo & Mapelli 2018; Shikauchi et al. 2022). Another CE parameter, the binding energy factor λ, is calculated by the process in Appendix A of Claeys et al. (2014). For the supernova (SN) model, we took the delayed SN engine model by Fryer et al. (2012). Finally, we used the model in Giacobbo & Mapelli (2020) to calculate the BH and NS natal kicks.

The whole process described in this section is shown in Fig. 13, which can be divided into two part. First, from Step 1 to Step 3 we generated the initial parameters of the binaries and construct the evolved sample of BH/NS-LCs (Sect. 4.3.1). Second, from Step 4 to Step 7 we scaled and put this evolved sample to the Milky Way, and applied the Ξ-relation (Sect. 4.3.2). We repeated the second part for 500 times and give the statistical result of the solvable BH/NS-LC sample (Sect. 4.5).

4.3.1 Evolved sample of BH/NS-LC

In Step 1 and Step 2 we evolved the binaries with MOBSE, which needs the input parameters, [m_p, m_s,P, e, Z, age]. The mass of the primary m_p is sampled from the power-law initial mass function (IMF) in Kroupa (2002). The mass of the secondary m_s, a less massive companion, is the product of m_p and mass ratio q, which is generated from a uniform distribution q ∈ [0.08/m_p, 1]. We referred to the study of Sana et al. (2012) to initialise the massive binary systems. We adopted a power-law period distribution log P∝ (log P)^−0.55 and a power-law eccentricity distribution e∝ e^−0.42. The eccentricity falls in the range [0, 0.9], while log P is in the range [0.15, 5.5], a larger range from Sana et al. (2014).

For the metallicity Z and the age we followed O20, who considered different Milky Way components, such as thin disk, thick disk, halo, and bulge. We ignored the globular clusters because their total mass only accounts for 0.005%–0.01% of the Milky Way stellar mass (e.g. Kravtsov & Gnedin 2005; Kruijssen & Mieske 2010). A short introduction of the distributions of Z and age are described here, based on O20.

Thin disk. The age of binaries in the thin disk is a uniform distribution in the range of [0, 10] Gyr; the metallicity decreases evenly from 1 Z_⊙ to 0.1 Z_⊙. The total mass of the thin disk is 90% of the mass of the whole stellar disk, which is 5.17 ± 1.11 × 10¹⁰ M_⊙ (Licquia & Newman 2015).

Thick disk. We followed O20, and considered that the thick disk formed uniformly between 9 Gyr and 11Gyr ago with a constant metallicity of 0.25 Z_⊙. The thick disk only accounts for 10% of the stellar disk.

Bulge. Referring to the figures in Kobayashi & Nakasato (2011), we assumed that the bulge can be devided into two parts. One part, containing 48% of the mass, formed from 12 Gyr to 10 Gyr ago with Z that increased from 0.1 Z_⊙ to 1 Z_⊙. The other part formed 10 Gyr ago, maintaining the metallicity at 1.5 Z_⊙.

Halo. Since most of the mass of the halo is dark matter, the stellar mass is only about 1.4 × 10⁹ Z_⊙ (out to 100 kpc) (Deason et al. 2019). As O20 did, we also set two specific metallicity values for the halo, Z = 0.01 Z_⊙, 0.02 Z_⊙. The binaries with 0.01 Z_⊙ formed between 12 Gyr and 11 Gyr ago, while the binaries with 0. 01 Z_⊙ formed between 11 Gyr and 10 Gyr ago.

To generate the evolved sample, which is used in Sect. 4.3.2, we generated 7.5 million to 17.5 million initial binaries for different Galactic components. We list the evolved results with the input initial systems in Table 4.

After the binary evolution, in Step 3, we added the absolute magnitude G to the luminous component in each system, where we chose the YBC database (Chen et al. 2019) for the bolometric correction. The method for obtaining G with the YBC website⁵ is described in Appendix A.

Fig. 13 Process of constructing the Gaia solvable BH/NS-LC population.

4.3.2 Galactic realisation

We used the following process in Step 4 to scale our evolved BH/NS-LCs to the whole Milky Way. According to O20, a sample of 2.5 million binaries with a primary from a 5–150 M_⊙ IMF accounts for a 3.1 × 10⁸ M_⊙ stellar population from a 0.08–150 M_⊙ IMF. Here the binary fraction is fixed at 50%. For a Galactic component c, which stands for [thin disk, thick disk, bulge,halo], we get N_sample,c, the number of BH/NS-LCs in component c, using the linear scaling relation $$\begin{array}{lll}{N}_{\text{binary},c}\hfill & =\hfill & 25\times {10}^6\cdot \frac{M_c}{3.1\times {10}^8{M}_{\odot }}\hfill \\ N_{\text{sample},c}\hfill & =\hfill & N_{\text{binary},c}\cdot \frac{N_{bh/ns-lc,s}}{N_{\text{ini},c}},\hfill \end{array}$$(23)

where M_c is the total mass of component c, and N_binary,c is the number of binary in component c, whose primary follows the 5–150 M_⊙ IMF; N_ini,c (Col. 4) and N_bh/ns–lc,c (Col. 2, Col. 3) are listed in Table 4. Here we used the same IMF to assign the mass for different Galactic components and do the sampling with replacements since we only evolved part of the potential binary systems in each component. We repeated the Galactic realisation 500 times to get a better estimation of the number of solvable sources and evaluate the uncertainty of the sampling.

In Step 5, we generated the position for each system using the following model. For the thin disk, thick disk, and bulge, we used a function MC_ amp in the Compact Object Synthesis and Monte Carlo Investigation Code (COSMIC) programme package (Breivik et al. 2020). For the MC_ amp we referred to McMillan (2011), who gives an axisymmetric model for the disks and bulge. For the halo stars, we used the dual stellar halo model in Kafle et al. (2014) with a break radius at 17.2 kpc.

Using the absolute magnitude from Step 3 and the Galactic position from Step 5, we calculated the apparent magnitude and get the σ_AL with Sect. 4.1 in Step 6. Then we applied our Ξ-relation and derived $\frac{\sigma \left(m_p\right)}{m_p}$ All the sources that satisfy our criteria in Step 7, $\frac{\sigma \left(m_p\right)}{m_p}$ <0.3, would be recorded in one realisation. The total sample of the 500 realisations makes up the sample of solvable BH/NS-LC population.

4.4 Galactic distribution of the solvable BH/NS-LCs

We select a random realisation with 48 BH-LCs and 102 NS-LCs, and plot them in Fig. 14 with a real Milky Way background⁶. It is clear that BH-LCs mainly fall on the disk⁷, while the NS-LCs disperse much more to all of the sky.

We also provide the solvable domain for the realistic simulated Galactic BH/NS-LC population in Fig. 15. In the top panel most of solvable sources fall in the m_G range of 8 to 14. In the bottom panel the Galactic solvable area is different for BH-LCs and NS-LCs. The NS-LCs and BH-LCs both mainly concentrate around the disk. The NS-LCs all fall in the area around the Solar System in 5 kpc, mostly in 3 kpc, while the BH-LCs distribute between the Solar System and the Galactic centre, clustering at the two ends. The BH-LCs could be observed in the dense Galactic centre at a similar frequency to the neighbourhood of Solar System because the BH-LCs in the centre usually have a larger a_s, a brighter secondary, and a higher number density.

Fig. 14 One Galactic realisation of our simulated BH/NS-LC population with a real Milky Way background. The blue points in the top panel are the 48 BH-LCs, while the red points in the bottom panel are the 102 NS-LCs. These are the typical numbers for Gaia solvable BH/NS-LCs in this work.

4.5 Discussion of the solvable BH/NS-LCs

In this work we use $\frac{\sigma \left(m_p\right)}{m_p}$ < 30% as the criteria for solvable H/NS-LCs, and ignore the possible obstruction in some dense areas. Table 5 gives the number of BH/NS-LCs that are used in the comparison in this section. The N1 column lists the predicted number of BH/NS-LCs with 30 or 10% precision in a 5 yr mission or a 10 yr mission, calculated by Eqs. (10) and (7). In a 5 yr Gaia mission, the number of solvable BH-LCs and NS-LCs are ${48}_{-7}^{+7}$ and ${102}_{-10}^{+11}$, respectively, of which 93% BH-LCs and 97% NS-LCs have periods P<5 yr. This result is consistent with the number range of 40–340 BH-LCs in Wiktorowicz et al. (2020). If the mission lasted for 10 yr, the number of solvable BH-LCs and NS-LCs would increase to 108 and 168. In both the 5 yr mission and the 10 yr mission, about 4% systems of the solvable BH-LCs have a precision $\frac{\sigma \left(m_p\right)}{m_p}<0.1$, while 7% of the solvable NS-LCs have this level of precision.

As a comparison, we calculated N₂ in the same way as N₁, but set the Φ₄(P, e, i) =1 in Eq. (10), which means N₂ is free from the influence of [e, i]. We find 66 BH-LCs and 168 NS-LCs in a 5 yr mission. Thus, our Φ₄(P, e, i) relation throws away 33% of the BH-LCs and 38% of the NS-LCs. This indicates that it is very important to do a complete simulation study of eccentricity and inclination before we can derive any statistical conclusion from the future observation result. For N3 we used Eq. (9), the result of A19, to get Ξ and calculated $\frac{\sigma \left(m_p\right)}{m_p}$. It is also free from different [e, i], so we compared it with N₂. In a 5 yr mission, N₃ predicts 2.3 times more BH/NS-LCs $\left(\frac{\sigma \left(m_p\right)}{m_p}<30\%\right)$, due to the difference between Eqs. (10) and (9). Here we ignored the sources with periods longer than 5 yr in N₂.

In addition, we applied the same criteria as A19 did, Ξ<30%, to the BH/NS-Giants in our evolved sample, and list the result in Table 6. We got 16±4 solvable BH-Giants and 25±5 solvable NS-Giants in our evolved sample, while A19 predicted 74 solvable BH-Giants and 190 NS-Giants. In the solvable sources of this work, only about 1/16 BH-Giants and 4/25 NS-Giants can reach Ξ<10%, while these fractions are much higher in A19, reaching 45/74 and 90/190, respectively.

We propose three reasons for this difference. The main reason might be that we used a different parameter distribution for the binary initialisation. Our values are different not only in the binary system parameters, such as period, eccentricity, mass ratio, and binary fraction, but also in metallicity and age. The second reason is that we used a Milky Way model combining different observation and simulation analyses, which are different from the simulated Milky Way-like galaxy m12i (Wetzel et al. 2016) used in A19. The third reason might be our different treatment of the bolometric correction and extinction, which have a great influence on the σ_AL. From the above comparison, there is a big difference in the number ratio of the precision at 30% and 10% between A19 and our comparison sample result. Thus, it seems that the third reason also plays an important role.

Fig. 15 Position and apparent magnitude distributions of BH-LCs and NS-LCs. Top: two-dimensional distribution of the solvable BH/NS-LCs in the mG – D plane. Bottom: galactic distribution of the solvable BH/NS-LCs in the $\sqrt{X^2+Y^2}-Z$ plane. The colour stands for the density of the solvable systems appearing in each bin.

5 Conclusion

We studied Gaia’s ability to derive the mass of the dark component in BH/NS-LC binary systems, using a realistic MCMC astrometry simulation. We extended the Ξ-relation from 5 yr to 12.5 yr, which is 2.5 times a 5 yr mission lifetime. A 10yr mission would improve the precision of mass measurement and obtain more BH/NS-LCs.

Assuming we obtained the mass m_s of the secondary, we connected the relative error of the primary m_p and real observable variables by an intermediate parameter KF, the Keplerian factor, in Eqs. (6) and (7).

Using a MCMC simulation method, we gave the Ξ-relation, which can be used to predict Ξ, the relative error of KF, with direct observation variables [a_s, σ_AL, N, P, e, i, β]. Assuming we obtained the mass of m_s by additional spectroscopic data, we were able to predict the $\frac{\sigma \left(m_p\right)}{m_p}$ with Eq. (7). The example in Sect. 3.10 shows that our Ξ-relation is more adaptable and precise than Eq. (9), the relation in A19. In addition, we especially explored the influence of the parameters [P, e, i, β], which could be used as the correction function for further statistical study of BH/NS-LCs with Gaia. The result in Sect. 4.5 shows that the eccentricity e and inclination i would influence about 38% of the solvable sources.

A Galactic BH/NS-LC population is generated by the evolution code MOBSE. In this process, we used a series of results from different observations or simulations. By applying our Ξ-relation to this BH/NS-LC population, we predict that ${48}_{-7}^{+7}$ BH-LCs and ${102}_{-10}^{+11}$ NS-LCs could be solved with a 5 yr Gaia mission, consistent with previous results (e.g. Wiktorowicz et al. 2020; Andrews et al. 2019). A comparison between the result of A19 and this work is given in Sect. 4.5, showing the influence of using different models of the Ξ-relation, binary population, Milky Way evolution and other details. In addition, we provide the distribution of these systems in the Milky Way as an instruction for future searching project.

Acknowledgement

The authors sincerely thank the anonymous referee for their useful comments. The authors are also grateful for the precious suggestion and help from Robert Soria. Y.W. is supported by National Science Foundation of China (NSFC) under grant numbers 11988101/11933004, National Key Research and Development Program of China (NKRDPC) under grant numbers 2019YFA0405504/2019YFA0405000, and Strategic Priority Program of the Chinese Academy of Sciences under grant number XDB41000000. S.L. is supported by the Youth Innovation Promotion Association CAS, the grants from the Natural Science Foundation of Shanghai through grant 21ZR1474100, and National Natural Science Foundation of China (NSFC) through grants 12173069, and 11703065. We acknowledge the science research grants from the China Manned Space Project with NO.CMS-CSST-2021-A12 and NO.CMS-CSST-2021-B10. A.O. is supported by the Polish National Science Center (NCN) grant Maestro (2018/30/A/ST9/00050). N.G. is supported by European Union’s H2020 ERC Starting grant No. 945155–GWmining, Leverhulme Trust grant No. RPG-2019-350, and Royal Society grant No. RGS-R2-202004.

Appendix A Absolute magnitude

In Sect. 4.3.1 we used the YBC database (Chen et al. 2019) to perform the bolometric correction and calculate the absolute magnitude of the luminous companion in our evolved BH/NS-LCs. To utilise the YBC database from the website, we sliced our evolved sample into different grids of [Mass, L, T_eff] (mass, luminosity, and effective temperature of the secondary), which contains 112,000 grids in total. We selected ten random companion stars in each grid, or we selected all of the stars if the grid contains fewer than ten stars. Some of the grids are empty if no star falls into them. This subsample, containing 21789 sources, is from the YBC website and we used the table for M_G. For each source i in the total sample, we selected a source j in the sub-sample that falls in the same grid of (Mass, L, T_eff) and has the closest T_eff to source i. We then calculated the absolute magnitude M_G,i as $$M_{G,i}=-2.5{\log }_{10}(\frac{L_i}{L_j})+M_{G,j},$$(A.1)

where the M_G,j is obtained from the YBC website.

References

All Tables

All Figures

Fig. 1 Schematic diagram of the sky-projected orbit of the luminous component and the Gaia epoch data. The central blue cross is the barycentre. The blue line shows the binary orbit of the luminous star. The black points are the mock Gaia epoch data, which show offsets from the orbit. These offsets are caused by the uncertainty in AL direction, indicated by the orientation of the error bar. All the points omit the proper motion and the annual motion.

In the text

	Fig. 2 Relative bias Δ(mp) at different P for the parameter combinations of [mp, ms, P, ϖ]. The colour of each point stands for a different distance D (see colour bar), which is 1 /ϖ. The size of the point gives the type of the primary, a 1.4 M⊙ NS or a 8 M⊙ BH. The grey horizontal line indicates the 0.1 limit, while the red dotted line is P = 50 d.
In the text

Fig. 3 Main relation Φ0 constructed by a set of fitted lines. The simulation results are shown as points of different colours for different periods, from 50 days to 4500 days (see legend at right), while a series of fitted lines in the log – log space are also plotted in the same colour. The thick blue line is the fitting result of P = 1 yr. For comparison, the dotted line shows the result of A19, calculated by Eq. (9).

In the text

	Fig. 4 Coefficients k and b of each line fitted in Fig. 3. The dashed grey lines show the periods of 1 yr and 5 yr. The outlier is the result at P = 1 yr, which was excluded from the study.
In the text

	Fig. 5 Value of Ξ affected by the 1 yr degeneracy. The blue squares are the median values of the 50 simulations of each period. The blue shadow is the area between the 16th percentile and the 84th percentile. The orange line is calculated with Eq. (14), with k and b obtained from the red line in Fig. 4, ignoring the outlier of P = 1 yr.
In the text

	Fig. 6 Values of Ξ normalised to ${\Xi }_{{\sigma }_{\text{AL}}=0.1}$. The blue points ares the median of all the values of $\Xi \text{/}{\Xi }_{{\sigma }_{\text{AL}}=0.1}$ with the same σAL. The red line is the Φ3 that equals $\frac{0.1}{{\sigma }_{\text{AL}}}$.
In the text

	Fig. 7 Values of Φ4 varied with i and e in different period ranges. For 50 d < P ≤ 1150 d a 2D nearest-neighbour interpolation function is used to fit Φ4. When P falls in the range of [1150 d,2700 d], [2700 d,4050 d] and [4050 d,12.5 yr], linear interpolation functions are used to fit Φ4.
In the text

	Fig. 8 Φ5 fitted to the ratio of Ξ∣β∣ and ${\Xi }_{\left|{\beta }_0\right|}$. Each point is the median of the ratio for all the periods in [100 d,4500 d] with the same ∣β∣. The Ξ drops to 0.7 ${\Xi }_{\left|{\beta }_0\right|}$ when ∣β∣ is larger than 62°.
In the text

	Fig. 9 Values of Ξ of a NS-dwarf system (1.4 M⊙ – 1 M⊙) with inclination angle i = 30°. Sources at [10, 1, 0.5] mas are plotted in red, yellow, and blue, respectively. As a comparison, the dashed line is plotted in the same colour as the Ξ of the same source with no RV data.
In the text

Fig. 10 One example of a simulated source and the Ξ-relation. (a) Median of Ξ (black points) calculated from the simulation result in Sect. 3.10 and the prediction of our model (red line) and the A19 model (blue dashed line). (b) Median of $\frac{\sigma \left(m_p\right)}{m_p}$ (black points) obtained from the simulation result directly, and the predicted $\frac{\sigma \left(m_p\right)}{m_p}$ by applying Eq. (7) to our predicted Ξ (red line) and that of A19 (blue dashed line). In both panels the grey vertical line indicates the 5 yr limit in the A19 model, and the grey horizontal dashed line shows the difference in Ξ and $\frac{\sigma \left(m_p\right)}{m_p}$. The interruption of our model around 1yr is determined by Φ1(P).

In the text

	Fig. 11 Process to explore the Gaia solvable area for the dark-luminous systems in Sect. 4.2.
In the text

	Fig. 12 Solvable domain for eight combinations of primary and secondary. The parameters and results are listed in Table 3. The luminous companion in panels [1,2,3,4,8] is a dwarf, and is a giant in panels [5,6,7]. The colour stands for the value of σ(mp)/mp, from 0 to 0.3. Different y limits are used in each panel, as indicated by ticks on the y-axis.
In the text

	Fig. 13 Process of constructing the Gaia solvable BH/NS-LC population.
In the text

	Fig. 14 One Galactic realisation of our simulated BH/NS-LC population with a real Milky Way background. The blue points in the top panel are the 48 BH-LCs, while the red points in the bottom panel are the 102 NS-LCs. These are the typical numbers for Gaia solvable BH/NS-LCs in this work.
In the text

	Fig. 15 Position and apparent magnitude distributions of BH-LCs and NS-LCs. Top: two-dimensional distribution of the solvable BH/NS-LCs in the mG – D plane. Bottom: galactic distribution of the solvable BH/NS-LCs in the $\sqrt{X^2+Y^2}-Z$ plane. The colour stands for the density of the solvable systems appearing in each bin.
In the text