Geometrical Distances of Extragalactic Binaries through Spectroastrometry

SA measures the variations of angular displacements of photocenters along the projected baseline with wavelengths, which can be derived from the surface brightness distributions of the object (R. Lachaume 2003; S. Rakshit et al. 2015). For a binary star system, the surface brightness can be modeled as

$$\begin{eqnarray}&&{ \mathcal O }({\boldsymbol{\alpha }},\lambda )=\displaystyle \sum _{i=1,2}\left[{{ \mathcal O }}_{{\rm{c}},i}({\boldsymbol{\alpha }}-{{\boldsymbol{\alpha }}}_{i})-{{ \mathcal O }}_{{\rm{\ell }},i}({\boldsymbol{\alpha }}-{{\boldsymbol{\alpha }}}_{i},\lambda )\right],\end{eqnarray} \tag{ 1 }$$

where α is the angular coordinates on the plane tangent to the celestial sphere, λ is the wavelength, and ${{ \mathcal O }}_{{\rm{c}},i}$ and ${{ \mathcal O }}_{{\rm{\ell }},i}$ are the surface brightness of the continuum and absorption line of the ith star, respectively. The angular displacement of the photocenter is therefore

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{\epsilon }}(\lambda ) & = & \frac{\displaystyle \int {\boldsymbol{\alpha }}{ \mathcal O }({\boldsymbol{\alpha }},\lambda ){{\rm{d}}}^{2}{\boldsymbol{\alpha }}}{\displaystyle \int { \mathcal O }({\boldsymbol{\alpha }},\lambda ){{\rm{d}}}^{2}{\boldsymbol{\alpha }}}\\ & = & \frac{{\sum }_{i=1,2}\left[{F}_{{\rm{c}},i}{{\boldsymbol{\epsilon }}}_{{\rm{c}},i}-{F}_{{\rm{\ell }},i}(\lambda ){{\boldsymbol{\epsilon }}}_{{\rm{\ell }},i}(\lambda )\right]}{{\sum }_{i=1,2}\left[{F}_{{\rm{c}},i}-{F}_{{\rm{\ell }},i}(\lambda )\right]}.\end{array}\end{eqnarray} \tag{ 2 }$$

Here, F_c,i and F_ℓ,i are the flux of the continuum and absorption line of the ith star:

$$\begin{eqnarray}\begin{array}{rcl}{F}_{{\rm{c}},i} & = & \displaystyle \int {{ \mathcal O }}_{{\rm{c}},i}({\boldsymbol{\alpha }}-{{\boldsymbol{\alpha }}}_{i}){{\rm{d}}}^{2}{\boldsymbol{\alpha }},\\ {F}_{{\rm{\ell }},i}(\lambda ) & = & \displaystyle \int {{ \mathcal O }}_{{\rm{\ell }},i}({\boldsymbol{\alpha }}-{{\boldsymbol{\alpha }}}_{i},\lambda ){{\rm{d}}}^{2}{\boldsymbol{\alpha }},\end{array}\end{eqnarray} \tag{ 3 }$$

and _c,i and _ℓ,i are the corresponding angular displacement of the photocenter of the continuum and absorption line:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\epsilon }}}_{{\rm{c}},i} & = & \displaystyle \int {\boldsymbol{\alpha }}{{ \mathcal O }}_{{\rm{c}},i}({\boldsymbol{\alpha }}-{{\boldsymbol{\alpha }}}_{i}){{\rm{d}}}^{2}{\boldsymbol{\alpha }},\\ {{\boldsymbol{\epsilon }}}_{{\rm{\ell }},i}(\lambda ) & = & \displaystyle \int {\boldsymbol{\alpha }}{{ \mathcal O }}_{{\rm{\ell }},i}({\boldsymbol{\alpha }}-{{\boldsymbol{\alpha }}}_{i},\lambda ){{\rm{d}}}^{2}{\boldsymbol{\alpha }}.\end{array}\end{eqnarray} \tag{ 4 }$$

For simplicity, the SA observation is conducted when there is no eclipse, and so _c,i = α_i. Defining the continuum flux ratio between the primary and secondary star as ℓ ≡ F_c,1/F_c,2 and the ratio between the line flux and total flux as

$$\begin{eqnarray}&&{f}_{{\rm{\ell }},i}(\lambda )\equiv \frac{{F}_{{\rm{\ell }},i}(\lambda )}{{\sum }_{i=1,2}\left[{F}_{{\rm{c}},i}-{F}_{{\rm{\ell }},i}(\lambda )\right]},\end{eqnarray} \tag{ 5 }$$

