Three direct imaging epochs could constrain the orbit of Earth 2.0 inside the habitable zone

TL;DR

This study assesses how multiple direct imaging observations can accurately determine the orbits of Earth-like exoplanets within habitable zones, informing future mission planning and target prioritization.

Contribution

It introduces a method using mock observations and MCMC fitting to evaluate the orbital parameter retrieval accuracy for exoplanets with realistic observational constraints.

Findings

Three imaging epochs can constrain semi-major axis within 5% for planets at 1 AU and 10 pc.

Optimal observation cadence is at least 90 days apart.

Orbital parameters can be estimated despite under-constrained data with semi-analytical methods.

Abstract

Space-based direct imaging missions (HabEx, LUVOIR) would observe reflected light from exoplanets in the habitable zones of Sun-like stars. The ultimate---but not sole---goal of these concept missions is to characterize such planets. Knowing an exoplanet's orbit would help two-fold: (i) its semi-major axis informs whether the planet might harbour surface liquid water, making it a priority target; and (ii) predicting the planet's future location would tell us where and when to look. The science yields of HabEx and LUVOIR depend on the number, cadence, and precision of observations required to establish a planet's orbit. We produce mock observations using realistic distributions for the six Keplerian orbital parameters, experimenting with both Beta and uniform eccentricity distributions, and accounting for imperfect astrometry ($\sigma$ = 3.5 mas) and obscuration due to the inner working angle of a high-contrast imaging system (IWA = 31 mas). Using Markov chain Monte Carlo methods, we fit the orbital parameters, and retrieve their average precisions and accuracies as functions of cadence, number of epochs, and distance to target. Given the time at which it was acquired, each image provides two data: the $x$ and $y$ position of the planet with respect to its star. Parameter retrieval based on one or two images is formally under-constrained; yet the semi-major axis posterior can be obtained semi-analytically. For a planet at 1 AU around a star at a distance of 10 pc, three epochs constrain the semi-major axis to within $\lesssim$5\%, if each image is taken at least 90 days apart.

Tables

Table 1. Table 1: The Keplerian orbital parameters and their priors. The prior range refers to the allowable parameter space explored in the Markov chain, while truth values are drawn from the input range.

Name	Symbol	Prior distribution	Prior range	Input
Semi-major axis	$a$	Uniform in natural log	[0.01, 50] AU	[0.95, 1.70] AU
Eccentricity	$e$	Beta ($\sigma$=0.081) or uniform	[0, 1)	[0, 1)
Inclination	$i$	Uniform in cosine	[0, $\pi$]	[0, $\pi /2$]
Argument of periapsis	${\omega }_p$	Uniform	[0, $2\pi$)	[0, $2\pi$)
Longitude of the ascending node	$\Omega$	Uniform	[0, $2\pi$)	[0, $2\pi$)
Mean anomaly at first epoch	$M_0$	Uniform	[0, $2\pi$)	[0, $2\pi$)

Table 2. Table 2: Z-score statistics for data in figures 4 and 5 . Cases (a) and (c) are fine, but (b) is dangerous.

	(a) Beta		(b) Mixed		(c) Uniform
	mean	stdev.	mean	stdev.	mean	stdev.
Epoch 3	-0.47	1.19	-1.02	2.56	-0.18	1.26
Epoch 4	-0.10	0.99	-0.28	2.59	-0.14	1.11
Epoch 5	0.15	1.12	0.23	4.11	-0.13	1.10

Equations

$$\ln L=-\frac{{\chi }^2}{2}-k(\ln \pi +2\ln {\sigma }_{xy}),$$

$${\chi }^2=\sum _{i=1}^k\frac{(x_{\mathrm{obv},i}-x_{\mathrm{model},i}{)}^2+(y_{\mathrm{obv},i}-y_{\mathrm{model},i}{)}^2}{{\sigma }_{xy}^2}.$$

