Astrometric mass measurement of compact companions in binary systems with Gaia

For binary systems with an unseen primary and a luminous secondary, the astrometric wobble of the secondary could be used to study the primary. With Gaia, it is possible to measure the mass of the black hole or neutron star with a luminous companion (hereafter BH/NS-LC). Our aim is to provide a method for predicting Gaia's ability in measuring the mass of BH/NS-LCs. We also tried to estimate the number of solvable BH/NS-LCs using Gaia. We used a realistic Markov chain Monte Carlo simulation of mock Gaia observations to obtain a relation between the uncertainty of mass measurement of the primary in BH/NS-LCs with the observable variables of the secondary astrometric orbit. Furthermore, we used the MOBSE code to evolve a Galactic BH/NS-LC sample with a combined Milky Way model. Our relation is applied to this sample to estimate the number of solvable BH/NS-LCs. We derived a good relation between the mass uncertainty and the binary parameters. For the first time, we show the quantitive influence of the period P, inclination i, eccentricity e, and ecliptic latitude $\beta$ to the mass measurement. Our results suggest that $48^{+7}_{-7}$ BH-LCs and $102^{+11}_{10}$ NS-LCs are solvable during a 5 yr Gaia mission. We also give the distribution of the distance and apparent magnitude of the Gaia solvable BH/NS-LCs. This solvable sample would be increased by additional spectroscopic data or a prolonged Gaia mission. The mass uncertainty relation could be used in future simulations of BH/NS-LCs observed by Gaia. The prediction of the solvable BH/NS-LCs is not only influenced by the process in generating the Galactic BH/NS-LC sample, but is also affected by our uncertainty relation. In particular, the relations of parameters such as $[P, e, i, \beta]$ are very useful to correct the selection effect in the statistic results of the future BH/NS-LC sample observed by Gaia.