we have

$$\begin{eqnarray}&&{\boldsymbol{\epsilon }}(\lambda )=(1+{f}_{{\rm{\ell }},1}+{f}_{{\rm{\ell }},2}){{\boldsymbol{\epsilon }}}_{{\rm{c}}}-{f}_{{\rm{\ell }},1}{{\boldsymbol{\epsilon }}}_{{\rm{\ell }},1}-{f}_{{\rm{\ell }},2}{{\boldsymbol{\epsilon }}}_{{\rm{\ell }},2},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{\epsilon }}}_{{\rm{c}}}=\frac{\ell {{\boldsymbol{\alpha }}}_{1}+{{\boldsymbol{\alpha }}}_{2}}{\ell +1}\end{eqnarray} \tag{ 7 }$$

is the photocenter of the continuum flux. At the reference wavelength λ_r, where there is no absorption line, the displacement of the photocenter (λ_r) is only determined by that of the continuum photons:

$$\begin{eqnarray}&&{\boldsymbol{\epsilon }}({\lambda }_{{\rm{r}}})={{\boldsymbol{\epsilon }}}_{{\rm{c}}}.\end{eqnarray} \tag{ 8 }$$

The differential phase curve, which measures the phase difference between the observed wavelength λ and the reference wavelength λ_r, is therefore

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}\phi (\lambda ) & = & -\frac{2\pi }{\lambda }{\boldsymbol{B}}\cdot \left[{\boldsymbol{\epsilon }}(\lambda )-{\boldsymbol{\epsilon }}({\lambda }_{{\rm{r}}})\right]\\ & = & -\frac{2\pi }{\lambda }{\boldsymbol{B}}\cdot \left\{{f}_{{\rm{\ell }},1}\left[{{\boldsymbol{\epsilon }}}_{{\rm{c}}}-{{\boldsymbol{\epsilon }}}_{{\rm{\ell }},1}(\lambda )\right]+{f}_{{\rm{\ell }},2}\left[{{\boldsymbol{\epsilon }}}_{{\rm{c}}}-{{\boldsymbol{\epsilon }}}_{{\rm{\ell }},2}(\lambda )\right]\right\}.\end{array}\end{eqnarray} \tag{ 9 }$$

As shown by Equation (9), the differential phase Δϕ(λ) depends on the angular displacement of the photocenter of the absorption line _ℓ,i(λ) relative to the photocenter of the continuum _c, weighted by the relative intensity of the line f_ℓ,i. The continuum photocenter is the mean position of the two stars weighted by their luminosity. If the luminosity ratio ℓ between the primary and secondary stars is equal to their mass ratio q, the continuum photocenter _c will coincide with the mass center of the binary and keep stationary. Otherwise, it will orbit around the mass center. The profiles _ℓ,1(λ) and _ℓ,2(λ) will be redshifted and blueshifted respectively or vice versa, and the directions of _ℓ,1 and _ℓ,2 are opposite. As a result, the differential phase is an S-shaped curve, whose amplitude and width are mainly determined by the angular separation and orbital velocity of the binary.

For a star in an elliptical orbit with eccentricity e and period P, the eccentric anomaly E at time t can be obtained by solving the following Kepler equation (W. D. Heintz 1978):

$$\begin{eqnarray}&&2\pi \left(\frac{t-{T}_{0}}{P}\right)=E-e\sin E,\end{eqnarray} \tag{ 10 }$$

where T₀ is the time passage through the periastron.

The orbit of the binary and its projection onto the celestial plane are illustrated by Figure 1. To project the true orbit onto the celestial sphere, we define the following Thiele–Innes elements as

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}A & = & a(\cos \omega \cos {\rm{\Omega }}-\sin \omega \sin {\rm{\Omega }}\cos i),\\ B & = & a(\cos \omega \sin {\rm{\Omega }}+\sin \omega \cos {\rm{\Omega }}\cos i),\\ F & = & a(-\sin \omega \cos {\rm{\Omega }}-\cos \omega \sin {\rm{\Omega }}\cos i),\\ G & = & a(-\sin \omega \sin {\rm{\Omega }}+\cos \omega \cos {\rm{\Omega }}\cos i),\end{array}\end{array}\end{eqnarray} \tag{ 11 }$$

