DIRECT IMAGING IN THE HABITABLE ZONE AND THE PROBLEM OF ORBITAL MOTION

The HZ is generally agreed to be the region around a star where a planet can have liquid water on its surface. This is far from simply related to the blackbody equilibrium temperature, as it depends on atmospheric composition and the action of the greenhouse effect (Kasting et al. 1993; Kopparapu et al. 2013), among other factors. For our purposes it is enough to assume that the HZ is generally located at about one AU from a star, scaled by the star's luminosity:

$$\begin{equation} a_{{\rm HZ}} \approx \sqrt{L*/L_{\odot }}{\rm AU}. \end{equation} \tag{ 1 }$$

Traub (2012) provided three widths for the HZ based on various considerations, and then used the first 136 days of data from the Kepler mission to estimate that the fraction of Sun-like stars (spectral types FGK) with an Earth-like planet in the HZ is η_⊕ ≈ 0.34. More generally, this analysis indicates that η_planet ≈ 1.2, implying that every Sun-like star is likely to have a planet in its HZ, and some will have more than one. While this exciting result is based on a very large extrapolation from the earliest Kepler results, it is currently one of our best estimates of planet frequency in the HZ.

This topic was recently brought to the fore with the announcement of α Cen Bb by Dumusque et al. (2012). Discovered using the RV technique, α Cen Bb is an msin i = 1.13 M_⊕ planet orbiting a K1 star at 0.04 AU. While certainly not in the HZ, this discovery has exciting implications for the presence of planets in the HZ of the nearest two Sun-like stars.

The above arguments hint that planets will be common in the HZ of Sun-like stars. We are about to enter a new era of exoplanet direct imaging. With the next generation of giant telescopes and high-performance spaced-based coronagraphs we will be searching for planets in this scientifically important region around nearby stars.

The typical search for exoplanets with direct imaging has used 2.4 m (Hubble Space Telescope, HST) to 10 m (Keck) telescopes. These surveys have mostly concentrated on young giant planets, which are expected to be self-luminous as they dissipate heat from their formation. This allows them to be detected at wider separations from their host stars, where reflected starlight would be too faint. This has also caused planet searches to typically work at H band (∼1.6 μm), with exposure times of ∼1 hr. Examples conforming to these archetypes include Lowrance et al. (2005) using HST/NICMOS, the Gemini Deep Planet Search (Lafrenière et al. 2007), the Simultaneous Differential Imaging survey using the Very Large Telescope and MMT (Biller et al. 2007), the Lyot Project at the Advanced Electro-Optical System telescope (Leconte et al. 2010), the International Deep Planet Survey (Vigan et al. 2012), and the Near Infrared Coronagraphic Imager at Gemini South (Liu et al. 2010).

These searches have had some success. Examples include the four planets orbiting the A5V star HR 8799 (Marois et al. 2008b, 2010), with projected separations of 68, 38, 24, and ∼15 AU. These correspond to orbital periods of ∼460, ∼190, ∼100, and ∼50 yr, respectively. The A5V star β Pic also has a planet (Lagrange et al. 2010) orbiting at ∼8.5 AU with a period of ∼20 yr (Chauvin et al. 2012). Another A star, Fomalhaut, has a candidate planet on an 872 yr (115 AU) orbit (Kalas et al. 2008). At these wide separations it takes months, or even years, to notice orbital motion.

In the much closer HZ, however, orbital periods will be on the order of one year. We show in some detail that this is fast enough to yield projected motions of significant fractions of the point-spread function (PSF) FWHM over the course of an integration. The resulting smeared out image of the planet will have a lower S/N, making our observations less sensitive.

In addition to HZ planets having higher orbital speeds than the current generation of imaged exoplanets, integration times required to detect them will be much longer. Direct imaging surveys to date have mostly worked in the infrared while attempting to detect young planets still cooling after formation. The coming campaigns to image planets in the HZ of nearby stars will focus on older planets, which will be less luminous in the near-infrared. In the HZ, starlight reflected from the planet will be more important. The result is integration times required to detect such planets will be tens of hours, rather than the ∼1 hr characteristic of current campaigns.

Consider the Exoplanet Imaging Camera and Spectrograph (EPICS), an instrument proposed for the E-ELT. Kasper et al. (2010) predicted that EPICS will be able to image the RV detected planet Gl 581d, which has a semi-major axis of 0.22 AU with a period of ∼67 days (Forveille et al. 2011; Vogt et al. 2012). This orbit places it on the outer edge of the HZ of its M2.5V star (von Braun et al. 2011). EPICS will be able to detect Gl 581d, at a planet/star contrast of 2.5 × 10⁻⁸, in 20 hr with S/N = 5 (Kasper et al. 2010). Since this is a ground-based instrument, a 20 hr integration will be broken up over at least two nights. Plausible observing scenarios could extend this to several nights, taking into account such things as the need for sky rotation. As we will show, the planet will move several FWHM on the EPICS detector during a multi-day observation.

More generally, Cavarroc et al. (2006) showed that when realistic non-common path wave front errors are taken into account, the integration times required to achieve the 10⁻⁹ to 10⁻¹⁰ contrast necessary to detect an Earth-like planet around a Sun-like star approach 100 hr on the ground, even on a 100 m telescope with extreme adaptive optics (AO) and a perfect coronagraph. One of several concerns about the feasibility of a 100 hr observation from the ground is that such a long observation will be broken up over many nights.

With net exposure times of 20–100 hr, and total elapsed times for ground-based observations of several to tens of days, HZ planets will move significantly over the course of a detection attempt. The focus of this investigation is the impact of the orbital motion of a potentially detectable planet on sensitivity.

