DIRECT EXOPLANET DETECTION WITH BINARY DIFFERENTIAL IMAGING^*

A fundamental limit to directly imaging planets on solar system scales is the "speckle floor." This is the region extending out to ∼15 from the star that is dominated by speckles: time-varying aberrations in the point spread function (PSF). Because the locations and brightnesses of the speckles change slowly with time, it is difficult to fully subtract them, even with the help of adaptive optics (AO). Speckle noise thus dominates over photon and sky noise close to the star (Racine et al. 1999), limiting our ability to directly detect faint planets at small separations.

Angular differential imaging (ADI; Marois et al. 2006), in which the target star serves as its own PSF reference, has helped improve contrast close to the star, but only to a certain level and usually at the cost of attenuating the fluxes of the off-axis companions. ADI also requires a sufficient amount of sky rotation to occur throughout the observations: the less the sky rotates, the more a companion's flux at a given separation is attenuated. This means that planets located in the speckle floor region will be more attenuated than those farther out, making their detection even more challenging.

Simultaneous differential imaging (SDI; Racine et al. 1999; Marois et al. 2003) can in theory completely remove speckles because the reference PSF is the target itself imaged simultaneously at different wavelengths, but it requires the faint companions to have unique and prominent spectral features over narrow wavelength regions compared to the stars—an assumption that is so far uncertain. Furthermore, the different Strehl ratios of the two simultaneously acquired images preclude true photon noise limited PSF subtraction. For these reasons, it has been difficult for current and past direct imaging surveys to detect low-mass planets at small separations.

At wider separations, a handful of planets have been directly imaged (e.g., HR 8799 bcde, Marois et al. 2008, 2010; Beta Pic b, Lagrange et al. 2009, 2010; HD 95086 b, Rameau et al. 2013). However, with such a small sample size, it is difficult to make statistical inferences on the general properties of giant planets, and even harder to constrain formation and evolution theories. To alleviate this problem, the number of imaged exoplanets must increase.

We know from several imaging surveys that giant planets are rarely found at wide separations (Nielsen & Close 2010; Janson et al. 2012; Biller et al. 2013; Wahhaj et al. 2013; Maire et al. 2014). We also know from radial velocity and transit surveys that giant planets do reside at small separations (semimajor axis a < 5 AU; Butler et al. 2006; Cumming et al. 2008; Borucki et al. 2011). Therefore future imaging surveys must improve detection sensitivity close to the star.

Detecting planets at close separations—without building new telescopes/instrumentation—calls for a different PSF subtraction technique. To avoid self-subtraction (which plagues ADI), the target star cannot be used as its own PSF reference. To avoid atmospheric requirements in the companion (which plague SDI), the target and PSF reference must be imaged at the same wavelength. "Classical" PSF subtraction, wherein observations of a reference PSF are interleaved with observations of the target, is inefficient (two or more stars imaged to search for planets around one) and suffers from a time-varying PSF. Therefore the only remaining option is to observe two (or more) stars at the same time at the same wavelength.

This is not a new idea. Speckle holography has been demonstrated as a viable technique to infer the presence of faint companions using the PSFs of multiple simultaneously imaged stars (Lohmann & Weigelt 1978, pp. 479–494; Weigelt 1978; Koresko 2002), as long as they are within the isoplanatic patch at the observed wavelength (Christou & Hege 1986). In space, the Hubble Space Telescope has been used to search for circumstellar dust in a binary by using each component's PSF (Low et al. 1999). From the ground with AO, McElwain (2007) used Keck/OSIRIS to search for companions in several close visual binaries (no companions were found).