where a is the semimajor axis of the relative orbit of the secondary star relative to the primary star, i is the orbital inclination, Ω is the position angle of the ascending node, and ω is the argument of periastron. So the relative angular position of the secondary to the primary α₀ = (α_x, α_y) is given by

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{x} & = & \frac{A}{D}\left(\cos E-e\right)+\frac{F}{D}\sqrt{1-{e}^{2}}\sin E,\\ {\alpha }_{y} & = & \frac{B}{D}(\cos E-e)+\frac{G}{D}\sqrt{1-{e}^{2}}\sin E,\end{array}\end{eqnarray} \tag{ 12 }$$

where D is the distance of the binary and α_x and α_y point to the directions of increasing decl. and R.A., respectively. Putting the origin at the center of mass of the binary, the angular positions of the two stars are

$$\begin{eqnarray}&&{{\boldsymbol{\alpha }}}_{1}=-\frac{1}{1+q}{{\boldsymbol{\alpha }}}_{0},{{\boldsymbol{\alpha }}}_{2}=\frac{q}{1+q}{{\boldsymbol{\alpha }}}_{0},\end{eqnarray} \tag{ 13 }$$

where q is the mass ratio between the primary and secondary.

Zoom In Zoom Out Reset image size

To obtain the RVs of the two stars, we further define

$$\begin{eqnarray}&&C=a\sin \omega \sin i,H=a\cos \omega \sin i.\end{eqnarray} \tag{ 14 }$$

The relative RV of the secondary to the primary is given by

$$\begin{eqnarray}&&{V}_{0}=\frac{2\pi }{P\left(1-e\cos E\right)}\left(-C\sin E+H\sqrt{1-{e}^{2}}\cos E\right),\end{eqnarray} \tag{ 15 }$$

where v₀ > 0 means the star moves away from the observer. So the RVs of the two stars are

$$\begin{eqnarray}&&{V}_{1}=-\frac{1}{1+q}{V}_{0},{V}_{2}=\frac{q}{1+q}{V}_{0}.\end{eqnarray} \tag{ 16 }$$

These projected velocities measured by spectrograph are more sensitive to the mass ratio of the binary stars, in particular the RV ratio.

For simplicity, we model the surface of the star as a uniform disk where there is no eclipse, i.e.,

$$\begin{eqnarray}&&{{ \mathcal O }}_{{\rm{c}},i}({{\boldsymbol{\beta }}}_{i})=\left\{\begin{array}{l}{F}_{{\rm{c}},i}/\pi {\theta }_{i}^{2}\,{\rm{if}}\,| {{\boldsymbol{\beta }}}_{i}| \leqslant {\theta }_{i}\\ 0\,{\rm{if}}\,| {{\boldsymbol{\beta }}}_{i}| \gt {\theta }_{i},\end{array}\right.\end{eqnarray} \tag{ 17 }$$

where θ_i is the angular radius of the ith star and β_i = α − α_i is the relative angular displacement to the disk center.

For the absorption line, we include the rotation of the star and the thermal broadening, i.e.,

$$\begin{eqnarray}&&{{ \mathcal O }}_{{\rm{\ell }},i}({{\boldsymbol{\beta }}}_{i},\lambda )=\frac{{F}_{{\rm{\ell }},i}}{\pi {\theta }_{i}^{2}}\frac{1}{\sqrt{2\pi }{\sigma }_{i}}\exp -\frac{{[\lambda -{\lambda }_{i}({{\boldsymbol{\beta }}}_{i})]}^{2}}{2{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&{\lambda }_{i}({{\boldsymbol{\beta }}}_{i})={\lambda }_{0}\left[1+\frac{{V}_{i}}{c}+\frac{{v}_{{\rm{rot}},i}| {{\boldsymbol{\beta }}}_{i}| \sin ({\phi }_{{\rm{rot}},i}-{\phi }_{\alpha })}{c{\theta }_{i}}\right],\end{eqnarray} \tag{ 19 }$$

with σ_i being the thermal broadening of the absorption line, λ₀ being the central wavelength, c being the speed of light, v_rot,i being the projected rotational velocity, and ϕ_rot,i and ϕ_α being the position angles of the projected rotational axis and relative angular displacement β_i, respectively.

The eclipse happens when the sum of two stars' radii is larger than the distance between them, i.e., ∣α₀∣ < θ₁ + θ₂. When the secondary star is obscured, the flux from it will be

$$\begin{eqnarray}&&{F}_{{\rm{c}},2}^{\prime} ={\iint }_{| {{\boldsymbol{\beta }}}_{2}+{{\boldsymbol{\alpha }}}_{0}| \gt {\theta }_{1}}{{ \mathcal O }}_{{\rm{c}},2}({{\boldsymbol{\beta }}}_{2}){{\rm{d}}}^{2}{{\boldsymbol{\beta }}}_{2}.\end{eqnarray} \tag{ 20 }$$

For a uniform disk, we can figure out that

$$\begin{eqnarray}\begin{array}{rcl}{F}_{{\rm{c}},2}^{\prime} & = & {F}_{{\rm{c}},2}\left(1-\frac{{{\rm{\Theta }}}_{2}-\cos {{\rm{\Theta }}}_{2}\sin {{\rm{\Theta }}}_{2}}{\pi }\right.\\ & & \left.-\frac{{\theta }_{1}^{2}}{{\theta }_{2}^{2}}\cdot \frac{{{\rm{\Theta }}}_{1}-\cos {{\rm{\Theta }}}_{1}\sin {{\rm{\Theta }}}_{1}}{\pi }\right),\end{array}\end{eqnarray} \tag{ 21 }$$

when ∣α₀∣ > ∣θ₁ − θ₂∣, where

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Theta }}}_{1} & = & \arccos \frac{| {{\boldsymbol{\alpha }}}_{0}{| }^{2}+{\theta }_{1}^{2}-{\theta }_{2}^{2}}{2| {{\boldsymbol{\alpha }}}_{0}| {\theta }_{1}},\\ {{\rm{\Theta }}}_{2} & = & \arccos \frac{| {{\boldsymbol{\alpha }}}_{0}{| }^{2}+{\theta }_{2}^{2}-{\theta }_{1}^{2}}{2| {{\boldsymbol{\alpha }}}_{0}| {\theta }_{2}}.\end{array}\end{eqnarray} \tag{ 22 }$$