Though it has not yet been a significant issue in direct imaging of exoplanets, orbital motion has been considered in several closely related contexts. Here we briefly review a select portion of the literature. A very similar problem has been addressed in the context of searching for objects in our solar system, such as Kuiper Belt Objects (KBOs), which can have proper motions on the order of 1'' to 6'' per hour (Chiang & Brown 1999). Blinking images to look for moving objects by eye is a well-established technique. A more computationally intensive form of blinking images proceeds by shifting-and-adding a series of short exposures along trial paths, usually assumed to be linear. This "digital tracking" makes it possible to detect KBOs too faint to appear in a single exposure. This has been done both from the ground (Chiang & Brown 1999; Yamamoto et al. 2008) and from space with HST (Bernstein et al. 2004). More recently Parker & Kavelaars (2010) have taken into account nonlinear motion and optimized selection of the search space, especially important given the large data sets that facilities such as the Large Synoptic Survey Telescope will produce.

Orbital motion is an important consideration when planning coronagraphic surveys of the HZs of nearby stars. Brown (2005) treats the problem of completeness extensively. Large parts of the HZ will be within the inner working angle of the Terrestrial Planet Finder-Coronagraph (TPF-C) and so undetectable during a single observation. Also discussed in Brown (2005) is photometric completeness—that is how long the TPF-C must integrate on a given star to detect an Earth-like planet in the HZ. Other work on this topic includes Brown & Soummer (2010) and Brown (2004). These analyses consider orbital motion only between observations, not during a single observation as we do here. In general, the scenarios considered for these studies involved space-based high-performance coronagraphs on medium to large telescopes. In such cases exposure times were short enough and continuous so that orbital motion should be negligible during a single observation.

The work most similar to our analysis here is the detection of Sirius B at 10 μm by Skemer & Close (2011); in fact, it was part of our motivation for the present study. Skemer & Close (2011) used the well-known orbit of the white dwarf companion to Sirius to de-orbit 4 yr worth of images. Before accounting for orbital motion, Sirius B appeared as only a low S/N streak, but after shifting based on its orbit it appears as a higher S/N point source from which photometry can be extracted. Similar to this method, we will analyze the prospects for de-orbiting sequences of images, only we consider the case with no prior information at all, and with orbital elements with significant uncertainties.

We begin by considering a focal plane detector working at a wavelength λ in μm. The FWHM of the PSF for a telescope of diameter D in m, neglecting the central obscuration, is

$$\begin{equation} {\rm FWHM} = 0.2063\frac{\lambda }{D}{\rm arcsec}. \end{equation} \tag{ 2 }$$

If we are observing a planet in a face-on circular (FOC) orbit with a semi-major axis of a in AU at distance d in pc, its angular separation will be a/d arcsec. At the focal plane the projected separation will then be

$$\begin{equation} \rho = 4.847\frac{a D}{\lambda d} {\rm in FWHM}. \end{equation} \tag{ 3 }$$

We note that it will occasionally be convenient to specify ρ in AU instead of FWHM. When it is not clear from the context we will use the notation ρ_au to denote this.

The orbital period is $P=365.25\sqrt{a^3/M_*}$ days around a star of mass M_* in M_☉. In one period, the planet will move a distance equal to the circumference of its orbit, 2πρ, so the speed of the motion in an FOC orbit will be¹

$$\begin{eqnarray} v_{{\rm FOC}} &=& 0.0834\left(\frac{D}{1{\rm m}}\right) \left(\frac{1\ \mu {\rm m}}{\lambda }\right) \left(\frac{1{\rm pc}}{d}\right) \nonumber\\ &&\times\sqrt{ \left(\frac{M_*}{1\,M_\odot }\right)\left(\frac{1{\rm AU}}{a}\right)} {\rm in FWHM day}^{-1}. \end{eqnarray} \tag{ 4 }$$

In the general case, the equations of motion in the focal plane are

$$\begin{eqnarray} \dot{x}&=&v_{{\rm FOC}}\sqrt{\frac{1}{1-e^2}}[e\sin (f)\,(\cos (\Omega)\cos (\omega +f) \nonumber \\ &&-\sin (\Omega)\sin (\omega +f)\cos (i)) \nonumber \\ &&-(1+e\cos (f))\,(\cos (\Omega)\sin (\omega +f)\nonumber \\ &&+\sin (\Omega)\cos (\omega +f)\cos (i))] \nonumber \\ \dot{y}&=&v_{{\rm FOC}}\sqrt{\frac{1}{1-e^2}}[e\sin {(f)}\, (\sin {(\Omega)}\cos (\omega +f)\nonumber \nonumber\\ &&+\cos (\Omega)\sin (\omega +f)\cos (i))\nonumber\\ &&-(1+e\cos (f))\,(\sin (\Omega)\sin (\omega +f)\nonumber \\ &&-\cos (\Omega)\cos (\omega +f)\cos (i))] \nonumber \\ v_{{\rm om}} &=& \sqrt{\dot{x}^2+\dot{y}^2}, \end{eqnarray} \tag{ 5 }$$

where Ω is the longitude of the ascending node, ω is the argument of pericenter, i is the inclination, and the true anomaly f depends on a, e, and the time of pericenter passage τ through Kepler's equation (Murray & Correia 2010).

In Figure 1 we show the variation in projected orbital speed for both circular orbits at several inclinations and face-on eccentric orbits (i = 0), for a planet orbiting a 1 M_☉ star at 1 AU. In the plots we normalized speed to 1, and provide v_FOC for several interesting cases. These various scenarios produce projected orbital speeds of appreciable fractions of an FWHM per day. We will later show that, especially for ground-based imaging, this causes a significant degradation in our sensitivity.