We can improve on these previous studies in three ways. First, we can take advantage of modern AO systems, which now offer high Strehl ratio imaging, especially at long wavelengths. This improves our ability to detect faint planets close to their host stars. Second, we can image near 4 μm. At this wavelength, gas giant planets are thought to be bright (Baraffe et al. 2003; Burrows et al. 2003), and the isoplanatic patch is large (∼10''–30''). A larger isoplanatic patch means binaries with larger projected separations will have nearly identical PSFs. Indeed, Heinze et al. (2010) imaged several binaries at 4–5 μm with the MMT deformable secondary AO system and showed that subtracting each PSF from the other could be more effective than ADI. This is possible because the target is no longer used as its own PSF reference, so the flux from the companion will never be attenuated by self-subtraction. Furthermore because the PSF is imaged at the same wavelength as the target, the companion is not required to have unique spectral features relative to the star (as with SDI).

Finally, we can take advantage of advanced post-processing algorithms like principal component analysis (PCA, Soummer et al. 2012), which has been shown to significantly improve PSF subtraction. This raises the possibility of completely removing the speckle floor from each star during the data reduction stage and thus reaching the photon noise limit. The photon noise limit is a critical threshold because, once reached, longer integrations result in increased sensitivity. This is not the case with ADI, where at small separations speckle noise dominates and adequate sky rotation, rather than integration, is preferred.

It is now clear that planets can and do form in binaries (Buchhave et al. 2011; Ngo et al. 2015). With the added benefit of potentially superior PSF subtraction, visual binaries should no longer be excluded from direct imaging planet surveys. In this paper, we show that Binary Differential Imaging⁶ (BDI), in which two stars are imaged simultaneously, at high Strehl ratio, at the same wavelength, and within the isoplanatic patch, both improves detection sensitivity at small separations by ∼ half a magnitude and is 2–4× more observationally efficient than current ground- and space-based direct imaging surveys, respectively. In Section 2, we describe the BDI technique. In Section 3, we present proof of concept data obtained at the Magellan Clay telescope using MagAO (Close et al. 2010) and the near-infrared (NIR) camera Clio-2 (Sivanandam et al. 2006). In Section 4 we compare the achieved contrast of ADI and BDI. In Section 5 we discuss the technical and scientific benefits of BDI, on both ground- and space-based telescopes. In Section 6 we summarize and conclude.

Two stars are simultaneously imaged on a detector at the same wavelength with AO. The separation of the stars is within the isoplanatic patch at that wavelength. Star A is brighter than Star B by $f={10}^{\frac{{M}_{A}-{M}_{B}}{-2.5}},$ where M is the magnitude of a given star at the chosen wavelength. Consider two pixels (Pixel A and Pixel B) each having the same coordinates (radius, position angle) relative to each star. Assume that a hypothetical planet resides around Star A and some of its flux falls in Pixel A; Star B has no planets. We will calculate the total signal-to-noise ratio (S/N) of the planet's flux that results from subtracting one pixel from the other. In real high-contrast imaging PSF subtraction, many pixels comprise an image that is subtracted from another image, but the case is more simply explained using a single pixel.

The total flux in Pixel A is given by

$$\begin{eqnarray}&&{T}_{A}={I}_{A}+{I}_{P}+S,\end{eqnarray} \tag{ 1 }$$

where I_A is the flux from Star A, I_P is the flux from the planet, and S is the flux from the sky. The sky flux is removed by subtracting blank sky (S₂) so that ${T}_{A}\to {I}_{A}+{I}_{P}+S-{S}_{2}.$ We will assume that the sky flux is completely removed so that the total flux in Pixel A is just the flux from the star and planet. The same is true for Pixel B (without planet flux). In both cases, however, both read noise (R) and photon noise from the star (σ_I) and sky (σ_sky) contribute to the uncertainty. We will assume that the photon noise from the planet (${\sigma }_{{I}_{P}}$) is negligible. Therefore the total noise in Pixel A is given by:

$$\begin{eqnarray}{\sigma }_{{T}_{A}} & = & \sqrt{{\displaystyle \frac{\partial {T}_{A}}{\partial {I}_{A}}}^{2}{\sigma }_{{I}_{A}}^{2}+{\displaystyle \frac{\partial {T}_{A}}{\partial S}}^{2}{\sigma }_{S}^{2}+{\displaystyle \frac{\partial {T}_{A}}{\partial {S}_{2}}}^{2}{\sigma }_{{S}_{2}}^{2}+2{R}^{2}}\\ & = & \sqrt{{\sigma }_{{I}_{A}}^{2}+2{\sigma }_{\mathrm{sky}}^{2}+2{R}^{2}}\\ & = & \sqrt{{I}_{A}+2S+2{R}^{2}},\end{eqnarray} \tag{ 2 }$$

where the stellar photon noise ${\sigma }_{{I}_{A}}=\sqrt{{I}_{A}},$ and we have assumed that ${\sigma }_{S}={\sigma }_{{S}_{2}}={\sigma }_{\mathrm{sky}}=\sqrt{S}$ and that the read noise is the same in the target and sky images. The total noise in Pixel B is similarly given by

$$\begin{eqnarray}{\sigma }_{{T}_{B}} & = & \sqrt{{I}_{B}+2S+2{R}^{2}}\\ & = & \sqrt{{I}_{A}/f+2S+2{R}^{2}},\end{eqnarray} \tag{ 3 }$$

where we have made the substitution ${I}_{B}={I}_{A}/f.$

Scaling and subtracting the Pixel B flux from Pixel A, we then have:

$$\begin{eqnarray}&&F={T}_{A}-{{fT}}_{B},\end{eqnarray} \tag{ 4 }$$

which = I_P if the PSFs are identical. However, this is an ideal case and is unrealistic. Instead, there will probably be some slight differences in the PSFs due to non-isoplanatism. Therefore the remaining flux will be F = I_P + i, where i is this (positive or negative) extra flux. We can define a new variable $G=F-i={T}_{A}-{{fT}}_{B}-i$ and calculate the total noise:

$$\begin{eqnarray}{\sigma }_{\mathrm{BDI}} & = & \sqrt{{\displaystyle \frac{\partial G}{\partial {T}_{A}}}^{2}{\sigma }_{{T}_{A}}^{2}+{\displaystyle \frac{\partial G}{\partial {T}_{B}}}^{2}{\sigma }_{{T}_{B}}^{2}+{\displaystyle \frac{\partial G}{\partial i}}^{2}{\sigma }_{i}^{2}}\\ & = & \sqrt{{I}_{A}+2S+2{R}^{2}+{f}^{2}({I}_{A}/f+2S+2{R}^{2})+{\sigma }_{i}^{2}}\\ & = & \sqrt{(1+f){I}_{A}+(1+{f}^{2})(2S+2{R}^{2})+{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 5 }$$

where we have assumed that the uncertainty in the scaling factor σ_f ≈ 0. In reality, this uncertainty will be non-zero, but it is typically much smaller than the other sources of noise and is therefore ignored here.

The total S/N of the planet in the BDI case is given by

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{\mathrm{BDI}}=\displaystyle \frac{c\;{I}_{A}}{\sqrt{(1+f){I}_{A}+(1+{f}^{2})(2S+2{R}^{2})+{\sigma }_{i}^{2}}}.\end{eqnarray} \tag{ 6 }$$

Here we have made the substitution ${I}_{P}=c\;{I}_{A},$ where c is the contrast of the planet relative to the host star at the location of the planet.

It is thus evident that as f increases beyond the ideal f = 1, S/N_BDI decreases. In the case of Pixel B minus Pixel A, we again have Equation (6) except ${I}_{A}\to {I}_{B}$ and $f\to 1/f.$ In this case, S/N_BDI is maximized as $f\to \infty .$ If the two PSFs are truly identical, then σ_i = 0 and the noise in BDI PSF subtraction is limited by the brightness ratio of the two imaged stars, with better PSF subtraction (in terms of noise) occurring when the brighter star PSF is subtracted from the fainter star PSF. If the PSFs are not identical, then the noise is limited by the amount of non-isoplanatism (σ_i).

2.1. Comparison to ADI

In the case of ADI PSF subtraction, Pixel B and Pixel A both correspond to Star A, but Pixel B is recorded on the detector some time after Pixel A. We will again assume that Pixel A contains the planet flux, and that Pixel B contains no planet flux. However, to account for the self-subtraction that is inherent to ADI, we will assume that ${I}_{P}\to a(r){I}_{P},$ where a is the attenuation factor that depends on the distance from the star r.

In theory, the brightness ratio of Pixel A and B (f) is unity, but the PSFs will never be identical. Therefore the noise is limited by the differences in the PSFs. To illustrate, we can derive σ_ADI in the case where we have the same star imaged at two different times resulting in two different total fluxes, T_A (which contains the planet flux) and T_B. The subtraction of the two pixels gives $F={T}_{A}-{T}_{B}=a(r){I}_{P}+k,$ where k is the residual flux that results from the non-identical PSF subtraction. As before we can define a new variable $H=F-k={T}_{A}-{T}_{B}-k$ and calculate the noise:

$$\begin{eqnarray}{\sigma }_{\mathrm{ADI}} & = & \sqrt{{\displaystyle \frac{\partial H}{\partial {T}_{A}}}^{2}{\sigma }_{{T}_{A}}^{2}+{\displaystyle \frac{\partial H}{\partial {T}_{B}}}^{2}{\sigma }_{{T}_{B}}^{2}+{\displaystyle \frac{\partial H}{\partial k}}^{2}{\sigma }_{k}^{2}}\\ & = & \sqrt{2{I}_{A}+2(2S+2{R}^{2})+{\sigma }_{k}^{2}}.\end{eqnarray} \tag{ 7 }$$

The total S/N of the planet in the ADI case is given by

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{\mathrm{ADI}}=\displaystyle \frac{a(r)\;c\;{I}_{A}}{\sqrt{2{I}_{A}+2(2S+2{R}^{2})+{\sigma }_{k}^{2}}}.\end{eqnarray} \tag{ 8 }$$

Equation (7) is identical to the ideal (f = 1) BDI noise (Equation (5)), except for the different noise term σ_k (compared to σ_i). This term (speckle noise) can in practice be estimated by calculating the standard deviation of pixels in an annulus at a given radius; it typically increases closer to the star and can be orders of magnitude larger than the other noise sources (Racine et al. 1999). Thus ADI will yield lower noise than BDI if σ_k < σ_i. However, S/N_ADI is also limited by the attenuation factor a(r). Figure 1 shows the theoretical S/N for BDI (ideal, f = 1) and ADI, for different values of σ_k, σ_i, and a(r). It will be our goal in this paper to determine with real data how these parameters affect the achieved S/N of each imaging technique.

Zoom In Zoom Out Reset image size

2.2. Comparison to Classical PSF Subtraction

We can also compare BDI PSF subtraction to classical reference star PSF subtraction (sometimes referred to as "Reference Differential Imaging," or RDI), which is typically employed for space-based telescopes. We will assume we are observing with a coronagraph on the the James Webb Space Telescope (JWST) NIRCam instrument. We will not compare to PSF "roll" subtraction (Schneider & Silverstone 2003) because JWST's limited roll angle (∼10°) will likely preclude exoplanet imaging using this technique, at least for close-in planets.

The noise due to BDI PSF subtraction with NIRCam (no coronagraph⁷ ) is again given by Equation (5). For simplicity, we will assume the ideal case of two identical binary stars (f = 1), no sky or read noise, no noise due to the instrument (i.e., distortion), and no isoplanatism effects σ_i = 0, so that ${\sigma }_{\mathrm{BDI}}=\sqrt{2{I}_{A}}.$ The signal in the BDI case is simply cI_A, so that:

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{\mathrm{BDI},\mathrm{space}}=\displaystyle \frac{c\;{I}_{A}}{\sqrt{2{I}_{A}}}.\end{eqnarray} \tag{ 9 }$$

The noise due to coronagraphic PSF subtraction is similar to the BDI case (again assuming an ideal f = 1), except for the extra noise due to PSF subtraction residuals, σ_k. Therefore we have: ${\sigma }_{\mathrm{coron}}=\sqrt{2{{gI}}_{A}+{g}^{2}{\sigma }_{k}^{2}},$ where g is the fractional decrease in flux that is caused by the coronagraph. The signal of the planet is ${{tcI}}_{A},$ where t is the throughput of the coronagraph. The total S/N is then:

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{\mathrm{coron},\mathrm{space}}=\displaystyle \frac{t\;c\;{I}_{A}}{\sqrt{2{{gI}}_{A}+{g}^{2}{\sigma }_{k}^{2}}}.\end{eqnarray} \tag{ 10 }$$

Comparing the S/N terms, we can see that ${\rm{S}}/{{\rm{N}}}_{\mathrm{BDI},\mathrm{space}}\;\gt {\rm{S}}/{{\rm{N}}}_{\mathrm{coron},\mathrm{space}}$ when ${\sigma }_{k}\gt {g}^{-1}\sqrt{2{I}_{A}({t}^{2}-g)}.$ This relation can also be expressed as ${\sigma }_{k}\gt {g}^{-1}\sqrt{2({t}^{2}-g)}{\sigma }_{{I}_{A}},$ where ${\sigma }_{{I}_{A}}$ is the stellar photon noise. For g ≈ 10⁻¹ and t = 0.81 (Krist et al. 2007), the S/N from BDI PSF subtraction is larger if ${\sigma }_{k}\gtrsim 10{\sigma }_{{I}_{A}}.$ It is not yet clear what σ_k will be for JWST/NIRCam, but this threshold is a useful diagnostic. Even if σ_k is small (a few× the stellar photon noise), BDI may still be preferred because it can reach better IWAs than coronagraphic imaging at a minor cost of slightly larger noise. Furthermore, because BDI is essentially photon noise limited, with enough integration it can always supersede classical PSF subtraction.

To test the feasibility of BDI—in particular, the similarity of two PSFs separated by several arcseconds on the sky—we carried out proof of concept observations of a binary with MagAO/Clio-2. We observed the young visual binary HD 37551, which consists of a G7 (primary) and K1 (Torres et al. 2006) separated by ∼4'' at a distance of 77 ± 9 pc (van Leeuwen 2007). The primary (henceforth Star A) has a V magnitude = 9.650 ± 0.017 and the secondary (henceforth Star B) has V = 10.538 ± 0.017 (Lampens et al. 2001). HD 37551 is thought to be a member of the AB Dor moving group (Torres et al. 2006; Malo et al. 2013). However, the Li depletion pattern of the low-mass AB Dor members and their color–magnitude diagram positions appear to be indistinguishable with the Pleiades (Luhman et al. 2005; Barenfeld et al. 2013), whose Li-depletion boundary age is 130 ± 20 Myr (Barrado y Navascués et al. 2004). We therefore adopt this age for HD 37551.

We observed the binary on the night of UT 2014 November 26 at the Magellan Clay telescope of the Las Campanas Observatory in Chile. We obtained diffraction-limited AO images (300 modes) of the binary in the narrowband 3.9 μm filter (Δ λ = 0.23 μm, λ_c = 3.95 μm; hereafter [3.9]) using the narrow camera (plate scale = 001585 per pixel⁸ ). Star A served as the guide star for the AO system. The instrument rotator was turned off to enable ADI; therefore Star B rotated with the sky around Star A. Observing conditions were good, with the seeing ∼075. Every few minutes, we nodded the binary off the detector to acquire images of blank sky. Neither star's PSF was saturated in any of the recorded images.

After discarding frames of poor quality (e.g., open loop), we acquired 189 images of the binary for a total of 31.5 minutes of integration. The parallactic angle (PA) changed by 1937 throughout the observations. While obviously not a large amount of sky rotation, this is sufficient for directly comparing the detection sensitivities of ADI and BDI.

All data reduction was performed using custom scripts in Matlab. We divided the images by the number of coadds, corrected them for nonlinearity, and divided them by the integration time to obtain units of detector counts/s. To remove the sky background, for each image of the binary we constructed an optimal sky image using all the sky images, such that the noise in a small rectangular box far from the stars was minimized. Next we determined the sub-pixel location of each star by calculating the center of light inside a 025 aperture centered on the stars. We then registered and cropped the images of each star. Finally, we corrected for bad pixels. Figure 2 shows the median PSFs for Star A and Star B, respectively. As expected, the PSFs (and speckles) are nearly identical, despite the different stellar brightnesses and 4'' separation.

Zoom In Zoom Out Reset image size

Next, we performed ADI and BDI PSF subtraction. To simplify the comparison of the two techniques, we performed PSF subtraction on only Star A. Note that according to Section 2, the PSF subtraction noise for BDI will be higher in this case, so our BDI results should be viewed as conservative.

3.1. ADI Reduction

3.2. BDI Reduction

To fully test the capabilities of BDI, we computed the overall contrast limits achieved by ADI and BDI. We accomplished this using a Monte Carlo approach similar to the one employed in Rodigas et al. (2015) whereby a single artificial planet with a random contrast (Δm) at a random separation (r) and position angle (θ) is repeatedly inserted into the raw data and recovered using both ADI and BDI. We inserted and recovered a total of 10,000 artificial planets, which were flux-scaled replicas of the Star A PSF. The allowed parameter ranges were $r\in [0\buildrel{\prime\prime}\over{.} 15,1\buildrel{\prime\prime}\over{.} 55]$, $\theta \in [0,2\pi )$, ${\rm{\Delta }}m\in [5,9]$. Figure 3 shows examples of artificial planets recovered by both BDI and ADI.

Zoom In Zoom Out Reset image size

We computed the S/Ns of the recovered artificial planets using the updated definition of S/N from Mawet et al. (2014), which is appropriate for speckle-dominated regions close to the star. Figure 4(a) shows the separation and contrast values of planets that were detected with S/N ≥ 5 by ADI and BDI. BDI generally detects planets closer to the star than ADI, being ∼0.5 mag better in terms of contrast within ∼1''. Another way of showing this is using a binned contrast map, wherein we counted the number of successful detections (S/N ≥ 5) by both ADI and BDI in 01 × 0.25 mag boxes. Figure 4(b) shows that BDI tends to detect a higher fraction of planets closer to the star than ADI. In terms of planet masses, using the COND mass–luminosity evolutionary tracks (Baraffe et al. 2003) for a 130 Myr old star, a half magnitude improvement translates into ∼1 Jupiter mass improvement in sensitivity for BDI over ADI.

Zoom In Zoom Out Reset image size

5.1. Benefits of BDI

5.2. Limitations

We have shown that binary stars should be considered viable future targets for large direct imaging surveys because of the potential for improved PSF subtraction close to the stars using BDI. This is important because planets are thought to be rare at wide separations. Based on the results of our proof of concept BDI data, we are conducting a new direct imaging search for exoplanets in young visual binaries with MagAO/Clio-2. This survey should be able to detect fainter planets closer to their host stars and will tell us about exoplanet demographics in stellar populations that have never been probed before. BDI on a space-based telescope like JWST/NIRcam would not be limited by isoplanatism effects and would therefore be an even more powerful technique for imaging and discovering exoplanets.

We thank the anonymous referee for helpful comments and suggestions. T.J.R. acknowledges support for Program number HST-HF2-51366.001-A, provided by NASA through a Hubble Fellowship grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. E.E.M. acknowledges support from NSF grant AST-1313029. J.R.M. was supported under contract with the California Institute of Technology (Caltech) funded by NASA through the Sagan Fellowship Program.