$$\mathrm{Prob}(a\mid a_{\mathrm{proj}},e=0)=\frac{1}{{\beta }^2\sqrt{{\beta }^2-1}},\beta \in [1,\infty ),$$

$$\left[\begin{array}{c}{x}_k\\ y_k\\ z_k\end{array}\right]=\left[\begin{array}{c}a\cos (M_k)-c\\ b\sin (M_k)\\ 0\end{array}\right],$$

$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right]=\left[\begin{array}{ccc}\cos \gamma & -\sin \gamma & 0\\ \sin \gamma & \cos \gamma & 0\\ 0& 0& 1\end{array}\right]\cdot \left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \beta & -\sin \beta \\ 0& \sin \beta & \cos \beta \end{array}\right]\cdot \left[\begin{array}{ccc}\cos \alpha & -\sin \alpha & 0\\ \sin \alpha & \cos \alpha & 0\\ 0& 0& 1\end{array}\right]\cdot \left[\begin{array}{c}x\\ y\\ z\end{array}\right],$$

$$M_k=M_0+\sqrt{\frac{\mu }{a^3}}k\delta t\frac{86400\sec }{\mathrm{day}},$$

Full text

planets and satellites: detection — planets and satellites: terrestrial planets — telescopes

††journal: AJ

\turnoffedit

1 Introduction

The first generation of space-based direct imaging missions would observe planets in the habitable zones of Sun-like stars. As recommended by the National Academies of Sciences, Engineering, and Medicine’s recent Exoplanet Science Strategy, “NASA should lead a large strategic direct imaging mission capable of measuring the reflected-light spectra of temperate terrestrial planets orbiting Sun-like stars.”111https://www.nap.edu/read/25187/ To acquire these spectra, we must first plan the missions cost-effectively.

A single direct imaging observation would distinguish Earths from similarly-bright sub-Neptunes with a one-in-five success rate, due to degeneracies between radius, albedo, and phase. If the planet’s semi-major axis is known, on the other hand, the success rate is one in two.

In this paper we help prepare for future direct imaging missions such as HabEx and LUVOIR by determining how many visits are required to each newfound planet before we can confidently know its semi-major axis. In doing so, we hope to quantify a key input parameter in more comprehensive mission simulations (see Stark et al., 2018).

1.1 Images to orbits

A single epoch of direct imaging provides the two-dimensional position of the planet in the sky plane relative to its host star: $x$ and $y$, where the star is at the origin. Six parameters are needed to uniquely describe the three-dimensional Keplerian orbit of a planet (listed in Table 1). Given the time at which it was taken, each image constrains $x$ and $y$, at most three measurements should be required to fit all six parameters.

The fitting methods employed in previous direct imaging orbit retrieval efforts come from either Markov chain Monte Carlo (MCMC) or Bayesian rejection sampling (Pueyo et al., 2015; De Rosa et al., 2015; Rameau et al., 2016; Wang et al., 2016, 2018; Blunt et al., 2017; Kosmo O’Neil et al., 2018). The latter type are carefully demonstrated in Blunt et al. (2017) for closely-spaced observations over tiny orbital coverage ($\sim$3%).

Mede & Brandt (2017) present a software package to simultaneously fit direct imaging observations and radial velocity observations. Such an approach may be fruitful for eventual Earth twins if high precision radial velocity can be sensitive to the 10 cm/s signals (Fischer et al., 2016).

As for future direct imaging, previous design reference missions for HabEx and LUVOIR have not always implemented an optimized number or cadence of observations. Stark et al. (2014) look only at single-visit yields; Stark et al. (2016, 2015) do not consider that a mission would revisit candidates for the strict purpose of establishing their orbits, assuming that most stars would be revisited regardless to increase the total yield.

In this paper, we anticipate space-based observations from HabEx and LUVOIR. The work presented here is fundamentally distinct from these earlier efforts in that we explicitly quantify the number, cadence, and precision of observations required to establish, for targets at any distance, the orbit of a planet within its star’s habitable zone.

