Nii: a Bayesian orbit retrieval code applied to differential astrometry

The orbit retrieval procedure in Nii for a planet revolving around a target star is based on Bayesian analysis (Bayes, 1763). Let $\theta$ be the unknown parameter set that we want to estimate with a specific model $M$. Then, the posterior probability distribution for $\theta$ after we obtain observation data $D$ is of the following form:

where $p(D|\theta,M)$ is the likelihood reflecting the probability of generating the particular observation data $D$ if the parameters in the model $M$ are equal to $\theta$ and $p(\theta|M)$ is the prior distribution of $\theta$ representing our beliefs across different values of the parameters before observing the data. The denominator ${p(D|M)}$ is the normalization constant ensuring that the posterior distribution integrates to one. Given observation data $D$, we can simply omit the calculation of the denominator and evaluate the posterior distribution as

The main goal of orbit retrieval is to estimate the posterior probability distribution in Equation 7 by sampling. The point estimate of each parameter is set as the posterior mean because it is representative of the central position of the posterior distribution and mathematically accounts for the measure from a measure-theoretic perspective. Nii implements the APT-MCMC method (Liu, 2001; Gregory, 2005a), which is introduced in detail in Section 3.4. The underlying sampling algorithm of Nii’s MCMC strategy is Metropolis–Hastings (Metropolis et al., 1953; Hastings, 1970).

The choice of priors is very important in Bayesian inference since improper priors can produce misleading results. We assume that the prior distribution of each parameter is independent of that of each other parameter and set the prior distribution in a similar way to the radial velocity detection method (Gregory, 2005a).

In the case that the mass of the host star is known, every additional planet needs 7 more parameters to determine its orbit. In addition to the parameters describing planetary orbits, there is also an $\epsilon_{x}$ that represents any additional unknown errors. Table 2 shows the choice of priors for each parameter and their boundaries. We choose uniform priors for $e$, $\omega$, $\Omega$, and $M_{0}$, and the inclination $i$ is set to be uniform in $cosi$. For $M_{\mathrm{p}}$, the planetary orbital period $P$ and $\epsilon_{x}$, a uniform prior would be inappropriate since it would strongly favour a parameter range with larger values. We choose the Jeffreys prior for $M_{\mathrm{p}}$ and $P$ and a modified Jeffreys prior for $\epsilon_{x}$ (Gregory, 2005a).

For a simple model where the target star hosts only one planet, the joint prior distribution for the model parameters, assuming independence, can be written as

One of the features of the Nii code is that it considers the effect of using multiple-reference stars in the analysis of the observation errors. At all moments in an observation time series, we use the relative position angles between a target star and all the reference stars to calculate the standard deviation of the Gaussian observation error at each instant of time. Then, we separately pair the target star and each reference star to perform the orbit retrieval process, in which all the measurements of the wobbles between the target star and each reference star in the observation time series are used. Depending on the fitting results obtained using different reference stars, an additional voting process on the final fitted parameters may be required. Thus, the effect of using multiple reference stars in the differential astrometry method is to restrict the standard deviation at a single measurement time point and to cross-validate the final fitting results.

Suppose there are a total of $N$ measurements in an observation period $t$, and each measurement yields $\epsilon_{\eta}(t_{i})$ and $\epsilon_{\delta}(t_{i})$ between the target star and every reference star. If there are $N_{\mathrm{ref}}$ reference stars in total, then both $\epsilon_{\eta}(t_{i})$ and $\epsilon_{\delta}(t_{i})$ have $N_{\mathrm{ref}}$ elements. Assuming at each instant of time $t_{i}$ the measurement errors for the $N_{\mathrm{ref}}$ reference stars follow a Gaussian distribution with a mean of 0, we can calculate the standard deviation $\sigma_{\eta}(t_{i})$ and $\sigma_{\delta}(t_{i})$ in the Gaussian distribution based on the $N_{\mathrm{ref}}$ elements for both $\epsilon_{\eta}(t_{i})$ and $\epsilon_{\delta}(t_{i})$.

Since in differential astrometry a measurement gives the offset angles in two directions, namely, right ascension and declination, there are a total of $2N$ measurements between the target star and each reference star in the observation period $t$. Nevertheless, this is only an ideal situation, and some data may be missing in real cases. The likelihood function for the $2N$ measurements of the offset angles between the target star and each reference star is given by (Gregory, 2005a; Bishop, 2006)

where