When ∣α₀∣ ≤ ∣θ₁ − θ₂∣, we have ${F}_{{\rm{c}},2}^{\prime} =0$ if θ₁ ≥ θ₂ and ${F}_{{\rm{c}},2}^{\prime} ={F}_{{\rm{c}},2}(1-{\theta }_{1}^{2}/{\theta }_{2}^{2})$ if θ₁ < θ₂. When the primary star is obscured, the calculation of the flux can be carried out similarly. If the star is described by quadratic or nonlinear limb darkening rather than a uniform disk, the exact analytic formulae for the eclipse can be found in K. Mandel & E. Agol (2002).

The model parameters (Θ) for calculating SA, RV, and LC are summarized in Table 1. Given SA, RV, and LC data (${\mathscr{D}}$), we can reconstruct the posterior probability distribution of the model parameters in the Bayesian framework:

$$\begin{eqnarray}&&P({\boldsymbol{\Theta }}| {\mathscr{D}})=\frac{P({\mathscr{D}}| {\boldsymbol{\Theta }})P({\boldsymbol{\Theta }})}{P({\mathscr{D}})},\end{eqnarray} \tag{ 23 }$$

where P(Θ) is the prior distribution of the model parameters, $P({\mathscr{D}})$ is a normalization factor, and $P({\mathscr{D}}| {\boldsymbol{\Theta }})$ is the likelihood for given data and parameters. We assume the prior probabilities of the model parameters are independent. Each parameter's prior range and distribution are presented in Table 1.

Table 1. Parameters Used in the Binary Star Model

Parameter	Unit	Implication	Fiducial Value	Range	Prior
Orbits
D	kpc	Distance	50	[1, 103]	Log-uniform
P	day	Orbital period	300	...	Fixed
T0	day	Time of periastron	0	...	Fixed
a	au	Semimajor axis	2	[0.1, 10]	Log-uniform
e	...	Eccentricity	0.3	[0, 0.95]	Uniform
i	deg	Orbital inclination	87.13	[0, 180]	$\mathrm{Cos}i$ uniform
Ω	deg	Position angle	45	[0, 360]	Uniform
ω	deg	Argument of periastron	30	[0, 360]	Uniform
q	...	Mass ratio	1.5	[1, 100]	Log-uniform
ℓ	...	Flux ratio	3	[0.1, 100]	Log-uniform
Stars
ri	a	Size	0.1, 0.08	[10−3, 1]	Log-uniform
vrot,i	km s–1	Projected rotational velocity	10, 5	[0.01, 100]	Log-uniform
ϕrot,i	deg	Position angle of rotational axis	120, 60	[0, 360]	Uniform
σi	km s–1	Thermal broadening	10, 10	[0.01, 100]	Log-uniform
EWi	Å	Equivalent width	1, 1	[0.1, 10]	Log-uniform