Zoom In Zoom Out Reset image size

Our main focus here is on planets in the HZ. Our simple definition of the HZ results in $a_{{\rm HZ}} \propto \sqrt{L_*}$. Now, on the main sequence mass and luminosity approximately follow scaling laws of the form $L_* \propto M_*^b$, where b > 2 except for very massive stars. So according to Equation (4) we expect v_FOC in the HZ to increase as M_* decreases, i.e., M stars will have faster HZ planets than G stars. For example, a planet in the HZ of α Cen B (M_* = 0.9 M_☉, L_* = 0.5 L_☉) will be moving roughly 20% faster than a planet in the HZ of α Cen A (M_* = 1.1 M_☉, L_* = 1.5 L_☉) (stellar parameters from Bruntt et al. 2010).

To provide a more concrete example we return to the 20 hr observation of Gl 581d by the E-ELT/EPICS proposed by Kasper et al. (2010). Using a wavelength of 0.75 μm with Equation (4) we find v_FOC = 0.82 FWHM day⁻¹, or a total of 0.68 FWHM for a continuous 20 hr observation. Since this is a ground-based observation the actual amount of motion to consider is ∼1.15 FWHM over the ∼1.4 days minimum it would take to integrate for 20 hr. Were this a face-on orbit, an eccentricity of 0.25 (Forveille et al. 2011) would increase the maximum orbital speed to as much as 1.05 FWHM day⁻¹, or 1.47 FWHM minimum for a 20 hr ground-based observation.

So what does the orbital motion calculated above do to our observations? To find out we consider a simple model of aperture photometry. Let us assume that we are conducting aperture photometry with a fixed radius r_ap, that the PSF is Gaussian, and that we are limited by Poisson noise from a photon flux N per unit area. With these assumptions, the optimum r_ap is 0.7 FWHM, but taking into account centroiding uncertainty r_ap ≈ 1 FWHM is typical. We will approximate orbital motion at speed v_om by substituting x → x − v_omt − x₀. Orbits are of course not linear, but this will be approximately valid over short periods of time. The parameter x₀ allows us to optimize the placement of the aperture to obtain the maximum signal, i.e., centering the aperture in the planet's smeared-out flux. Note that with the exception of this centering parameter, this model appears quite naive in that we are not adapting the aperture radius and are pretending that we will not notice a smeared-out streak in our images.

Now the S/N in the fixed-size aperture after time Δt will be

$$\begin{eqnarray} &&{\rm S/N}_{{\rm fix}}\nonumber\\ &&\quad = \frac{\int _{0}^{r_{{\rm ap}}} \int _{0}^{2\pi } \int _0^{\Delta t} I_0 e^{\left(-4 \ln 2 ((r\cos \theta - v_{{\rm om}}t - x_0)^2 + r^2\sin ^2\theta)\right)} dt d\theta rdr}{\sqrt{N\pi r_{{\rm ap}}^2 \Delta t_{{\rm int}}}},\nonumber\\ \end{eqnarray} \tag{ 6 }$$

where I₀ is the peak value of the PSF. In the case of no orbital motion v_om = 0 and aperture r_ap = 1 FWHM, so we have

$$\begin{equation} {\rm S/N}_o = \frac{0.6 I_0 \sqrt{\Delta t}}{\sqrt{N}}. \end{equation} \tag{ 7 }$$

As a simple alternative to a fixed size aperture, we also consider allowing our photometric aperture to expand along with the motion of the planet. This aperture will collect the same signal as in S/N_o, but the noise increases with the area as 2r_apv_omΔt, so we have

$$\begin{equation} {\rm S/N}_{{\rm exp}} = \frac{0.6 I_0 \sqrt{\Delta t}}{\sqrt{N\left(1+(2/\pi)v_{{\rm om}}\Delta t\right)}}. \end{equation} \tag{ 8 }$$

A convenient scaling is to multiply top and bottom by $\sqrt{v_{{\rm om}}}$ and work in normalized S/N units of $I_o/\sqrt{Nv_{{\rm om}}}$. This puts time in terms of FWHM of motion, = v_omΔt, and allows comparisons without specifying v_om.

In Figure 2 we plot the normalized S/N versus time (measured in terms of FWHM of motion) with and without orbital motion and for both the fixed and expanding aperture cases. For the fixed aperture, after ∼2 FWHM of orbital motion a maximum of 0.69 is reached, and from there noise is added faster than signal. This means that further integration only degrades the observation.

Zoom In Zoom Out Reset image size

The expanding aperture S/N_exp exceeds the maximum of S/N_fix after about 8 FWHM of motion, and

$$\begin{equation} \displaystyle \lim _{x \rightarrow \infty } 0.6 \frac{\sqrt{x}}{\sqrt{1+(2/\pi)x}} = 0.6\sqrt{\frac{\pi }{2}} \approx 0.75. \end{equation} \tag{ 9 }$$

So if we integrate four times longer, adjusting the aperture size would allow us to gather a little more S/N, but only to a point. Given this large increase in telescope time for a relatively small improvement in S/N (only ∼9% even if we integrate forever), and its better performance for smaller amounts of motion, the fixed-radius aperture will be our baseline for further analysis—keeping in mind that in some cases it may not be the true optimum.

The peak in S/N_fix (Equation (6)) sets the maximum nominal integration time before orbital motion will prevent us from achieving the science goal. That is Δt_max = (S/N_max/0.6)². If the observation of a stationary planet would require an integration time longer than Δt_max, then we can't achieve the desired S/N on an orbiting planet. This also sets the maximum orbital motion _max = v_omΔt_max. From Figure 2 we find that _max = 1.3 FWHM. If more than 1.3 FWHM of motion occurs during an observation, we will not achieve the required S/N.