Equation 9 demonstrates that the $N_{\mathrm{ref}}$ reference stars are used together only to determine the standard deviation of the distribution of the measurement errors at each instant of time $t_{i}$. In the orbit retrieval process, we separately fit the orbital parameters using the offset angles between the target star and each reference star in the entire time series $t$. This approach can reduce a sizable number of parameters in $\theta$ because each reference star has a different unknown error represented by $\epsilon_{x}$. Finally, we summarize the fitting results obtained from different reference stars to evaluate the reliability of the inferred orbital parameters.

The posterior distribution given by Equation 7 is highly multimodal, and an MCMC process can become stuck in many local optimal solutions and fail to fully explore the global distribution. One solution to this problem is referred to as simulated tempering (Geyer & Thompson, 1995; Gregory, 2005b), in which flatter versions of the target posterior distributions are generated using a temperature parameter $\beta$:

where $\beta$ varies from 1, corresponding to the coldest original posterior distribution, to 0, corresponding to the hottest distribution that is completely flat. Using a flat posterior distribution described by Equation 11, an MCMC sampler with a suitable $\beta$ can easily escape from local optimal solutions.

An attractive technique, parallel tempering (PT), runs several Markov chains with different $\beta$ values in parallel (Liu, 2001; Gregory, 2005b). In a PT-MCMC run, one chain with $\beta=1$ corresponds to the original posterior distribution, and other chains with different smaller $\beta$ values form a set of ladders leading to different higher temperature distributions. At random intervals, e.g., a state $s$, a pair of adjacent chains with temperature parameters $\beta_{i}$ and $\beta_{i+1}$ are chosen at random, and their current parameter states, $\theta_{s,i}$ and $\theta_{s,i+1}$, are interchanged with a probability of

Such swaps of parameter states allow parallel chains to exchange their sampling areas frequently, facilitating the convergence of global sampling.

In our orbit retrieval model, 8 MCMC chains are run in parallel with a set of $\beta=\{0.01,0.02,0.05,0.1,0.25,0.5,0.75,1.0\}$. A proposed swap is randomly evaluated using Equation 12 every 10 to 30 iterations.

One key factor that determines the completeness of MCMC sampling is the Metropolis–Hastings step size of each parameter $\theta_{i}$ (characterized by the standard deviation of the Gaussian proposal distribution $\sigma_{\theta_{i}}$). If the step sizes are too small, it is difficult for the MCMC sampler to cover the entire parameter space and find the optimal solution since the obtained posterior density is highly dependent on the starting location of the chain. In contrast, if the step sizes are too large, the MCMC sampler rejects the majority of proposals, and only a few effective samples are obtained. For a large number of parameters, an ideal combination of step sizes leads to an acceptance rate of $\sim$ 23.4% (Gelman et al., 1997).

In a high-dimensional orbit fitting model, it is difficult to find an ideal combination of sampling step sizes for all PT Markov chains. First, the ideal step sizes for different parameters $\theta_{i}$ are unique and change under different $\beta$ values. Second, when exploring different regions of the target posterior distribution, the ideal set of step sizes that is compatible with the shape of the local posterior distribution should also be adjusted accordingly. Therefore, in the case of the simplest single-planet model with 8 parameters, a PT-MCMC program with 8 chains with different $\beta$ values needs to tune 64 different step sizes in real time to reach the optimal acceptance rate, which is almost impossible to perform manually.

The Nii code adopts a different automatic control system from the existing model (Gregory, 2005a, b). In our APT-MCMC scheme, the step sizes of all the parameters of all the PT Markov chains are adjusted dynamically in real time after an initial burn-in stage. At the initial burn-in stage, we set the standard deviation of the Gaussian proposal of each parameter, $\sigma_{\theta_{i}}$ , to $10\%$ of the parameter’s prior range to ensure that the sampler can access all the high-density regions of the posterior distribution. However, these initial step sizes result in extremely low acceptance rates for Markov chains. Thus, after the burn-in stage, the Nii code monitors the acceptance rate of each Markov chain in real time and dynamically adjusts the standard deviations of the Gaussian proposals of all the parameters for the chains with a low acceptance rate. This process is implemented by dividing the entire Markov chains into many short cycles and adding many tuning stages for the chains with poor acceptance rates after each short cycle. In these tuning stages, the program separately adjusts the $\sigma_{\theta_{i}}$ of each parameter with a scaling array and observes the corresponding changes in the acceptance rate. Then, a new combination of step sizes is determined by combining two factors: the sensitivity of the acceptance rate to each $\sigma_{\theta_{i}}$ and the deviations between the ideal acceptance rate and the acceptance rates obtained by all the scaled parameters in the tuning stage. The strategy here is to select from the tuned parameter matrix the parameter set that can result in an acceptance rate that is closest to the ideal rate of $23.4\%$ in the local posterior distribution being sampled. Although this method cannot obtain the ideal acceptance rate of $\sim 23.4\%$, it is sufficient for efficient convergence. In our test runs, an APT-MCMC sampler can converge on the posterior distribution within 500,000 iterations.