Note. Fiducial values follow the sequence of the first star and second star.

Download table as:  ASCIITypeset image

Suppose all measurements are independent and errors are Gaussian-distributed. The likelihood of the data is given by

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}P({\mathscr{D}}| {\boldsymbol{\Theta }}) & = & \displaystyle \prod _{i=1}^{{N}_{{\rm{t}},{\rm{sa}}}}{\prod }_{j=1}^{{N}_{\lambda }}{\prod }_{k=1}^{{N}_{{\rm{b}}}}\frac{1}{\sqrt{2\pi }{\sigma }_{\phi ,ijk}}\exp \left[-\frac{{({\phi }_{ijk}-{\phi }_{ijk}^{{\rm{mod}}})}^{2}}{2{\sigma }_{\phi ,ijk}}\right]\\ & & \times {\prod }_{i=1}^{{N}_{{\rm{t}},{\rm{sa}}}}{\prod }_{j=1}^{{N}_{\lambda }}\frac{1}{\sqrt{2\pi }{\sigma }_{{\rm{\ell }},ij}}\exp \left[-\frac{{({F}_{{\rm{\ell }},ij}-{F}_{{\rm{\ell }},ij}^{{\rm{mod}}})}^{2}}{2{\sigma }_{{\rm{\ell }},ij}}\right]\\ & & \times {\prod }_{i=1}^{{N}_{{\rm{t}},{\rm{rv}}}}{\prod }_{j=1}^{2}\frac{1}{\sqrt{2\pi }{\sigma }_{{\rm{v}},ij}}\exp \left[-\frac{{({V}_{ij}-{V}_{ij}^{{\rm{mod}}})}^{2}}{2{\sigma }_{{\rm{v}},ij}}\right]\\ & & \times {\prod }_{i=1}^{{N}_{{\rm{t}},{\rm{lc}}}}\frac{1}{\sqrt{2\pi }{\sigma }_{{\rm{c}},i}}\exp \left[-\frac{{({F}_{{\rm{c}},i}-{F}_{{\rm{c}},i}^{{\rm{mod}}})}^{2}}{2{\sigma }_{{\rm{c}},i}}\right],\end{array}\end{array}\end{eqnarray} \tag{ 24 }$$

where ϕ_ijk, ${\phi }_{ijk}^{{\rm{mod}}}$, and σ_ϕ,ijk are the measured value, predicted value, and measurement uncertainty of the differential phase at the ith time slot, jth wavelength channel, and kth baseline; F_ℓ,ij, ${F}_{\ell ,ij}^{{\rm{mod}}}$, and σ_ℓ,ij are the measured value, predicted value, and measurement uncertainty of the absorption line at the ith time slot and jth wavelength channel; V_ij, ${V}_{ij}^{{\rm{mod}}}$, and σ_v,ij are the measured value, predicted value, and measurement uncertainty of the RV of the jth star at the ith time slot; and F_c,i, ${F}_{{\rm{c}},i}^{{\rm{mod}}}$, and σ_c,i are the measured value, predicted value, and measurement uncertainty of the eclipsing LC at the ith time slot.

After establishing the formulations for the prior distribution of model parameters and the likelihood function, we can utilize appropriate Monte Carlo techniques (S. Sharma 2017) to sample the posterior probability distribution of the parameters. Specifically, we employ the Diffusive Nested Sampling algorithm (B. J. Brewer et al. 2011) for this purpose. Noteworthy for its ability to navigate intricate scenarios, including high-dimensional parameter spaces and multimodal distributions, this algorithm also facilitates the computation of model evidence essential for model contrasts and comparisons.

In order to study the impact of data quality and model parameters on distance measurement, we must construct a fiducial case as a benchmark for comparison. Table 1 lists the model parameters for simulating the fiducial data set. The resulting projected orbits of both stars are shown in the left panel of Figure 2.

Zoom In Zoom Out Reset image size