We also show the fractional reduction in S/N in Figure 2. Almost no degradation occurs until after ∼0.2 FWHM of motion has occurred. S/N is reduced by ∼1% after 0.5 FWHM of motion, ∼5% after 1.0 FWHM, and by ∼19% after 2.0 FWHM of motion. We must now decide how much S/N loss we can accept in our observation.

The above analysis assumes a continuous integration. On a ground-based telescope one must consider that the maximum continuous integration time is ≲ 12 hr, and in practice will likely be much shorter when performing high contrast AO corrected imaging. For instance, an exposure of 20 hr might have to be broken up over 4 or 5 or more nights, when considering the vagaries of seeing (required AO performance), airmass (either through transmission or r₀ requirements), rotation rate (for angular differential imaging (ADI)), and weather. We can adapt the calculations for a ground-based integration as follows:

$$\begin{eqnarray} &&{\rm S/N}_{{\rm gnd}} \nonumber\\ &&\quad = \frac{\int _{0}^{r_{{\rm ap}}} \int _{0}^{2\pi } \left[ \sum _{j=1}^{j=M} \int _{t_j}^{t_j + \Delta t_j} I_0 e^{(-4 \ln 2 ((r\cos \theta - v_{{\rm om}}t - x_0)^2 + r^2\sin ^2\theta))} dt \right] d\theta rdr}{\sqrt{N \pi r_{{\rm ap}}^2 \Delta t_{{\rm int}}}}.\nonumber\\ \end{eqnarray} \tag{ 10 }$$

In this expression we have broken the observation up into M integration sets which start at times t_j and have lengths Δt_j. The total integration time is $\Delta t_{{\rm int}} = \sum _{j=1}^{j=M} \Delta t_j$ and the total elapsed time of the observation is Δt_tot = t_M + Δt_M − t₁.

We plot the results for a few ground-based scenarios in Figure 3. As one can see, observations of planets with orbital motion will be significantly degraded from the ground. This problem, which has been negligible in the high contrast planet searches to date, only becomes worse as we consider larger telescopes and improvements in AO technology which allow searches at shorter wavelengths. We next analyze how this reduction in S/N will affect our ability to detect exoplanets by increasing the rate at which spurious detections occur.

Zoom In Zoom Out Reset image size

Now we turn to the problem of detecting a planet of a given brightness. A planet is considered detected if its flux is above some threshold S/N_t, which is chosen for statistical significance. The goal in choosing this threshold is to detect faint planets while minimizing the number of false alarms. For the purposes of this analysis we assume Gaussian statistics, in which case the false alarm probability (P_FA) per trial is

$$\begin{equation} P_{{\rm FA}} = \frac{1}{2}{\rm erfc}\left(\frac{{\rm S/N}}{\sqrt{2}}\right). \end{equation} \tag{ 11 }$$

Typically, planet hunters use a threshold of S/N = 5, which gives P_FA = 2.9 × 10⁻⁷. The number of false alarms per star, the false alarm rate (FAR), is then

$$\begin{equation} {\rm FAR} = P_{{\rm FA}}\times N_{{\rm trials}}, \end{equation} \tag{ 12 }$$

where N_trials is the number of statistical trials per star. Following Marois et al. (2008a), for a stationary planet N_trials is just the number of photometric apertures in the image. A typical Nyquist sampled detector of size 1024 × 1024 pixels has N_trials ∼ 8 × 10⁴. Thus, an S/N = 5 threshold will result in FAR ∼ 0.02—about 1 false alarm for every 50 observations. In the speckle limited case with non-Gaussian statistics, FAR will be worse than this for the same S/N (Marois et al. 2008a). In any case, the FAR is the statistic which determines the efficiency of a search for exoplanets with direct imaging. A high FAR will cause us to waste telescope time following up spurious detections, while raising the S/N threshold to counter this limits the number of real planets we will detect.

The reduction of S/N caused by orbital motion confronts us with three options. Option I is to maintain the detection threshold constant and accept the loss of sensitivity. Option II is to lower the detection threshold to maintain sensitivity, accepting the increase in FAR. Option III is to correct for orbital motion, which, as we will show, also causes an increase in FAR.

3.3.1. Option I: Do Nothing

The default option is to do nothing, keeping our detection threshold set as if orbital motion is not significant. The drawback to this is that we will detect fewer planets. To quantify this we use the concept of completeness, that is, the fraction of planets of a given brightness we detect. For Gaussian statistics and detection threshold S/N_t = 5, the search completeness is given by

$$\begin{equation} C(\epsilon) = 1-\frac{1}{2}{\rm erfc}\left(\frac{{\rm S/N}(\epsilon) - 5}{\sqrt{2}}\right), \end{equation} \tag{ 13 }$$

where = v_omΔt is the amount of motion. In Figure 4 (top) we show the impact of orbital motion on search completeness. Maintaining the detection threshold lowers completeness. How much depends on the completeness level, with brighter planets being less affected. For planets bright enough to yield 95% completeness with no motion, significant reduction in the number of detections begins after ∼1 FWHM of motion. For 99.7% completeness the impact becomes significant after ∼1.5 FWHM.

Zoom In Zoom Out Reset image size

3.3.2. Option II: Lower Threshold

3.3.3. Option III: De-orbit

