7.2.5 Orbital solutions
Calculation of orbital solutions
The orbit of a binary system is generally described by seven parameters: the first three are the period, $P$, the epoch of periastron passage, ${T}_{0}$, and the eccentricity, $e$. Then come the semimajor axis of the orbit described by the photocentre, ${a}_{0}$, the inclination of the orbital plane with respect to the sky, $i$, the position angle of the ascending node, $\mathrm{\Omega}$, and the periastron longitude measured from the ascending node, $\omega $. The precise definition of these parameters can be found in any astronomical handbook dealing with double stars (see, e.g. Binnendijk 1960; Heintz 1978). These four parameters are called Campbell’s elements, and each measurement of the position of the photocentre relative to the barycentre of the system is related to these elements by a completely nonlinear relationship.
The calculation of the orbital solutions is made partly linear by replacing the Campbell elements by the ThieleInnes elements. These, denoted by the letters $A$, $B$, $F$ and $G$, are related to the Campbell elements by the following relations:
$$\begin{array}{cc}A=\hfill & {a}_{0}(\mathrm{cos}\omega \mathrm{cos}\mathrm{\Omega}\mathrm{sin}\omega \mathrm{sin}\mathrm{\Omega}\mathrm{cos}i)\hfill \\ B=\hfill & {a}_{0}(\mathrm{cos}\omega \mathrm{sin}\mathrm{\Omega}+\mathrm{sin}\omega \mathrm{cos}\mathrm{\Omega}\mathrm{cos}i)\hfill \\ F=\hfill & {a}_{0}(\mathrm{sin}\omega \mathrm{cos}\mathrm{\Omega}+\mathrm{cos}\omega \mathrm{sin}\mathrm{\Omega}\mathrm{cos}i)\hfill \\ G=\hfill & {a}_{0}(\mathrm{sin}\omega \mathrm{sin}\mathrm{\Omega}\mathrm{cos}\omega \mathrm{cos}\mathrm{\Omega}\mathrm{cos}i)\hfill \end{array}$$  (7.12) 
and the alongscan abscissa of the photocentre is then given by the equation:
$$\begin{array}{cc}w=\hfill & {w}_{\mathrm{ss}}+\frac{\partial w}{\partial d}(\mathrm{cos}Ee)\times A+\frac{\partial w}{\partial a}(\mathrm{cos}Ee)\times B\hfill \\ & \\ & +\frac{\partial w}{\partial d}\sqrt{1{e}^{2}}\mathrm{sin}E\times F+\frac{\partial w}{\partial a}\sqrt{1{e}^{2}}\mathrm{sin}E\times G\hfill \end{array}$$  (7.13) 
where ${w}_{\mathrm{ss}}$, $\frac{\partial w}{\partial a}$ and $\frac{\partial w}{\partial d}$ are the same as is Section 7.2.4; $E$ is the eccentric anomaly, which depends on the epoch $t$ by the relation:
$$Ee\mathrm{sin}E=2\pi \frac{t{T}_{0}}{P}$$  (7.14) 
As ${w}_{\mathrm{ss}}$ depends on five unknowns, the orbital solution finally has twelve unknowns, nine of which enter into a linear relationship. The solution is therefore sought by exploring the space defined by $P$, ${T}_{0}$ and $e$ with the LevenbergMarquardt method. For each star, the period interval at the beginning of the calculation goes from 10 days to the duration covered by the observations divided by 0.6.
Significance of orbital solutions
The significance of the orbital solutions is defined as ${a}_{0}/{\sigma}_{{a}_{0}}$, since a semimajor axis as small as its uncertainty corresponds to an orbit that may well not exist. However, ${a}_{0}$ and its uncertainty must be deduced from a solution that has been calculated in ThieleInnes elements. For that purpose, the formulae provided by Binnendijk (1960) has been used:
$$\begin{array}{c}{a}_{0}=\sqrt{u+\sqrt{{u}^{2}{v}^{2}}}\hfill \\ \\ \mathrm{with}\{\begin{array}{cc}u\hfill & =\frac{1}{2}\left({A}^{2}+{B}^{2}+{F}^{2}+{G}^{2}\right)\hfill \\ v\hfill & =AGBF\hfill \end{array}\hfill \end{array}$$  (7.15) 
The uncertainties of the Campbell’s elements may be derived from $A$, $B$, $F$, $G$ and from the variance–covariance matrix (Halbwachs et al. 2023); for the semimajor axis, one obtains:
$$\begin{array}{c}{\sigma}_{{a}_{0}}=\frac{1}{2{a}_{0}}\sqrt{{\alpha}^{2}{\sigma}_{A}^{2}+{\beta}^{2}{\sigma}_{B}^{2}+{\varphi}^{2}{\sigma}_{F}^{2}+{\gamma}^{2}{\sigma}_{G}^{2}+2\alpha \beta {\mathrm{cov}}_{AB}+2\alpha \varphi {\mathrm{cov}}_{AF}+..+2\varphi \gamma {\mathrm{cov}}_{FG}}\hfill \\ \\ \mathrm{with}\{\begin{array}{cc}\alpha \hfill & =A\left(1+\frac{u}{\sqrt{{u}^{2}{v}^{2}}}\right)G\frac{v}{\sqrt{{u}^{2}{v}^{2}}}\hfill \\ \beta \hfill & =B\left(1+\frac{u}{\sqrt{{u}^{2}{v}^{2}}}\right)+F\frac{v}{\sqrt{{u}^{2}{v}^{2}}}\hfill \\ \varphi \hfill & =F\left(1+\frac{u}{\sqrt{{u}^{2}{v}^{2}}}\right)+B\frac{v}{\sqrt{{u}^{2}{v}^{2}}}\hfill \\ \gamma \hfill & =G\left(1+\frac{u}{\sqrt{{u}^{2}{v}^{2}}}\right)A\frac{v}{\sqrt{{u}^{2}{v}^{2}}}\hfill \end{array}\hfill \end{array}$$  (7.16) 
where ${\mathrm{cov}}_{AB}$ is the term of the variance–covariance matrix corresponding to the ($A$, $B$) pair, ie ${\mathrm{cov}}_{AB}={\sigma}_{A}{\sigma}_{B}{\rho}_{AB}$. The significance of the orbital solution is then estimated.
From small eccentricity orbits to pseudocircular orbits
Orbits with small eccentricities are a problem because the uncertainties of the ThieleInnes elements are much larger than the elements themselves. In addition, the uncertainty on ${T}_{0}$ can greatly exceed the value of the period. These defects come in fact from the indeterminacy of the periastron, which affects ${T}_{0}$, but also $\omega $. However, if the uncertainties of the Campbell elements are derived from the solution obtained with ThieleInnes elements, one may note that the problem is primarily one of presentation, since ${a}_{0}$, $i$ and $\mathrm{\Omega}$ actually have very acceptable uncertainties. Nevertheless, it was decided to transform the orbits with eccentricity $$, not by circularising them, but only by fixing the eccentricity and by shifting the periastron to the ascending node (obviously, such an operation is only possible because the uncertainty of the eccentricity is always much greater than 0.0005). Orbits transformed in this way are called pseudocircular orbits hereafter. This transformation was achieved with several steps:

•
The Campbell elements were derived from the ThieleInnes elements, and the ThieleInnes elements were computed again, assuming $\omega =0$ in Equation 7.12.

•
After this operation, it is easy to see that $G=AF/B$ and $F=BG/A$. Therefore, only 3 of the 4 ThieleInnes elements are considered as unknowns. $G$ has been arbitrarily discarded, checking first that $B$ is always larger than ${10}^{5}$; otherwise, $F$ would have been chosen instead of $G$.

•
The periastron epoch has been shifted by subtracting the time $P\times \omega /(2\pi )$ so as to correspond to the passage at the ascending node.

•
The variancecovariance matrix has been recalculated for the remaining ten unknowns, since $G$ was discarded and $e$ is considered as fixed to its initial value.

•
The resulting $10\times 10$ matrix is transcribed into a $12\times 12$ matrix, with the terms corresponding to the discarded unknowns set to zero.
It should be noted that, in reality, the uncertainties and correlations of $e$ and $G$ are not zero nor NaN, although they appear so in the delivered solution. It is up to the users to identify the solutions treated in this way (identified by the condition eccentricity_error is NaN), and to use them appropriately. For example, the uncertainty on $G$, ${\sigma}_{G}$ (resp. $F$, ${\sigma}_{F}$) may be estimated from the propagation of the uncertainties of the three other elements, thanks to the equation:
$${\sigma}_{G}=G\sqrt{{\left(\frac{{\sigma}_{A}}{A}\right)}^{2}+{\left(\frac{{\sigma}_{B}}{B}\right)}^{2}+{\left(\frac{{\sigma}_{F}}{F}\right)}^{2}+2\left(\frac{{\rho}_{AB}{\sigma}_{A}{\sigma}_{B}}{AB}+\frac{{\rho}_{AF}{\sigma}_{A}{\sigma}_{F}}{AF}+\frac{{\rho}_{BF}{\sigma}_{B}{\sigma}_{F}}{BF}\right)}$$  (7.17) 
A user wishing to calculate the uncertainties of the Campbell elements using the equations of Halbwachs et al. (2023) will still need the correlation coefficients relating $G$ to the other ThieleInnes elements. Since $G=AF/B$, one has:
$$\begin{array}{cc}\hfill {\rho}_{AG}=& +1\hfill \\ \hfill {\rho}_{BG}=& 1\hfill \\ \hfill {\rho}_{FG}=& +1\hfill \end{array}$$  (7.18) 
These recommendations also apply to the few calculated circular orbits (there are only 8 at the end of the calculation), which are presented in the same way.
Initial selection of orbital solutions
As with the acceleration solutions, a dynamic criterion was used to ensure the plausibility of the solution. We considered the astrometric mass function, defined as:
$${f}_{\mathcal{M}}=\frac{{a}_{0}^{3}\times {365.25}^{2}}{{P}^{2}{\varpi}^{3}}$$  (7.19) 
where ${a}_{0}$ is derived from Equation 7.15, $\varpi $ is in the same unit as ${a}_{0}$ and $P$ is the period in days to have ${f}_{\mathcal{M}}$ in solar masses, ${\mathcal{M}}_{\odot}$. Since ${a}_{0}$ corresponds to the orbit of the photocentre, ${f}_{\mathcal{M}}$ is a function of the masses of the components, ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$, and also their fluxes, ${\mathcal{F}}_{1}$ and ${\mathcal{F}}_{2})$, as shown in the equation:
$${f}_{\mathcal{M}}=\frac{{\parallel {\mathcal{F}}_{1}{\mathcal{M}}_{2}{\mathcal{F}}_{2}{\mathcal{M}}_{1}\parallel}^{3}}{{({\mathcal{F}}_{1}+{\mathcal{F}}_{2})}^{3}{({\mathcal{M}}_{1}+{\mathcal{M}}_{2})}^{2}}$$  (7.20) 
Since the more massive component is usually much brighter than the other, the index “1” is assigned to the brightest component, and ${\mathcal{F}}_{2}$ is neglected. The maximum value of ${f}_{\mathcal{M}}$ is then:
$${f}_{\mathcal{M}}({\mathcal{F}}_{2}\ll {\mathcal{F}}_{1})=\frac{{\mathcal{M}}_{1}{q}^{3}}{{(1+q)}^{2}}$$  (7.21) 
where $q={\mathcal{M}}_{2}/{\mathcal{M}}_{1}$ is the mass ratio. The largest possible value for ${f}_{\mathcal{M}}$, obtained when $q=1$ is then ${\mathcal{M}}_{1}/4$. This limit has been adopted, keeping in mind that it is a bit optimistic since components of the same mass should have identical fluxes, and therefore a zero semimajor axis.
Some of the orbital solutions derived in a trial processing are shown in Figure 7.6. It can be seen that the reliable solutions have a mass function of less than 0.3${\mathcal{M}}_{\odot}$, and that, as expected, periods close to oneyear are almost absent. For their part, the doubtful orbits populate very specific periods, such as 1 year, 6 months, or harmonics of the year and the precession period of the satellite (63 days), cf. Holl et al. (2023a).
Figure 7.6(b) shows the same solutions in the period vs. “significance of the parallax” (i.e. $\varpi /{\sigma}_{\varpi}$) diagram. It can be seen that the dubious solutions, in red, encroach on the domain of valid solutions, but that it is nevertheless possible to avoid them with the condition:
$$\varpi /{\sigma}_{\varpi}>\frac{\mathrm{20\hspace{0.17em}000}}{{P}_{\mathrm{days}}}$$  (7.22) 
Such a selection means in fact accepting shortperiod solutions only when they concern nearby stars. It has the advantage of effectively eliminating doubtful solutions (including those with small mass functions, but which are identifiable by their concentration on particular periods) without radically preventing the detection of binaries with a truly large mass function. For these reasons, it has been included among the acceptance criteria considered during processing, along with the basic criteria presented in Section 7.2.3.
Final selection of the orbital solutions
After processing, the alternative solutions that had been accepted with $F2>25$ were rejected, because the proportion of mass functions greater than 0.3${\mathcal{M}}_{\odot}$ was significantly higher than the value found for $$: 1.75 % instead of 1.05 %. After that, the rate of large mass functions could hardly be improved, but there remained many solutions concentrated around specific periods. Considering the diagrams linking the period to the other parameters, it appeared from Figure 7.7 that the following two selection conditions avoid most of these doubtful orbits. These are:
$$  (7.23) 
and a filter based on the significance:
$${a}_{0}/{\sigma}_{{a}_{0}}>158/\sqrt{{P}_{\mathrm{days}}}$$  (7.24) 
After this last filtering, 166 500 astrometric orbits were retained, including about 1% of mass function exceeding 0.3${\mathcal{M}}_{\odot}$.