We use the Gelman-Rubin criterion, $R_{c}$, to assess whether the APT-MCMC sampling process has converged (Gelman & Rubin, 1992). Gelman-Rubin diagnostics use an analysis of both the between-chain variance and within-chain variance of several independent Markov chains to assess convergence. Suppose there are $N$ Markov chains of length $L$. Let $\theta_{ij}$ denote the $i$th of the $L$ iterations of $\theta$ in the $j$th simulated chain, $\bar{\theta}_{j}$ denote the sample posterior mean of chain $j$, and $\bar{\theta}$ be the overall sample posterior mean of the $N$ chains. The between-chains variance, $B$, is given by

The within-chain variance, $W$, is given by

Then, we can obtain an unbiased estimator of the marginal posterior variance of $\theta$, $var(\theta)$, by taking the weighted average of $B$ and $W$:

Since the initial value of $\theta$ is set randomly during MCMC sampling, $var(\theta)$ should overestimate the true marginal posterior variance. On the other hand, $W$ tends to underestimate the within-chain variance in the early phase of MCMC sampling. As the sampling process gradually converges, both $var(\theta)$ and $W$ stabilize at the true variance of $\theta$. Thus, the Gelman-Rubin criterion uses the potential scale reduction factor, $R_{c}$, as a convergence diagnostic

In practice, the values of the Gelman-Rubin criterion should be less than 1.2 for model parameters to declare convergence in MCMC sampling, in more rigorous diagnostics, the values of $R_{c}$ are required to be less than 1.1 (Gelman & Rubin, 1992; Brooks & Gelman, 1998).

To assess the convergence of the model parameters in our Bayesian orbit retrieval framework, we run a batch of independent APT-MCMC sampling tests to calculate the Gelman-Rubin criterion. For simplicity, we carry out this process for only one reference star. Because all the reference stars are modelled as stable systems with no observable strong perturbations, the simulated differential astrometric signals derived from different reference stars all give the relative displacement of the target star. Only the random error is different. In this case, the $R_{c}$ value obtained from multiple independent APT-MCMC sampling runs using the simulated observation of one reference star can help us evaluate the convergence of the sampling process using the data of other reference stars.

The simplest case is to infer the orbital parameters of a system with only one planet. We simulate the differential astrometric signals of a planet with a mass of 3.1$M_{\oplus}$ and an orbital period of 541 days. The observation strategy of our mission is to make 300 measurements of the offset angles between the target star and its reference stars over a period of 5 years. We randomly remove several consecutive clusters from the entire time series to better simulate real observations, which generally contain discontinuities because the targets are not always visible from the satellite. As a consequence, we generate 8 pairs of $\sim$ 250-point time series of the angular wobbles in right ascension and declination using 8 relatively stable reference stars. The standard deviation of the simulated Gaussian measurement error is 1 $\mu\mbox{as}$.

Using the Nii code that implements the APT-MCMC sampling strategy described in Section 3, we fit the mass and orbital parameters of the simulated single-planet system using the differential angular measurements of the 8 reference stars. Our MCMC chains contain 1,000,000 iterations and throw the first 300,000 iterations as burn-in. To ensure the convergence of the sampling parameters, we select one reference star, perform independent APT-MCMC runs, and calculate the Gelman-Rubin convergence diagnostic criterion $R_{c}$. Table 3 gives the posterior means and the standard deviations of 8 independent Markov chains using the same observation data of one reference star and the $R_{c}$ values calculated from this set of independent APT-MCMC samplings. The results indicate that the convergence is very good for $P$, $M_{\mathrm{p}}$, $e$, $cosi$, $M_{0}$, and $\epsilon_{0}$, and that the $R_{c}$ values for these 6 parameters are all $\lesssim$ 1.01. In addition, $\omega$ and $\Omega$ do not converge separately, as they both have an $R_{c}$ that is greater than 1.2. However, the summation of $\omega$ and $\Omega$ converges with $R_{c}$ $\approx$ 1.01. The fitted means, standard deviations and $R_{c}$ values of the summation of $\omega$ and $\Omega$ in these 8 independent runs are also given in Table 3. By comparing the theoretical differential astrometry curves corresponding to a variable $\omega$, a variable $\Omega$, and different combinations of $\omega$ and $\Omega$, we find that similar observation curves can be obtained as long as the sum of the two variables remains the same. The argument of periapsis $\omega$ is related to the orbital eccentricity, and thus, its value has no meaning in the case of circular orbits. Likewise, the longitude of the ascending node $\Omega$ is related to the orbital inclination, and thus, its value is also meaningless for circular orbits. For the magnitudes of the orbital eccentricity and inclination in our simulated model, the differences between the differential astrometry curves corresponding to a wide variety of combinations of $\omega$ and $\Omega$ are minuscule when the summation of these two variables remains unchanged. Therefore, in our APT-MCMC sampling process, $\omega$ and $\Omega$ yield different combinations with similar sums, thereby causing these two parameters to fail to converge individually.