Option III is to correct for orbital motion, hoping to maintain sensitivity while limiting the increase in P_FA. The essence of any such technique will be calculating the position of the planet during the observation, and de-orbiting in some way, say shift-and-add (SAA) on a sequence of images. The drawback of this approach is that it will produce more false alarms per observed star due to the increased number of trials, similar to lowering the detection threshold. If the orbit were precisely known, we could proceed with almost no impact on FAR. However, in the presence of uncertainties in orbital parameters or in a completely blind search we will have to consider many trial orbits. For now we can perform a "back-of-the-envelope" estimate of the number of possible orbits to understand how much FAR will increase. To do so, we begin by placing bounds on the problem.

We can first establish where on the detector we must consider orbital motion. At any separation r from the star, the slowest unbound orbit will have the escape velocity. Since we know that physical separation is greater than or equal to projected separation, r ⩾ ρ, and that maximum projected speed will occur for inclination i = 0, we know that

$$\begin{equation} v_{{\rm esc}} = \sqrt{2}v_{{\rm FOC}}(a \rightarrow \rho) \end{equation} \tag{ 14 }$$

sets the upper limit on the projected focal plane speed of an object in a bound orbit. We can also set an upper limit on the amount of motion _max we can tolerate over the duration Δt_tot of the observation based on the S/N degradation it would cause. So we only need consider orbital motion when

$$\begin{equation} \sqrt{2}v_{{\rm FOC}}(\rho)\Delta t_{{\rm tot}} > \epsilon _{{\rm max}}. \end{equation} \tag{ 15 }$$

From here we determine the upper limit on projected separation from the star for considering this problem:

$$\begin{equation} \rho _{{\rm max}} = 0.0136\,M_* \left(\frac{D}{\lambda d} \frac{\Delta t_{{\rm tot}}}{\epsilon _{{\rm max}}}\right)^2 {\rm AU.} \end{equation} \tag{ 16 }$$

By the same logic, for any point closer than ρ_max the maximum possible change in position is

$$\begin{equation} \Delta \rho _{{\rm max}} \approx \sqrt{2}v_{{\rm FOC}}(\rho)\Delta t_{{\rm tot}} {\rm in FWHM}. \end{equation} \tag{ 17 }$$

Then we must evaluate possible orbits ending anywhere in an area of π(Δρ_max)² FWHM² around an initial position.

These two limits set the statistical sensitivity of an attempt to de-orbit an observation. The number of different orbits, N_orb, will be determined by the area of the detector where orbital motion is non-negligible, and the size of the region around each point that we consider. That is,

$$\begin{equation} N_{{\rm orb}} \propto \int _{0}^{\rho _{{\rm max}}} \Delta \rho _{{\rm max}}^2 \rho d\rho, \end{equation} \tag{ 18 }$$

so

$$\begin{equation} N_{{\rm orb}} \propto \left(\frac{M_*}{\epsilon }\right)^2\left(\frac{D}{\lambda d}\right)^4 \Delta t_{{\rm tot}}^4. \end{equation} \tag{ 19 }$$

In general N_trials∝N_orb, so FAR∝P_FA × N_orb. Larger D, shorter λ, closer d, and smaller acceptable orbital motion will then all increase FAR.² Perhaps the most important feature of this result is that $N_{{\rm orb}} \propto \Delta t_{{\rm tot}}^4$—increasing integration time rapidly increases the FAR of a blind search. Note that this is still less severe than the exponential increase in P_FA found for merely lowering the threshold. In the next section we will test these relationships after fully applying orbital mechanics, and see that they hold.

Here we derive limits on the semi-major axis and eccentricity of trial orbits to consider. These limits are based only on the amount of orbital motion tolerable for the science case, and do not represent physical limits on possible orbits around the star.

It is always true that r ⩾ ρ. This implies that, for any orbit, the separation of apocenter must obey r_a ⩾ ρ. This allows us to set a lower bound on a, a_min, given a choice of e through

$$\begin{equation} \rho _{au} \le a_{{\rm min}}(1+e) \end{equation} \tag{ 20 }$$

which gives

$$\begin{equation} a_{{\rm min}} = 0.2063\frac{\lambda d \rho }{D(1+e)}. \end{equation} \tag{ 21 }$$

The fastest speed in a bound planet's orbit will occur at pericenter, and using the maximum tolerable motion _max during our observation of total elapsed time Δt_tot we can set an upper bound on a by noting that

$$\begin{equation} v_{{\rm FOC}}(a_{{\rm max}})\sqrt{\frac{1+e}{1-e}}\Delta t_{{\rm tot}} \le \epsilon _{{\rm max}} \end{equation} \tag{ 22 }$$

which leads to

$$\begin{equation} a_{{\rm max}} = \left(0.0834\frac{D}{\lambda d}\right)^2\frac{1+e}{1-e}M_*\left(\frac{\Delta t_{{\rm tot}}}{\epsilon _{{\rm max}}}\right)^2. \end{equation} \tag{ 23 }$$

Using the GMT example: for e = 0.0, a_max = 3.9 AU; and for e = 0.5, a_max = 11.8 AU. Using Equation (16) we have a projected separation limit of ρ_max = 7.7 AU, so it is possible for these definitions to produce a_max < a_min for certain choices of e at a given ρ. This condition tells us that at such a value of e no orbits can move fast enough to warrant consideration. Thus we can set a lower limit on e at projected separation ρ:

$$\begin{equation} e_{{\rm min}} = \frac{1}{2}\sqrt{\xi ^2+8\xi } - 1 - \frac{\xi }{2}, \end{equation} \tag{ 24 }$$

where we have simplified by pulling out

$$\begin{equation} \xi = 29.66\frac{\rho }{M_*}\left(\frac{\epsilon }{\Delta t}\right)^2\left(\frac{\lambda d}{D}\right)^3. \end{equation} \tag{ 25 }$$