1.2 A note on habitable zones

The circumstellar habitable zone (HZ) is a theoretical shell around a star within which planets can harbour liquid water at their surfaces. The inner and outer edges of the HZ have been variously modeled based on assumptions about climatic systems. Early habitable zone estimates (Hart, 1979) were physically driven by the destabilizing ice-albedo and water vapour feedbacks and were narrower than in current conceptions.

The HZ estimates we use come from physics gleaned from the Earth. Assuming the surface of a planet has some exposed silicate rock, CaSiO3 minerals in the rock will react to consume atmospheric CO2. The rate of this reaction increases exponentially with temperature, so the planetary effect of warming the atmosphere is to sequester more carbon and weaken the greenhouse effect (Walker et al., 1981; Kasting, 1988; Kasting & Toon, 1989). This is the well-known silicate weathering feedback, and is the basis of the HZs proposed in Kasting et al. (1993) and updated in Kopparapu et al. (2013).

This conception of the HZ has yet to be empirically validated (but it could through statistical study of exoplanet properties; Bean et al., 2017). If a planet does not have exposed surface rock, for example, then the HZs given by a model based on the silicate weathering feedback might not apply (e.g. Abbot et al., 2012). We therefore base our analysis on a variety of “hypothetical” HZ widths. The semi-major axis precision that would satisfy us depends on the sizes of these HZs relative to our semi-major axis estimate.

2 The orbit-retrieval model

This work quantifies, statistically, the orbit retrieval accuracy and precision for many simulated HZ planets. We simulate planets with each orbital element randomized. For each epoch of each planet, we calculate the planet’s $(x,y)$ position relative to its star, given a fixed cadence in days and with added Gaussian noise—here we adopt $\sigma_{\theta}=3.5$ mas as the baseline astrometric precision (HabEx Team, 2018; LUVOIR Team, 2018). Then we retrieve the posteriors on orbital elements from our synthetic dataset using MCMC. We repeat this numerical experiment under varying assumptions, described below.

Eccentricity distribution

Since the eccentricity distribution of HZ terrestrial exoplanets is unconstrained, we repeat our experiment assuming different eccentricity distributions—a beta distribution with parameters $a=0.867$, $b=3.03$ (Kipping, 2013; Nielsen et al., 2008), or a uniform distribution—for both the underlying true distribution, and the prior distribution used in the retrieval.

Cadence

The retrieval experiments are repeated, using the same synthetic orbits, for cadences of 30, 90, 180, and 270 days. This reveals which cadence gives the best fit, for the same number of epochs. Our result is optimized for orbital periods within the G2V habitable zone: we expect the best cadence to be set by how much the planet has moved, which is some fraction of the period. In principle, a mission simulation would consider data from previous epochs to choose when to revisit a star. We take the simpler approach of a fixed cadence to obtain a broad result, which could conveniently be applied to mission prognoses.

Distance

The distance of the target star will impact orbit retrieval through astrometric error and the inner working angle (IWA) of starlight suppression, both of which correspond to greater projected separations at greater distances. We take a baseline distance of 10 pc, and repeat our experiment for distances of 5 pc and 20 pc.

Each scenario considers 100 planets with orbital parameters $a$, $i$, $e$, $\omega_{p}$, $\Omega$, and $M_{0}$, sampled from the underlying distributions (Table 1). These are the prior distributions described in section 2.4—excepting $e$, which we allow to have different underlying and prior distributions, as detailed later.

To retrieve posterior probability distributions of the orbital parameters, we use emcee, an MCMC ensemble sampler (Foreman-Mackey et al., 2013). The goal of any MCMC implementation is to evaluate the posterior Prob(model|observation, $\sigma$) $\propto$ Prob(observation|model, $\sigma$) $\times$ Prob(model), where $\sigma$ is the measurement precision. Given estimates of the right hand side, we can estimate the left hand side. Thus setting up emcee requires:

Estimates of the prior probability distribution, Prob(parameter), for each orbital parameter, 2. 2.

A forward model that calculates the planet positions $(x_{0},y_{0}),(x_{1},y_{1}),...,(x_{k},y_{k})$ up to the $k$th epoch given orbital parameters, 3. 3.

A likelihood function, Pr(observation|parameters, $\sigma$), that calculates the probability of the observed $(x,y)$ position given the forward model.

The MCMC uses 30 walkers, randomly initialized in a Gaussian ball around the best-fit parmeters from a quick scipy.optimize likelihood maximization. The walkers have a burn-in time of 1000 steps and run for up to $5\times 10^{5}$ steps. The rest of the MCMC setup is detailed in subsections 2.1 through 2.4.

2.1 Likelihood function for detections

The likelihood of the observed position of the planet at some epoch is a normal distribution centred on the true position with width $\sigma_{xy}$. If we have $k$ images in which the planet is detected, then the log likelihood is

[TABLE]

where

[TABLE]

2.2 Likelihood function for non-detections

The first epoch is—by definition—a detection, as we are concerned with the orbit retrieval of actual planets rather than the odds of spotting a planet in the first place. Subsequent epochs may yield non-detections, however.

A planet may go undetected in an image for two reasons: (i) it is imaged at small projected separation and is occluded along with the starlight, or (ii) it does not reflect enough light. We ignore the latter case: the brightness of the planet depends on the geometric albedo and the scattering phase function, about which we presume nothing. Titan, for example, appears brighter at crescent phase, due to a distinctly non-Lambertian phase function (García Muñoz et al., 2017).

The relevant criterion in (i) is the IWA of the telescope, defined as the point where photon transmission through the instrument has decreased by 50%. We adopt an IWA of 31 mas. For a planet of given size, albedo, and phase function, and given a contrast floor, there will be a unique IWA at which it becomes undetectable at gibbous phase, and another where it becomes undetectable at crescent phase. We ignore brightness information and hence adopt a single hard-cutoff IWA.

A non-detection therefore still provides a constraint; it means that the angular separation of the planet and the star is less than the IWA. In this way, a nondetection is analogous to a measurement centered at the origin with an uncertainty of $\theta_{\rm IWA}$. If the planet is not detected at some epoch, then any set of parameters that would place it outside the IWA are assigned zero likelihood.

2.3 Posterior probability for less than three detections

Given only one or two detections, we expect MCMC to converge slowly, if at all, since the problem is formally under-constrained. Thus we eschew Kepler’s laws for a naive, semi-analytic approach, as long as the number of detections is less than three.

That is: for a given detection, the posterior probability distribution of $\ln(a)$ is a strong function of the projected separation, $a_{\rm proj}$. For circular orbits we can express this analytically:

[TABLE]

where $\beta=a/a_{\rm proj}$. Nonzero eccentricities complicate this analytic distribution, but we can model it numerically for given values of $a_{\rm proj}$ and $a_{\rm IWA}$ by drawing 1000 random orbits from the priors in Table 1 with $a_{\rm proj}$ fixed at the observed value. As figure 2 shows, the distribution of semi-major axes of these orbits peaks at the true $a$, and is hence a reasonable proxy for the posterior, for a given measurement of $a_{\rm proj}$. If a planet is not detected at some epoch, we also simulate random orbits, but with $a_{\rm proj}\leavevmode\nobreak\ \in\leavevmode\nobreak\ (0,a_{\rm IWA}$].

For two or more epochs, we construct “independent” posteriors for each $a_{\rm proj}$ observation, and multiply them together. The resulting joint probability distribution usually has one peak. Of course these measurements would not be truly independent, since they are linked by a Keplerian orbit. Thus our semi-analytic posterior distributions give an upper limit on the true uncertainty.

2.4 Prior distributions