Introduction
The first space astrometric mission, Hipparcos, observed about 120 thousand sources with an effective distance up to 1 kpc in 1989 . Of these sources, 18 thousand are recorded as non-single stars, 235 of which come with a binary orbit solution . 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. 2020).
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 Corresponding author: Jian Gao Corresponding author: Jifeng Liu for astrometric binaries from such astrometric solution by using Renormalised Unit Weight Error (RUWE) and excess_noise 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 astrometric 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 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 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 3,800-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 to 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;Chawla et al. 2021;Shikauchi et al. 2020Shikauchi et al. , 2022Breivik et al. 2019). 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 NS-Gs with a relative precision better than 0.3, which are somewhat solvable in Gaia's vision. This relative precision is defined by Equation 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 Section 5 we present our conclusions.

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.

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 et al. 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) 1 . In Fig. 1 we provide a schematic diagram of the mock Gaia observation of the luminous star in a fictitious NS-LC system 2 . We followed the instructions in Appendix A of Perryman et al. (2014) that describes the stellar motion in the AL direction 3 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: 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, [∆α 0 cos δ, ∆δ 0 , µ α cos δ, µ δ , ], to describe the system motion. 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 et al. 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 −1 . 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.

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: A = a s · (cos ω · cos Ω − sin ω · sin Ω · cos i), B = a s · (cos ω · sin Ω + sin ω · cos Ω · cos i), F = a s · (− sin ω · cos Ω − cos ω · sin Ω · cos i), (2) G = a s · (− sin ω · sin Ω + cos ω · cos Ω · cos i), 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 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 Then the RV can be calculated as where f is the true anomaly:

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: 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 = a s = m p m p + m s · (R p + R s ).
We then define a new parameter KF: In Eq. 6 the expression m 3 p M −2 tot is proportional to a 3 s 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 m p , so its relative error is influenced strongly by the mass ratio m s m p . Thus, using the error propagation method and Eq. 6, we derive the relative error of m p , 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 Ξ ≡ σ(KF)

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 .

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 [∆α 0 cos δ, ∆δ 0 , ω, Ω], which are fixed at [0.1 mas, 0.2 mas, 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 −1 , respectively, rather than the uniform 10 mas · yr −1 used in A19. The different components of proper motion [µ α cos δ, µ δ ] are [µ· sin π 6 , µ· cos π 6 ]. In order to use Eq. 1 to generate the AL direction data, we use GOST to get the observation information at [α 0 , δ 0 ]= [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−10 mas of σ AL (Lindegren et al. 2020), 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, 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 −1 .

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 Ξ, where P ≤ 5 yr and N 0 = 75. Here N is the number of Gaia observations, while N 0 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, β]: In Sects. 3.3, 3.5, and 3.6 we study the parameters [a s , N, σ AL ] and obtain the functions Φ 0 , Φ 2 , and Φ 3 . In Sect. 3.7 and Sect. 3.8, we derive the functions Φ 4 and Φ 5 for parameters [e, i, β]. We show the 1yr degeneracy in Sect. 3.4 and use the function Φ 1 to avoid this problem.

Bayesian model
For N mock Gaia observations, the chi-square χ 2 is defined by the equation 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 Equation 12 in A19, where ln(prior) is the prior of the parameters. A non-informative flat prior is used for [∆α 0 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 χ 2 will be 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.

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 m p,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.
3.3. Φ 0 : The main relation of a s and P We studied the mass combinations in Sect. 3.2, with P from 50 days to 4,500 days, and thus the semi-major axis a s varies from 10 −1 10 0 10 1 10 2 a s (mas) 0.02 mas 4 mas a s 10 −3 ∼100 mas 10 −4 mas 10a s mas P 50d∼12.5 yr 0 10P for different P. In the simulation of the main relation Φ 0 , 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 Ξ = k · log 10 a s + b.
Therefore, we obtained Φ 0 (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 Combining the results in Fig. 4 and ignoring the special case around 1yr, we find that Ξ is indeed proportional to a k s 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 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 Article number, page 5 of 15 A&A proofs: manuscript no. wyl2021a 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. 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 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 , ] 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=320 d 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 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. Φ 2 : 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. 2018Lindegren et al. , 2020. 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(0 d, 20 d). We calculate the value of Ξ N=124 /Ξ N=62 , which has a median ratio 0.70 +0.29 −0.18 . If the mission extended to ten years rather than five years, the longer observation time would also accumulate more observations. Thus, Lindegren et al. (2020) 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 Ξ: If extending the observation mission were the same as increasing the observation number, Eq. 17 would be the same as the following equation: Our scale is 0.70 +0.29 −0.18 , very close to 1/ √ 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. The uncertainty of Gaia epoch data σ AL is a direct variable that affects the mass measurement, which appears in Eq. 9: Ξ∼ a s σ AL . In the study of astrometric planets, Perryman et al. (2014) defined the astrometric signal-to-noise ratio, S /N≡ α σ FOV = a s σ 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 2 3 . So we chose different periods for sources at [0.5, 1, 10] mas, of which the S NR falls in the range [1,100]. In Fig. 6 we normalised the Ξ to Ξ(σ AL =0.1), and we got a quite good relation of Φ 3 (σ AL ): 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 Ξ.  We normalised the Ξ of each point to Ξ(e = 0.01, i = 30 • ) and find that it is hard to give a simple equation for Φ 4 (P, e, i). Therefore, we divided the period into four parts and use interpolation functions to describe Φ 4 (see Fig. 7).
There is a total trend in the four panels of Fig. 7 that Φ 4 increases when i gets larger. At high inclination, Φ 4 increases to at high eccentricity in panel a, while it has an opposite relation in panels c and d. In panel b, Φ 4 increases at both the higher and lower end of e, and decreases in the middle. The ratio of the maximum to the minimum of Φ 4 in the four panels could reach 4.5 to 8. Thus, Φ 4 would influence the final distribution of eccentricity, which should not be omitted in statistical studies. In Fig. 8

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 −1 . 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 a 3 s 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 K σ 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 (≡ a s ) 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.

Summary of the Ξ-relation
The final result of our simulation, the Ξ-relation, is a function containing six subfunctions. In summary, Φ 0 (a s , P) and Φ 1 (P) describe the overall profile of the full relation; Φ 2 (N) and Φ 3 (σ AL ) give us the result of the variation in the Gaia observations; and Φ 4 (P, e, i) and Φ 5 (β) 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. (2020). Therefore, we limit the m G 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. Article number, page 7 of 15 A&A proofs: manuscript no. wyl2021a    Figure 10.a not only gives the value of Ξ, but it also tells us the uncertainty of this error. Since we omitted the period around feasible for us to get the Ξ-relation with only six repeats of each simulated data point mentioned in Sect. 3.3. The interruption of our model around 1yr is determined by Φ 1 (P).
In Fig. 10.a 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 10.b shows the σ(m p ) 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.

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 σ(m p ) m p from Ξ, in which we assume σ(m s ), the error of m s , is 0.1m s obtained by other method.

Apparent magnitude
To obtain σ AL we used the apparent magnitude m G to interpolate from Fig. A.1 in Lindegren et al. (2020). Given that the absolute magnitude of the secondary is M G , the apparent magnitude m G is determined by m G = M G + 5 · log 10 D 10 pc 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), 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 = 0.54 mag kpc −1 . 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. (2020) with m G .

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 σ(m p ) m p <30% as the criterion for being astrometrically solvable. Fig. 11. Process to explore the Gaia solvable area for the dark-luminous systems in Sect. 4.2. 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 [α 0 , δ 0 ], which is the location of LB-1. Although the Ξrelation works well for nearby sources, we still set 10 pc 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.  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 σ(m p )/m p , from 0 to 0.3. Different y limits are used in each panel, as indicated by ticks on the y-axis.   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.40 kpc, respectively. The farthest limit of this 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.

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 4 , 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. , 2019. In MOBSE we adopted the following models. For the stellar winds we applied a metalicity-dependent model (see Belczynski et al. 2010 andChen et al. 2015). During the common envelope (CE) phase, the CE efficiency parameter α, which usually has a great uncertainty, is set to 3 for simplicity (e.g. 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).

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 powerlaw 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 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 9 M (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  for the bolometric correction. The method for obtaining G with the YBC website 5 is described in Appendix A.

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 8 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 N binary,c = 2.5 × 10 6 · M c 3.1 × 10 8 M N sample,c = N binary,c · N bh/ns−lc,c N ini,c , 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_samp in the Compact Object Synthesis and Monte Carlo Investigation Code (COSMIC) programme package A&A proofs: manuscript no. wyl2021a (Breivik et al. 2020). For the MC_samp 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 σ(m p ) m p . All the sources that satisfy our criteria in Step 7, σ(m p ) 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. 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 6 . It is clear that BH-LCs mainly fall on the disk 7 , while the NS-LCs disperse much more to all of the sky.

Galactic distribution of the solvable BH/NS-LCs
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.

Discussion of the solvable BH/NS-LCs
In this work we use σ(m p ) m p <30% as the criteria for solvable BH/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 N 1 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 Eq. 10 and Eq. 7. In a 5 yr Gaia mission, the number of solvable BH-LCs and NS-LCs are 48 +7 −7 and 102 +11 −10 , 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 σ(m p ) m p <0.1, while 7% of the solvable NS-LCs have this level of precision.
As a comparison, we calculated N 2 in the same way as N 1 , but set the Φ 4 (P, e, i)=1 in Eq. 10, which means N 2 is free from the influence of [e, i]. We find 66 BH-LCs and 168 NS-LCs in a 5 yr mission. Thus, our Φ 4 (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 N 3 we used Eq. 9, the result of A19, to get Ξ and calculated σ(m p ) m p . It is also free from different [e, i], so we compared it with N 2 . In a 5 yr mission, N 3 predicts 2.3 times more BH/NS-LCs ( σ(m p ) m p <30%), due to the difference between Eq. 10 and Eq. 9. Here we ignored the sources with periods longer than 5 yr in N 2 .
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.   Table 5. Predicted number of Galactic solvable BH/NS-LCs with Gaia. The statistical result is discussed in Sect. 4.5. Columns 1-3 list the basic information of the sample, from which we calculate [N 1 , N 2 , N 3 ] with three different methods. N 1 is our result, while N 2 and N 3 are for comparison. All three methods use Eq. 7 to obtain σ(m p )/m p , but different Ξ values. For N 1 , the Ξ is calculated from Eq. 10, while we set the Φ 4 (P, e, i) = 1 when we calculate the Ξ for N 2 . For N 3 , we use Eq. 9 to get the Ξ, and thus the binary period and the observation time are limited within 5 yr.

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 10 yr 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 Eq. 6 and Eq. 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