In practice, we might consider eccentricity ranges with e_max less than 1, thus improving our sensitivity. Inputs to our choice of e_max could include some prior distribution of eccentricities, or dynamical stability considerations in binary star systems and systems with known outer companions.

Now we describe an algorithm for sampling the possible trial orbits over a set of M sequential images. For now, we assume no prior knowledge of orbital parameters. We will employ a simple grid search through the parameter space bounded as described above.

• 1.  
  Determine the region around the star to consider using Equation (16).
• 2.  
  Identify regions of interest. In the best cases the orbital motion will be small enough that we will be able to stack the images and search the result for regions with higher S/N (e.g., S/N > 4) and limit further analysis to those areas. In the worst cases orbital motion will be large enough that we will need to blindly apply this algorithm at each pixel within the bounding region identified in the previous step. In the present GMT–αCen example we are in the former case.
• 3.  
  For each region, choose a size, perhaps based on v_esc (as in Equation (17)).
• 4.  
  Choose a starting point (x₁, y₁), with $\rho _1 = \sqrt{\vphantom{A^A}\smash{{{x_1^2 + y_1^2}}}}$. If we are proceeding pixel by pixel, then (x₁, y₁) describes the current pixel.
• 5.  
  Choose e ∈ e_min(ρ₁)...e_max using Equation (24) and assumptions about e_max.
• 6.  
  Choose a ∈ a_min(ρ₁, e)...a_max(e) using Equations (21) and (23).
• 7.  
  Choose time of pericenter τ ∈ t₁ − P(M_*, a)...t₁ where P is the orbital period and t₁ is the time of the first image. Now calculate the true anomaly f(t₁; a, e, τ, P) using Kepler's equation and physical separation using

  $$\begin{equation} r = \frac{a(1-e)}{1+e\cos (f)}. \end{equation} \tag{ 26 }$$
• 8.  
  if e ≠ 0: Choose ω ∈ 0...2πif e = 0: set ω = 0.
• 9.  
  if sin (ω + f) > 0:
• 10.  
  if sin (ω + f) = 0, we do not have a unique solution for inclination. This is the special case where the planet is passing through the plane of the sky.

  • (a)  
    for ω + f = 0 calculate Ω:

    $$\begin{equation} \sin \Omega = \frac{y}{r} \end{equation} \tag{ 30 }$$

    $$\begin{equation} \cos \Omega = \frac{x}{r} \end{equation} \tag{ 31 }$$

    or for ω + f = π calculate Ω:

    $$\begin{equation} \sin \Omega = \frac{-y}{r} \end{equation} \tag{ 32 }$$

    $$\begin{equation} \cos \Omega = \frac{-x}{r}, \end{equation} \tag{ 33 }$$

    determining Ω in the correct quadrant.
  • (b)  
    Choose i ∈ 0...π.
  • (c)  
    We now have a complete set of elements, and so can SAA as in step 9b above.
  • (d)  
    Repeat steps 10b to 10c until all i chosen.
• 11.  
  Repeat the above steps until the parameters ω, τ, a, and e are sufficiently sampled for each starting point.

The algorithm just described will produce a large number of trial orbits, many of which will be very similar. The information content of our image is set by the resolution of the telescope, so we can take advantage of this similarity to greatly reduce the number of statistical trials. This is done by grouping similar orbits into sequences of whole-pixel shift sequences, where the pixels are at least as small as FWHM/2. As we will see, we typically will want to oversample, to say FWHM/3, to ensure adequate S/N recovery.

We calculate the pixel-shift sequence for each orbit by determining which pixel the trial planet (or rather, the center of its PSF) lands on at each time step. Many orbits end up producing the same sequences of pixel-shifts, and we will keep only the unique ones for use in de-orbiting the observation. In Figure 5 we illustrate the outcome of the pixel-shift algorithm, showing two unique sequences and a few of the orbits that produced them.

Zoom In Zoom Out Reset image size

To test the above algorithm and the pixel-shift technique, we used our GMT α Cen A example and determined the trial orbits for various separations and Δt values. We set _max = 0.5 based on our earlier analysis of S/N. The results are summarized in Figure 6. The problem is generally well constrained in that we only have a finite search space for any initial point. The data used to construct Figure 6 are provided in Table 1. Comparing N_orb to N_shifts, note the large reduction in the number of trials (∼10⁸ to ∼10²) due to combining similar orbits.

Zoom In Zoom Out Reset image size

In Figure 7 we plot the area of the detector which contains the possible trial orbits at ρ₁ = 1.0 AU versus the total elapsed time Δt_tot. We conclude from this plot that the area around a given starting point is proportional to $\Delta t_{{\rm tot}}^2$. Also in Figure 7 we plot area versus separation from the star, and conclude that area is proportional to 1/ρ₁. Taken together these results give confidence that the $N_{{\rm orb}} \propto \Delta t_{{\rm tot}}^4$ scaling derived earlier holds when we fully apply orbital mechanics rather than the escape velocity approximation.

Zoom In Zoom Out Reset image size

Things are a bit more complicated when we consider the scaling of the number of whole-pixel shift sequences. We conducted two sets of trials at ρ₁ = 1.0 AU. In the first, the number of observations and their relative spacing were held constant regardless of Δt_tot. In the second set, the number of observations scaled with Δt_tot. As shown in Figure 7, when the number of observations is constant, the number of shifts scales as $\Delta t_{{\rm tot}}^2$, but when the number of observations grows with Δt_tot, the number of shifts scales as roughly $\Delta t_{{\rm tot}}^{3.6}$. Figure 7(d) shows that the number of shifts scales as 1/ρ₁. Taken together, we see that for a constant number of observations the pixel-shift technique will follow the $N_{{\rm orb}} \propto \Delta t_{{\rm tot}}^4$ scaling. However, if the number of observations also scales with Δt_tot, then our results imply that $N_{{\rm orb}} \propto \Delta t_{{\rm tot}}^{5.6}$. The value of the exponent likely depends on the details of the observation sequence, but this has important implications for observation planning.