As summarized in Table 1, parameterizations are chosen such that all priors save for eccentricity are flat. The log-uniform prior on semi-major axis is roughly consistent with recent literature (Petigura et al., 2013; Foreman-Mackey et al., 2014; Silburt et al., 2015; Burke et al., 2015; Christiansen et al., 2015; Kopparapu et al., 2018). We adopt a wide prior range of $a\in[0.1,50]$ AU.

The mean anomaly at the time of the first epoch is uniformly distributed in $[0,2\pi)$, as are arguments of periapsis and longitudes of the ascending node. Inclination is uniform in $\cos{i}\in[-1,1]$.

Eccentricity is drawn from a beta distribution with parameters $a=0.867,b=3.03$, as determined for radial velocity planets by Kipping (2013) and in agreement with Nielsen et al. (2008). Because this parameter’s underlying distribution is especially hard to constrain, we repeat our experiment with a uniform prior and uniformly-distributed true eccentricities, and again with a beta prior and uniformly-distributed true eccentricities (where we would be over-confident).

3 Results

Figure 1 illustrates the results of our orbit retrieval for five epochs of a planet at 20 pc with a 90-day cadence. This particular planet goes undetected in the second epoch, yet we still see an increased similarity between the retrieved and true orbits. In this example, the best-fit orbit after each epoch is the Markov chain link with the maximum posterior probability.

3.1 Semi-major axis retrieval under different eccentricity scenarios

The three panels in figures 4 and 5 show the accuracy and precision—respectively—of semi-major axis retrieval, for different prior assumptions and underlying distributions of orbital eccentricity. The distributions are either both beta (realistic), both uniform (pessimistic), or beta-distributed in the prior and uniform in the underlying truths (i.e., a pathological scenario which may underestimate error).

Figure 4 shows the accuracy of the retrieved $a$ estimates as a function of epochs for a sample of 100 planets. Errors are retrieved from semi-analytic posteriors of $a$ for the underconstrained first two epochs, while MCMC is used for the third epoch and up. “Accuracy” is the difference between the median of the posterior and the true semi-major axis.

Figure 5 is complementary to figure 4 and presents the precisions on the retrieved semi-major axes, along with the sample median of these precisions. Again, MCMC is used for the third epoch on. “Precision” is the half-width of the 1$\sigma$ confidence interval, and it shrinks with each additional observation. This is as expected by degrees of freedom: the inflection point near the third epoch indicates that here we begin to gain less precision with additional measurements.

Our baseline case is shown in the left panel of figures 4 and 5. Here, by the second epoch we achieve quite good precision. At $1\sigma$ confidence, we can constrain $a$ to <25% with two epochs, and to $\lesssim$10% with three epochs.

The central subplots in figures 4 and 5 demonstrate a case where our precision is overstated; the prior eccentricity distribution is narrower than the underlying distribution. As shown in Table 2, the standard deviation of the z-scores for this scenario of amiss “mixed” distributions is much greater than unity, which confirms the overstatement.222The z-score of a fit, defined here as the ratio between our accuracy and precision, has a standard deviation of about unity if the fit precision reliably encompasses the true value, and the retrieved values are Gaussian. Therefore, if we were to assume a narrower eccentricity distribution than nature provides, we would not be able to believe our retrieved semi-major axis precisions.

The other endpoint draws eccentricity guesses from a uniform distribution (right panel of figures 4 and 5). Despite this pessimistic prior, we can still constrain the semi-major axis to 25% by the third epoch, and by the fifth epoch, retrieval accuracies are not much different from the best-case scenario.

3.2 Optimal cadence

We ran our experiment for different observational cadences (30, 90, 180 and 270 days) to see which retrieves the most accurate orbits in the least number of epochs. Here, host stars are all at 10 pc and planets have beta-distributed eccentricities. Our “optimal cadence” pertains to planets with periods in or near the G2V habitable zone.