The SA observation is taken at a cadence of 41 days, lasting for 1 yr. Eight observations were generated after removing the observations at the time of the two eclipsing events. For simplicity, we keep the six projected baselines fixed in the UV plane during the observations, as shown in the right panel of Figure 2. The lengths of the baselines range from 100 m to 200 m. The spectral line for SA observation is Brγ, centered at 2.166μ m in the rest frame. The simulated line profiles and differential phase curves are all convolved by a Gaussian with FWHM = 30 km s^–1 (spectral resolution ${ \mathcal R }=10,000$) to account for instrument broadening. The bin size of the phase curve and spectra is half of the standard deviation of the Gaussian broadening function, i.e., 0.46 Å in this case. The uncertainties of the phase and flux measurements in each wavelength channel are 0.1∘ and 1%, respectively. The RV observation is taken at a cadence of 15 days, lasting for 1 yr. Two observations during the eclipse are also removed. The uncertainty of the velocity measurement is 2 km s^–1 for each star. The LC observation is taken at a cadence of 0.5 days, lasting for 1 yr, and the uncertainty of the flux measurement is 0.2%. In practice, we can always increase the cadence and precision of the observation by period folding.

We then apply Diffusive Nested Sampling to fit the binary model to the mock data and obtain the probability distribution of the model parameters. In practice, the orbital period and time passage through periastron can always be determined precisely from long-term monitoring of the binary and so will be fixed here, for simplicity. The simulated and reconstructed line profiles, differential phase curves, RVs, and LCs are shown in Figures 3, 4, and 5, respectively. As we can see, all mock data have been fitted very well within the margin of error. Posterior distributions of the model parameters are shown in Figure 6. First, the relative uncertainty of the distance is about 6.5% and the input value is within the 1σ range of the distribution. The degeneracy between the distance and other parameters is weak, and so the uncertainty can be further reduced by improving the quality of the data without introducing significant bias. Other orbital parameters are determined within 1% for dimensional parameters and 1∘ for angles. For star parameters, the relative sizes, thermal velocities, and equivalent widths are well determined, but the velocities and directions of rotation are poorly constrained, due to their degeneracies with thermal broadening under the current spectral resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After analyzing the fiducial case, we will vary the quality of different data to study their impact on the uncertainties of the model parameters, especially the uncertainty of the binary distance. The simulation results are summarized in Figure 7. As shown by the first row of Figure 7, the relative uncertainty of the distance (σ_D) is approximately proportional to the error of the differential phases (σ_ϕ) and the absorption-line flux (σ_ℓ), inversely proportional to the spectral resolution of the SA observation (${ \mathcal R }$), but less dependent on the error of the RV (σ_v), LC (σ_c), and the cadence of the LC (Δt_c). As σ_ϕ increases from 0.05∘ to 0.5∘, σ_D also rises from 4.6% to 12.1%. Similarly, σ_D changes from 4.2% to 15.3% when σ_ℓ varies from 0.5% to 5%. Conversely, as ${ \mathcal R }$ improves from 5000 to 20,000, σ_D decreases from 10.9% to 2.8%. From Equations (7) and (9), it follows that the amplitude of the differential phase depends on ℓ, which is determined by the analysis of eclipsing LCs. Approximately, we have

$$\begin{eqnarray}&&{\sigma }_{D}\unicode{x0007E}\sqrt{{\sigma }_{\ell }^{2}+{\sigma }_{\phi }^{2}+\,\rm{other\, terms}\,}.\end{eqnarray} \tag{ 25 }$$

In our fiducial scenario, we have σ_ℓ ~ 0.6% and σ_D ~ 6%. Consequently, the contribution of σ_ℓ to σ_D is minimal, indicating a lack of strong correlation between D and ℓ, as illustrated in Figure 6. This further elucidates why σ_D exhibits lesser dependence on the quality of the LC data.

Zoom In Zoom Out Reset image size

To further explore the error budget of the measurement of the distance, we also include parameters whose uncertainties have a similar dependence on data quality as the distance's uncertainties in Figure 7. The relative uncertainty of the semimajor axis (σ_D) has a weaker dependence on σ_ϕ and σ_ℓ but a stronger dependence on ${ \mathcal R }$, since the physical size is obtained through the integration of velocity information over time. Its contribution to σ_D can be seen in D ~ a/Δθ, where Δθ is the angular size of the binary measured by the differential phase. The uncertainties of the inclination angle (${\sigma }_{\cos i}$) and the position angle mainly depend on σ_ϕ and ${ \mathcal R }$. They affect the measurement of the distance by changing the projection of the binary orbit on the celestial sphere. The uncertainties of the mass ratio (q) increase slightly with larger σ_ϕ and σ_ℓ but decrease significantly with higher ${ \mathcal R }$. They can affect the distances of both stars to the mass center and their RVs. The uncertainties of the equivalent widths (EW₁ and EW₂) are very sensitive to σ_ℓ and ${ \mathcal R }$, and they can change the angular position of the photocenter by altering the weights from each star in Equation (9).