We next consider whether de-orbiting by whole-pixels adequately recovers S/N. To test this we "orbited" a Gaussian PSF on face-on orbits with various eccentricities, starting from pericenter. We then calculated shifts for detector samplings of 2, 3, and 4 pixels/FWHM, and then de-orbited by these shifts. The results are summarized in Table 2. On a critically sampled detector we only recover a 5σ planet to ∼4.9σ, a 2% loss of S/N. At 3 pixels/FWHM we do much better, recovering S/N to 4.97 for low eccentricities, and 4.95 for higher eccentricities. Performance for 4 pixels/FHWM sampling is similar. A 2% loss of S/N nearly doubles P_FA, so it appears that we should oversample to at least 3 pixels/FWHM, either optically or by re-sampling images during data reduction. In our analysis we have assumed that the limiting noise source is background photons (PSF halo or sky), so we ignore the increased readout noise expected from oversampling.

As we have noted several times, the main impact of orbital motion is to reduce S/N, which in turn reduces our statistical sensitivity. If we attempt to de-orbit an observation in order to recover S/N, we do so at the cost of a large increase in the number of trials. Worst case, this results in a proportional increase in FAR since nominally FAR = P_FA × N_orb. However, we expect significant correlation between trials of neighboring orbits and whole-pixel shifts. To investigate this, we performed a series of Monte Carlo (MC) experiments. A sequence of images with Gaussian noise was generated, and first stacked without shifting, hereafter called the naive-add. The same sequence was then shifted by each possible whole-pixel shift, assuming a 1 AU initial separation around α Cen A. This experiment was conducted for observations with total elapsed times Δt_tot of 4.2, 6.2, 8.2, and 10.2 days, with samplings of 2, 3, and 4 pixels/FWHM.

We performed several tests on each sequence. The first was a simple threshold test on the naive-add, with the threshold set for the worst-case orbital motion given by Equation (10) with v_om = v_esc. We performed simple aperture photometry, with an r_ap = 1 FWHM. As expected the resultant P_FA1 is as predicted by Equation (11). The next test was to apply a 5σ threshold after de-orbiting by whole-pixel shifts and adding. If all shifts were completely uncorrelated, then we would expect FAR = (2.9 × 10⁻⁷) × N_shifts, but as we predicted, shifts are correlated and P_FA2 is lower than this.

The final test performed was to apply both thresholds in sequence, such that a detection is made only if the naive-add results in S/N greater than the threshold for worst-case orbital motion, and the de-orbited SAA results in S/N > 5. This P_FA3 is lower than either P_FA1 or P_FA2, but still higher than if no orbital motion occurred.

The results of each trial are present in Table 3. Applying both threshold tests results in significant improvement over the naive-add in terms of FAR. Another interesting result is that sampling has only a minor impact on P_FA3. This makes some sense as we expect the correlation of neighboring shifts to be set by the FWHM, not the sampling. So even though the accuracy of S/N recovery is improved, and quite a few more shifts are required, these shifts remain correlated across the same spatial scale, resulting in little change in the overall FAR.

There is still an impact on completeness, however, because we are now conducting two trials instead of one. This lowers the true positive probability (P_TP). Consider a 5σ planet on the worst-case fastest possible orbit, for the 10.2 day elapsed time case. The threshold for the naive-add is 2.625. We have a 50% probability of detecting this planet after the naive-add. If it is detected on the first test, there is then some probability P_TP < 1 of detecting at S/N ⩾ 5 after de-orbiting. Worst case, this will be 50%, resulting in a net P_TP of 25%. In reality, it will be better than this as the two trials will be strongly correlated.

Even if this worst case of 25% were realized, this is still significant improvement over Option I. A 2.6σ signal would only be detected 10% of the time with a 5σ threshold. Given the reduction in P_FA from 4.3 × 10⁻³ to 2.2 × 10⁻⁵, likewise an improvement over Option II at 2.6σ, it is clear that de-orbiting by whole-pixel shifts does improve our ability to detect an orbiting planet. The situation will be even better for slower planets, and most of the area searched will not be subject to the worst-case orbital speed. We leave a complete analysis of the impact on search completeness for future work. One can also imagine adjusting the thresholds to optimize completeness at the expense of worse P_FA.

We end this section by concluding that a blind search when orbital motion is significant is tractable. Orbital motion will make such a search less sensitive, both in terms of number of false alarms and in terms of completeness, but Keplerian mechanics gives us enough tools to bound the problem. As we have shown, de-orbiting a sequence of observations can recover S/N to its nominal value, and we can do so while controlling the impact on statistical sensitivity. For the Δt_tot = 6.2 day observation, P_FA3 was roughly a factor of 10 higher than if no orbital motion occurred. This increase only occurs over a bounded region around the star, so the net effect on FAR will be contained. Using this factor of 10 as the mean value over the 7.7 AU = 69.1 FWHM radius region around α Cen A where orbital motion is significant, the FAR in this area will have gone from ∼1/1000 to ∼1/100 in our GMT/10 μm example. The key, though, appears to be to limit the elapsed time of the observation as the number of trials increases—decreasing sensitivity—proportionally to at least $\Delta t_{{\rm tot}}^4$ in a blind search.

The main caveat at this point in our analysis is that we have drawn the conclusion of tractability using Gaussian statistics. It is well known that speckle noise, which will often be the limiting noise source for high contrast imaging in the HZ, is not Gaussian and results in much higher P_FA for a given S/N (Marois et al. 2008a). Future work on this problem will need to take this into account.