Results are shown in figures 6 and 7—by the third epoch, $a$ can be constrained to within 10% of its true value for 68% of samples only for cadences 90 days or greater. In half a period, the planet will move to the furthest possible position from its original location. Revisiting the planet at its antipodal point means that we are less likely to miss the planet inside the IWA (e.g., compared to revisiting a quarter-period later), although this is tied to us ensuring the first epoch is a detection.

180- and 270-day cadences result in the best accuracy on semi-major axis at epoch 3, but 180-day cadences for near-360-day orbits may not break formal degeneracies between $\omega_{p}$ and $\Omega$. Hence the optimal cadence is not necessarily the same for $a$ as for other orbital parameters. While the 270-day cadence also performs well, in practice this long wait could risk losing track of the planet, as discussed below.

3.3 Effects of distance and astrometric precision

We have adopted a distance of 10 pc and a measurement error of $\sigma_{\theta}=3.5$ mas for our numerical experiments thus far. However, the number of epochs needed to constrain an orbit is distance-dependent. More-distant targets impede orbit retrievals in two ways: the corresponding measurement error on planet position increases (by a factor of four from 5–20 pc); and the larger projected separation of the IWA is more likely to obscure targets (nearly 0% nondetections at 5 pc to $\sim$70% at 20 pc).

To quantify the effect of distance, we repeated the retrieval experiments for the same planets at 5 pc and 20 pc, assuming 90-day cadences. At our assumed IWA and range of true semi-major axes, HZ planets will be obscured 100% of the time by 55 pc. This experiment is equivalent to increasing the IWA and astrometric uncertainty. Particular sources of astrometric error and their impact on orbit retrievals are discussed in Pueyo et al. (2015).

Figure 8 illustrates the change in semi-major axis retrieval accuracy with distance. As expected, the width of the 1$\sigma$ confidence interval increases for more distant targets, for all epochs.

The encouraging performance of the fit at 5 pc is attributed to the small astrometric error. For these nearest targets we would be capable of measuring $a$ to within <2% by the fourth epoch.

We are only showing this distance dependence of accuracy for one cadence. The optimal cadence is distance-dependent; with increasing distance, obscuration dominates measurement error in worsening the accuracy, preferring 180-day (antipodal) cadences.

4 Discussion

The minimum number of epochs we need depends on the science question we are asking.

4.1 Placing a the planet in the habitable zone

We first consider hypothetical HZs defined only by how narrow they are. Compare the widths of these hypothetical HZs with the semi-major axis retrievals of planets on 1-AU orbits (figure 9, left panel): at 1$\sigma$ confidence, for an optimistic, 25%-wide HZ centered on 1 AU, a single epoch of direct imaging generally suffices to nail the semi-major axis accurately. For a pessimistic 1%-wide HZ, four or five epochs are needed.

To compare these results with theoretical HZs from the literature: the right panel in figure 9 shows the same semi-major axis constraints, to scale with both a LUVOIR-esque inner working angle and with the theoretical HZ inner limits for solar twins from Kopparapu et al. (2013). By the fourth epoch we are 68% sure a planet is beyond the Moist Greenhouse inner limit. To constrain an orbit beyond the more optimistic Recent Venus inner limit at 0.75 AU at the same confidence, about one epoch is required.

We are implicitly assuming that the HZ is a function of semi-major axis only. While high eccentricities may negatively affect the stability of surface water (Bolmont et al., 2016), the high thermal inertia of extant oceans can buffer against transitions to “snowball” states at apoastron (Dressing et al., 2010). Either way, the long-term stability of a planetary climate depends, to the first order, on the average incident flux over the entire orbit (Williams & Pollard, 2002), which is proportional to $(1-e^{2})^{-1/2}$.

4.2 Predicting a planet’s future location