Now we will explore the impact of the model parameters' values on their uncertainties. Most dimensional parameters, such as D, a, and P, will only adjust the scale of the data curves, and so their effects are equivalent to that of the data quality. We will focus on those dimensionless parameters that can change the shapes of the data curves. The simulation results are displayed in Figure 8. Note that when we vary the relative size of the secondary star (r₂), we keep the size of the primary (r₁) fixed. But when we change r₁, we keep the ratio r₂/r₁ fixed.

Zoom In Zoom Out Reset image size

As shown by the first row of Figure 8, σ_D does not vary a lot as we alter the values of some model parameters and only increases slightly with smaller inclination (larger $\cos i$), higher light ratios (ℓ), larger r₂/r₁, and smaller r₁. Generally, it ranges from 4.3% to 9.3% for different values of these parameters.

For comparative analysis, we also illustrate in Figure 8 the correlation between the uncertainties of other parameters and their respective input values. It is observed that the uncertainties of $\cos i$, q, r₁, and r₂ exhibit constancy within a defined range of input values. Notably, the uncertainty in ℓ showcases a similar dependency on the input values of $\cos i$, ℓ, r₂/r₁, and r₁ as observed for D. A decrease in the inclination angle of the stellar size weakens the eclipsing effect, consequently leading to a reduction in the precision of the ℓ measurements, which in turn increases the uncertainty of the distance measurement.

Our geometric method for determining the distances of binaries is independent of empirical relations, eliminating bias from calibrations and sample selections. In this paper, we briefly outline the major sources of errors rather than accurately quantify them. The systematic uncertainties associated with this method primarily originate from the model bias related to the geometry of the components, including factors such as shapes and limb-darkening effects when fitting observational data. Limb darkening can produce complex yet regular brightness distributions that vary depending on stellar type, necessitating careful parameterization during the analysis of actual data. This phenomenon will also influence the stellar angular diameters derived from interferometric measurements. Normally, the mean disparity between a uniform disk and a limb-darkening disk's size is about 6% for  ~ 1000 stars listed in the JMMC Measured Stellar Diameters Catalog (G. Duvert 2016). Given that the stellar angular diameter comprises just 10% of the orbital semimajor axis for binaries with  ~ 1 au separations, the distance measurement error from this phenomenon should be below 0.6%, significantly lower than the statistical error. Nonetheless, as future observational precision advances, this phenomenon will gain prominence, necessitating its consideration.

On the other hand, the intrinsic spectral line profiles of individual stars may also play a crucial role in influencing the values of the orbital parameters derived from SA observations. Within our model, the individual line profile is intricately shaped by the combined effects of thermal broadening and the rotational characteristics of the star. Figure 6 illustrates that the thermal broadening velocities, projected rotational velocities, and the orientations of the rotation axes cannot be precisely constrained by the data, suggesting that the distance measurement is not highly reliant on the model employed for line profiles. However, when stars are sufficiently separated and the orbital velocity falls below the width of an individual star's line profile, the impact of the line shape becomes significantly pronounced. Additionally, the presence of extended stellar winds in massive stars introduces further complexity, as the inclusion of emission lines from these winds in SA observations may escalate the uncertainties in distance measurements (I. Waisberg et al. 2017).

To more accurately assess the magnitude of the systematic errors introduced by model selection, we propose applying the method to nearby binary star systems with known distances of various types, encompassing a range of characteristics, such as colors, shapes, metallicity, and other relevant attributes, for interferometric geometric measurements. Through this comparative analysis, we aim to identify the type of binary star that minimizes systematic errors, thereby optimizing measurements for future extragalactic binary distance determinations.

As demonstrated in Equation (9), the angular displacement of the photocenter is intricately linked to the luminosity ratio existing between the two stars in a binary system. Consequently, it is imperative to integrate eclipsing LCs into the data set for the precise determination of the luminosity ratio. During practical observations, if the stellar size is much smaller compared to the orbital semimajor axis or if the line of sight's inclination angle is more face-on, the weakening of the eclipsing effect can introduce heightened uncertainty in measuring the luminosity ratio. In such scenarios, the inclusion of the spectral energy distribution (SED) of the binary system becomes essential. Through the fitting and decomposition of these SEDs, the individual temperatures of the stars and their corresponding luminosity ratios can be derived (V. V. Jadhav et al. 2021; B. A. Thompson et al. 2021). While the absolute flux calibration of the SED holds lesser significance in this approach, a meticulous examination of the extinction across different wavelengths remains crucial for obtaining dependable measurements of the luminosity ratio.