Next we consider a more strongly bounded scenario, where we have significant prior information about the orbit of the planet from RV surveys.

In order to minimize the number of trial orbits to consider, we can apply various constraints taking advantage of the information we have from the RV detection.

In the case of a multi-planet system dynamical analysis can place constraints on the inclination based on system stability. For Gl 581, Mayor et al. (2009) found the system was stable for i > 30. We can also make use of the geometric prior for inclination, where we expect P_i = sin (i) in a population of randomly oriented systems.

Since this is a reflected light observation, the orbital phase and its impact on the brightness of the planet must be considered. The planet's reflected flux is given by

$$\begin{equation} F_p(\alpha) = F_*\left(\frac{R_p}{r}\right)^2 A_g(\lambda) \Phi (\alpha), \end{equation} \tag{ 34 }$$

where F_* is the stellar flux, R_p is the planet's radius, r is its separation, A_g(λ) is the wavelength-dependent geometric albedo, and Φ is the phase function at phase angle α. The phase angle is given by

$$\begin{equation} \cos (\alpha) = \sin (f + \omega)\sin (i). \end{equation} \tag{ 35 }$$

In general, determining the quantity A_g(λ)Φ(α) requires atmospheric modeling (Cahoy et al. 2010). For now, we assume that Φ follows the Lambert phase function

$$\begin{equation} \Phi (\alpha) = \frac{1}{\pi }\left[\sin (\alpha) + (\pi - \alpha)\cos (\alpha)\right]. \end{equation} \tag{ 36 }$$

We assume that the prediction of Kasper et al. (2010) was made for the planet at quadrature, α = π/2, where Φ = 0.318. We then require that the mean value of Φ during the observation be greater than this value—that is, the planet is as bright or brighter than it is at quadrature.

An important consideration in an RV-cued observation will be when to begin. As a first approximation, we assume that maximizing planet–star separation will maximize our sensitivity. This may not be true when working in reflected light due to the phase and separation dependent brightness of the planet in this regime. Proceeding with the approximation for now, we expect to plan this observation to be as close to apocenter as possible. In this case we will begin integrating 3.1 days before t₀ + P/2.

To understand the area where we will be searching for Gl 581d, we first conducted an MC experiment to calculate the possible positions of the planet at t = t₀ + P/2 − 3.1 days. To do so, we drew random values of a, e, w, and t₀ from Gaussian distributions with the parameters of Table 4. We drew a random value of i from the sin (i) distribution, and rejected any value of i ⩽ 30 based on the dynamical prior. Finally Ω was drawn from a uniform distribution in 0...2π. This process was repeated 10⁹ times, and the frequency at which starting points occur in the area around the star was recorded. The results are shown in Figure 8 for the V12 circular model and for the F11 eccentric model. The figure shows the area which must be searched to obtain various completeness. For instance, if we desire 95% completeness in the V12 model, we must consider an area of 71 apertures. Since this S/N = 5 detection is broken up into five distinct integrations, our first attempt will have S/N = 2.24, giving an FAR = 0.89 for the first 4 hr integration. In other words, we should expect a false alarm in addition to a real detection.

Zoom In Zoom Out Reset image size

Now we assume that we have an initial detection at S/N ∼ 2.24 within the highest probability regions.³ In order to follow up this detection over subsequent nights, we must determine the possible locations of the planet, constrained by the RV-derived orbital elements.

We proceed by choosing a, e, ω, and t₀ from Gaussian distributions as above. Now as long as r > ρ we will have a unique solution for i and Ω given the randomly chosen parameters (see the blind search algorithm above). We take into account dynamical stability by rejecting any orbit which has i ⩽ 30. The orbit determined in this fashion was then projected 6.2 days into the future and the frequency of these final points was recorded. We show the result for the V11 model in Figure 9, top panel. Using the RV determined parameters and their uncertainties allows us to determine the probability density of orbit endpoints, and determine how much of the search space we must consider for a given completeness. The whole-pixel shifts were also calculated using a sampling of FWHM/3, and are shown in the legend. We also applied the blind search algorithm to this observation from the same starting point, and show the results for comparison in the bottom panel of Figure 9. As expected the RV priors significantly reduce the search space—we have 942 trial shift-sequences to consider instead of 12,000.

Zoom In Zoom Out Reset image size

Another important consideration here is that our initial 2.24σ detection will have a large position uncertainty, which we estimate by $\sigma _{\rho _0} = {\rm FWHM}/{\rm S/N}$. We added a random draw for the starting position, and repeated the MC experiment for F11 and also conducted a run for the V12 parameters. The results are shown in Figure 10. The number of shift sequences is much higher due to the uncertainty in the starting position caused by our low S/N initial detection, but we expect correlations to come to the rescue as in our α Cen example. To compare to Figure 9 keep in mind that the blind search would have to be applied to all 5500 pixels in the search space indicated by Figure 8.

Zoom In Zoom Out Reset image size

As in the GMT/α Cen example, we leave for future work a complete analysis of sensitivity and completeness. The large number of trial shifts calculated when we include uncertainty in the starting position motivates us to suggest that we will ultimately turn this analysis over to a much more robust optimization strategy, such as a Markov Chain Monte Carlo (MCMC) routine. Once an area of the image was identified with a high post-shift S/N, an MCMC analysis could determine the very best orbit and assign robust measures of significance to the result.

We also note that these results likely overestimate the number of trial orbits since we have assumed uncorrelated errors. In reality the RV best-fit parameters are likely strongly correlated, which should act to reduce the number of orbits to consider.