In some cases, we would need to constrain all six orbital elements. For example, knowing the entire orbit would help us predict the exoplanetary ephemerides. The best theoretical precision on extrapolating a planet’s position is set by the detector pixel scale, although in practice, the precision we achieve depends on the post-processing. Using three epochs, we have predicted the position of a planet at the fourth epoch to within a given number of pixels on average. For 180-day cadences, we find the position to be constrained to within one pixel; for 90-day cadences, within two pixels, with five epochs necessary to get down to a one-pixel accuracy. The trade-off is that waiting longer between observations provides better accuracies for the Keplerian parameters, but any errors are amplified by the greater distance the planet will have travelled.

On a related note, we have only considered single-planet systems in our study. However, we expect a substantial fraction of planets to occur in multi-planet systems (Zhu et al., 2018). This presents an interesting tension with respect to the optimal cadence: for multi-planet systems, we may prefer shorter observational cadences because the planet will move less, so there is less risk of confusing it with other planets. Nevertheless, planet disambiguation in multi-planet systems will be made a bit easier by the fact that we can expect these targets to have near-circular orbits (Van Eylen et al., 2018).

The question this work has asked is: how well can one locate a planet in the habitable zone? Under the assumption of one planet per star, we find that 180- and 270-day cadences have the best precision and accuracy. In a realistic mission, however, these long cadences would likely be poor at resolving confusion between multiple planets in the same system. A 90-day cadence is a good compromise.

Even if the third image placed a planet in the habitable zone at 95% confidence, we would still need more epochs to establish its orbit before pursuing expensive spectroscopy. Yet the point is that we could begin prioritizing targets after only one or two epochs. This quantifies a key parameter in Design Reference Missions for future direct imaging concept missions HabEx and LUVOIR. In the LUVOIR study, constraining orbits within the habitable zone is the third step to identifying a habitable planet, after (1) establishing the target star list and (2) performing multi-colour point source photometry to rule out background objects (LUVOIR Team, 2018). Our work finds that the orbit-constraining step can be done more efficiently than before, reducing the minimum number of observations from six to three. This extra time could be spent characterizing more planets.

This work is supported by the McGill Space Institute, the TEPS training program, and an NSERC Discovery Grant. The authors wish to thank the HabEx AEIWG for their useful conversation: especially Leslie Rogers, Eric Nielsen, and Scott Gaudi. Comments from an anonymous reviewer also greatly improved the quality of this manuscript. All research was conducted on the territory of the Kanien’kehá:ka, the keepers of the Eastern Door of the Haudenosaunee Confederacy.

Appendix A Forward model

We assume that all planets are on bound Keplerian orbits.

If we hover with the planet’s orbital plane below us and our right ears towards the +$x$ reference direction, we will “see” an ellipse. Its shape can be described parametrically as the first two dimensions of a three-dimensional matrix:

[TABLE]

where $a$ is the semi-major axis, $b=a\sqrt{1-e^{2}}$ is the semi-minor axis, $c=\sqrt{a^{2}-b^{2}}$ is an $x$-intercept, $M_{k}$ is the mean anomaly at the $k$th epoch, and the orbital plane is defined by $\bm{x}$-$\bm{y}$.

The ellipse we see in the detector plane, $\bm{x^{\prime}}$-$\bm{y^{\prime}}$, has been rotated through three angles in the order $\Omega$, $i$, $\omega_{p}$. There is a Euclidean rotation matrix corresponding to this series of rotations:

[TABLE]

where $\alpha=\Omega$, $\beta=i$, and $\gamma=\omega_{p}$.

The projected separation in the image plane is the Euclidean distance, $\sqrt{(x^{\prime})^{2}+(y^{\prime})^{2}}$.

A.1 Planets move

Both the cadence, $\delta t$, and the walker’s current “guesses” of $a$ and $M_{0}$ control the planet’s position at the $k$th epoch (i.e., set by the instantaneous mean anomaly $M_{k}$):

[TABLE]

where $\sqrt{{\mu}/{a^{3}}}$ is the mean motion of the planet (in radians per second) and $\mu$ is the standard gravitational parameter.