In this subsection, we will discuss the feasibility of using existing interferometric equipment, particularly the next-generation beam-combiner GRAVITY+ on VLTI (GRAVITY Collaboration et al. 2022), to measure the distances of eclipsing binary systems in the LMC. For this purpose, we choose eclipsing binaries of the EA type in the LMC from the OGLE Collection of Variable Stars (A. Udalski et al. 2015), with I-band magnitudes brighter than 18 and orbital periods spanning from 100 to 1000 days. Subsequent cross-referencing with the GAIA database (Ó. Jiménez-Arranz et al. 2023) revealed a total of 265 binaries. Of these, 163 binaries possess K-band magnitudes obtained from the Two Micron All Sky Survey (M. F. Skrutskie et al. 2006).

Table 2 presents 17 binaries in the sample with RV observations (D. Graczyk et al. 2018; B. Hoyman et al. 2020). These systems are composed of two red giants, most of which are G- or early K-type giant stars. The K-band magnitudes of these systems range from 13 to 15, and their orbital parameters are consistent with those of our mock data, making them possible candidates for future interferometric observations.

We also show the colors and the G-band absolute magnitudes of the 163 selected binaries from the LMC in the left panel of Figure 9. 17 binaries with known orbital parameters are marked by black stars and located at the red giant branch of the H-R diagram. For the remaining targets, the red circles represent redder ones, with BP – RP > 0.7. They may have similar orbital parameters as those 17 targets and be appropriate targets for GRAVITY+ if they are bright enough. The blue circles are targets with BP – RP < 0.7. They are likely composed of O stars or B stars, and we need further spectroscopy data to justify whether they can be used for distance measurements. In the right panel of Figure 9, we present the distribution of the K-band magnitude for selected targets with BP – RP > 0.7. There are 22 targets with K < 14 and 74 targets with K < 15.

Zoom In Zoom Out Reset image size

The number of photons received by the interferometer per spectral channel for a target can be estimated by (H. Campins et al. 1985)

$$\begin{eqnarray}\begin{array}{rcl}{N}_{{\rm{phot}}} & \approx & 1.2\times 1{0}^{5-0.4(K-13)}{\left(\frac{D}{\,\rm{8.2}\,\,\,\rm{m}\,}\right)}^{2}\left(\frac{{\rm{\Delta }}t}{\,\rm{10}\,\,\,\rm{hr}\,}\right)\\ & & \times {\left(\frac{{ \mathcal R }}{1\times 1{0}^{4}}\right)}^{-1}\left(\frac{\eta }{\,\rm{1}\,\,\,\rm{\%}\,}\right),\end{array}\end{eqnarray} \tag{ 26 }$$

where K is the K-band magnitude of the target, D is the diameter of the aperture, Δt is the duration of the exposure, ${ \mathcal R }$ is the spectral resolution, and η is the overall throughput of the instrument. Under the best conditions, the measurement uncertainty of the phase at each spectral channel will be about (M. Shao et al. 1988)

$$\begin{eqnarray}&&{\sigma }_{\phi }={N}_{{\rm{phot}}}^{-1/2}\,\rm{rad}\,\approx 0.17{\left(\frac{{N}_{{\rm{phot}}}}{1.2\times 1{0}^{5}}\right)}^{-1/2}{\rm{\deg }}.\end{eqnarray} \tag{ 27 }$$

Ideally, the differential phase curve of the brightest target in Table 2 will have an amplitude of 1 − 2 ° and can be obtained using GRAVITY+ on VLTI with a signal-to-noise ratio (SNR) of 5 − 10. However, we acknowledge that additional noise sources, such as atmospheric pistons and mechanical vibrations, pose significant challenges to observations. As a first step, we suggest selecting long-period eclipsing binary systems with K ~ 10 and D ~ 5 kpc within the Milky Way to verify the purely geometric method for distance, as their theoretical SNRs for interferometric phases can exceed 100. Then we will decide whether to observe eclipsing binaries in the LMC using existing interferometers or to wait for future instruments with longer baselines, larger apertures, or higher throughputs, such as unprecedented interferometers combining Extremely Large Telescope and VLTI (~ 10 km baselines) in the next decade.