Figure 1 shows the corner plots (Foreman-Mackey, 2016) of all the model parameters, among which $\omega$ and $\Omega$ are displayed in the form of their summation, obtained using the differential astrometric signals of one of the reference stars in the single-planet model. We can see from the marginal posterior distributions that all the orbital parameters converge well. The corner plots also show that $M_{0}$ is strongly correlated with the sum of $\omega$ and $\Omega$. $M_{\mbox{p}}$ and $cosi$ are also correlated since both of these parameters determine the magnitude of the disturbance in differential astrometric signals.

Figure 2 compares the simulated differential astrometric signals corresponding to the reference star used in the fitting of Figure 1 and the theoretical signals generated by the marginal-posterior mean values of $M_{\mbox{p}}$, $a$, $e$, $cosi$, $P$, and $M_{0}$. The values of $\omega$ and $\Omega$ used to generate the theoretical signals are set to the values in the final iteration of the MCMC chain with $\beta=1$. These theoretical signals can accurately reproduce the angular wobbles in the simulated observations. Figure 2 further shows the residual plots of the fitted right ascension and declination curves. The mean values of the residuals in the right ascension and declination directions are -0.0671 and 0.0035 $\mu\mbox{as}$, respectively, and the corresponding standard deviations are 0.9795 and 0.9536 $\mu$as. These values are is in good agreement with the injected Gaussian observation error with a mean value of 0 $\mu$as and a standard deviation of 1 $\mu$as.

The retrieval processes using the simulated signals of the other 7 reference stars are also successful, and all the fitted results are summarized in Table 4. Among all the parameters, the orbital period $P$ achieves the best-fit: the largest error among the 8 cases is $<$ 1.5 days, and all cases have a one standard deviation of $\sim$ 1 day. The planetary mass $M_{\mbox{p}}$ is also accurately estimated with a maximum relative error of $\sim$ 4.8%. The other five orbital parameters are also reasonably derived. Because none of the 8 simulated reference stars is accompanied by an additional astrometric disturbance, the fitted unknown measurement errors $\epsilon_{x}$ in all of the 8 cases are less than 0.5 $\mu\mbox{as}$, which is less than the standard deviation of the injected Gaussian measurement error of 1 $\mu\mbox{as}$.

Retrieving the orbital parameters of a multiplanet system is far more complicated than retrieving those of a single-planet system, as the astrometric signals of the planets constructively or destructively interfere with one another due to their mutual gravitational interactions. Moreover, as shown in Equation 3, the strength of an astrometric signal is proportional to the product of a planet’s mass and its orbital semimajor axis. Therefore, when there are massive planets at distant orbits, the disturbances caused by low-mass planets at closer orbits become difficult to detect.

Here, we simulate a dual-planet system with two planets featuring different masses. The mass and orbital parameters of the low-mass planet, as well as the characteristics of the target star, are the same as those in the single-planet system case described in Section 4.1. Further, we add a massive planet of 25$M_{\oplus}$ at a larger distance into the dual-planet model; the orbital parameters of this high-mass planet are given in Table 5. In this system, the strength of the astrometric signal generated by the low-mass planet is approximately 10% of that of the massive planet. To include the mutual gravitational interaction between the two planets, the differential astrometric signal of the dual-planet system is simulated with the N-body package MERCURY6 (Chambers, 1999).

One correct approach to conduct the orbit retrieval of a dual-planet system is to simultaneously fit both planets’ orbital elements at the same time. Such an approach requires the combination of our APT-MCMC code with an N-body integration package. In the Nii code, instead of a more direct approach of using an N-body integration package to simulate the dual-planet system, we implement a simplified version of this method: we simultaneously fit both planets’ orbital elements using a linear superposition of the two planets’ Keplerian orbits. This process of simultaneously fitting the orbital elements of both planets involves 15 parameters, i.e., two pairs of the elements $(a,e,i,\omega,\Omega,M_{0},M_{\mbox{p}})$ and one parameter $\epsilon_{x}$ to account for additional measurement errors. However, we find that our APT-MCMC sampling strategy cannot converge for such a 15-parameter fitting process. Even if we set the initial values of the parameters near their true value, only the parameters of the massive planet can converge after 1,000,000 iterations. The main reason may be that the disturbance of the low-mass planet is too small compared to that of the massive planet. Therefore, it is difficult for our automatic control system to dynamically determine proper combinations of the sampling step sizes for the orbital parameters of the two planets corresponding to their different astrometric signal strengths. A more efficient sampling strategy is needed to achieve convergence of the two-planet model. We will investigate this in subsequent work.

Ultimately, we adopt a two-stage fitting process to retrieve the orbital elements of the dual-planet system. In the first round, we fit the two planets’ signals as though they originate from a single-planet system. As a result, the fitted unknown measurement errors $\epsilon_{x}$ are significantly larger than the standard deviation of the Gaussian measurement error. This implies that there may be other planets within the system. In the second round, we subtract the signals corresponding to the fitted massive planet with the parameters obtained in the first round of fitting and perform another round of fitting on the residuals. In both rounds, our MCMC chains contain 1,000,000 iterations and throw the first 300,000 iterations as burn-in.

Table 5 gives the results of these two rounds of fitting for all 8 reference stars. The orbital period $P$, $M_{0}$, and the sum of $\omega$ and $\Omega$ are well fitted in the first round. For the mass $M_{\mbox{p}}$, the fitted values are approximately 1.5% less than the injected value. The 8 fitting processes using different reference stars all yield a larger $e$ and a smaller $cosi$. An important indicator is the fitted unknown measurement errors $\epsilon_{x}$, which is $\sim$ 2.5 $\mu\mbox{as}$ for all 8 cases. Since the fitted values of $\epsilon_{x}$ are much larger than the standard deviation of the simulated Gaussian measurement error of 1 $\mu\mbox{as}$, other signals are inferred to exist within the system. In the second round of fitting, the main parameters of the low-mass planet, i.e., its planetary mass and orbital period, are well fitted, but the relative errors of the fitted values are larger than those in the case of the single-planet system in Section 4.1. Compared with the results of the single-planet system in Table 4, the other parameters are not well fitted, especially $e$ and $cosi$. We thus derive a much larger $e$ and a much smaller $cosi$. Because the fitted values of $\epsilon_{x}$ in the second round of fitting in the 8 cases are all less than 1.3 $\mu\mbox{as}$, similar to the standard deviation of the injected Gaussian measurement error, our fitting process ends after the second round.

The reason for the insufficient fits of $e$ and $cosi$ in the dual-planet system can be explained from Figures 3 and 4, which present the corner plots of all the model parameters. Once again, $\omega$ and $\Omega$ are displayed as their summation obtained in the first and second rounds of fitting using the differential astrometric signals of one of the reference stars in the dual-planet model. As shown in Figure 3, $M_{\mbox{p}}$ and $e$ are clearly correlated. In the first round, we seek to fit the astrometric signals of a dual-planet system using the disturbance of a single massive planet; consequently, the disturbance caused by the low-mass planet is initially explained by adjusting the eccentricity of the orbit. As a result, we obtain a larger orbital eccentricity and a smaller planetary mass in the first round of fitting. These inaccurate fitting values cause the residual after the first round of fitting to contain some distorted information, which further contributes to the deviations in the orbital eccentricity and orbital inclination in the second round of fitting. This outcome suggests that there should be a better way to fit the orbital elements of all the planets in a multiplanet system simultaneously. Unfortunately, our Nii code fails follow this approach, i.e., to fit all 15 orbital elements of the considered dual-planet system at the same time. Nevertheless, since our two-stage fitting approach can accurately obtain the masses and orbital periods of the two planets, which is the main information we seek, in this work, we continue to employ this imperfect method.

Figure 5 compares the simulated differential astrometric signals of the dual-planet system corresponding to the reference star used in the fitting of Figure 3 and 4 and the theoretical signals generated by the marginal posterior mean values of $M_{\mbox{p}}$, $a$, $e$, $cosi$, $M_{0}$, and $P$ in combination with $\omega$ and $\Omega$ in the final iteration in the MCMC chain. The mean values of the residuals in the right ascension and declination directions are -0.5072 and 0.1995 $\mu$as, respectively, and the corresponding standard deviations are 1.5358 and 1.1545 $\mu$as. We can see from the residuals that although the simulated observations can be accurately explained by the fitted curves overall, clear local features remain that cannot be reproduced by the fitted curves. These local features in the residuals originate from two aspects: one is the mutual gravitational action between the two planets that our fitting approach cannot explain, and the other is our inaccurate fits of the orbital eccentricity and inclination of the two planets.