Reference star differential imaging of close-in companions and circumstellar disks with the NIRC2 vortex coronagraph at W.M. Keck Observatory

TL;DR

This paper demonstrates that reference star differential imaging (RDI) with the Keck/NIRC2 vortex coronagraph enhances detection sensitivity for close-in companions and circumstellar disks, outperforming ADI in certain scenarios.

Contribution

The study shows that RDI improves detection sensitivity and accuracy in high contrast imaging, especially for close-in companions and circumstellar disks, compared to traditional ADI methods.

Findings

RDI increases detection significance of point sources by up to 5 times compared to ADI.

Pre-selection of reference frames improves detection significance by up to 3 times.

RDI enables accurate mapping of circumstellar disks without self-subtraction artifacts.

Abstract

Reference star differential imaging (RDI) is a powerful strategy for high contrast imaging. Using example observations taken with the vortex coronagraph mode of Keck/NIRC2 in $L^\prime$ band, we demonstrate that RDI provides improved sensitivity to point sources at small angular separations compared to angular differential imaging (ADI). Applying RDI to images of the low-mass stellar companions HIP 79124 C (192 mas separation, $\Delta L^\prime$=4.01) and HIP 78233 B (141 mas separation, $\Delta L^\prime$=4.78), the latter a first imaging detection, increases the significance of their detections by up to a factor of 5 with respect to ADI. We compare methods for reference frames selection and find that pre-selection of frames improves detection significance of point sources by up to a factor of 3. In addition, we use observations of the circumstellar disks around MWC 758 and 2MASS J16042165-2130284 to show that RDI allows for accurate mapping of scattered light distributions without self-subtraction artifacts.

Equations

$$I=\sum _{k=1}^K⟨X,Z^{(k)}⟩Z^{(k)},$$

$${\mathrm{MSE}}^{(k)}=\frac{1}{N_{\mathrm{pix}}}\sum _{i=1}^{N_{\mathrm{pix}}}{\left(X_i^{(k)}-M_i\right)}^2,$$

$${\mathrm{PCC}}^{(k)}=\frac{\mathrm{cov}(X^{(k)},M)}{\mathrm{std}(X^{(k)})\mathrm{std}(M)},$$

$$\mathrm{cov}(P,Q)=\frac{1}{N_{\mathrm{pix}}-1}\sum _{i=1}^{N_{\mathrm{pix}}}(P_i-\bar{P})(Q_i-\bar{Q}),$$

$$\mathrm{std}(P)=\sqrt{\frac{1}{N_{\mathrm{pix}}-1}\sum _{i=1}^{N_{\mathrm{pix}}}(P_i-\bar{P}{)}^2}.$$

$${\mathrm{SSIM}}^{(k)}=\frac{1}{N_{\mathrm{pix}}}\sum _{i=1}^{N_{\mathrm{pix}}}{L}_i^{(k)}{C}_i^{(k)}{S}_i^{(k)},$$

$$L_i=\frac{2\bar{X}\bar{M}+c_1}{{\bar{X}}^2+{\bar{M}}^2+c_1},$$

$$C_i=\frac{2\mathrm{std}(X)\mathrm{std}(M)+c_2}{\mathrm{std}(X{)}^2+\mathrm{std}(M{)}^2+c_2},$$

$$S_i=\frac{\mathrm{cov}(X,M)+c_3}{\mathrm{std}(X)\mathrm{std}(M)+c_3}.$$

Full text

Reference star differential imaging of close-in companions and circumstellar disks with the NIRC2 vortex coronagraph at W.M. Keck Observatory

Garreth Ruane

NSF Astronomy and Astrophysics Postdoctoral Fellow

Department of Astronomy, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

Henry Ngo

NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, British Columbia, Canada

Dimitri Mawet

Department of Astronomy, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA

Olivier Absil

F.R.S.-FNRS Research Associate

Space Sciences, Technologies, and Astrophysics Research (STAR) Institute, Université de Liège, Liège, Belgium

Élodie Choquet

Hubble Fellow

Department of Astronomy, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

Aix Marseille Univ., CNRS, CNES, LAM, Marseille, France

Therese Cook

Department of Astronomy, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA

Carlos Gomez Gonzalez

Elsa Huby

LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité,

5 place Jules Janssen, 92195 Meudon, France

Keith Matthews

Department of Astronomy, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

Tiffany Meshkat

IPAC, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

Maddalena Reggiani

Space Sciences, Technologies, and Astrophysics Research (STAR) Institute, Université de Liège, Liège, Belgium

Eugene Serabyn

Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA

Nicole Wallack

Division of Geological and Planetary Sciences, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

W. Jerry Xuan

Department of Astronomy, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA

Department of Physics and Astronomy, Pomona College, 333 N. College Way, Claremont, CA 91711, USA

Université˝́ Grenoble Alpes, IPAG, F-38000 Grenoble, France

1 Introduction

Adaptive optics (AO) imaging surveys determine crucial occurrence statistics for giant exoplanet populations with masses $\gtrsim$1$~{}M_{\mathrm{Jup}}$ at orbital separations $>$10 au and provide unique insight into their formation and migration histories (Bowler, 2016; Meshkat et al., 2017). Furthermore, resolving low-mass companions from their host stars reduces measurement noise allowing for in-depth spectral characterization, which yields a wealth of information about the physical properties of their atmospheres (Biller & Bonnefoy, 2017). AO observations also allow for the detailed mapping of scattered light from dust in circumstellar disks providing context for planet formation theories.

Current and future infrared, high-contrast imaging programs seek to (1) bridge the gap between the search completeness of direct imaging and radial velocity surveys, (2) characterize the atmospheres of lower mass planets, and (3) understand the interaction between planets and the circumstellar disks from which they form. All of the above science cases benefit from minimizing the angular separations from the host star at which instruments are sensitive enough to detect and characterize the planets and disks of interest.

Reducing the inner working angle (i.e. the separation at which the transmission is 0.5) of ground-based, high-contrast imaging systems to their theoretical limits (i.e. $\sim\lambda/D$, where $\lambda$ is the wavelength and $D$ is the telescope diameter; Guyon et al., 2006) opens a critical search space at solar system scales (1-30 au) for stars within a few hundred parsecs. From a technical point of view, it will also pave the way to studies of temperate exoplanets in reflected light with future large aperture telescopes ($D\sim$30 m), especially in the habitable zones of nearby, late-type stars where the planet-to-star flux ratio is expected to be $\sim$10*-8* for planet radii of $\sim 1R_{\earth}$ (see e.g. Guyon, 2018). Looking further into the future, these developments will also help maximize the number of Earth-like planets orbiting solar-type stars available for study with future exoplanet imaging space missions, such as the Habitable Exoplanet Observatory (HabEx; Gaudi et al., 2018) and Large UV/Optical/IR Surveyor (LUVOIR; The LUVOIR Team, 2018) concepts.

The vortex coronagraph mode (Foo et al., 2005; Mawet et al., 2005) of the NIRC2 instrument at W.M. Keck Observatory (Serabyn et al., 2017; Mawet et al., 2017; Ruane et al., 2017; Huby et al., 2017; Guidi et al., 2018; Reggiani et al., 2018; Xuan et al., 2018) has the optical throughput needed to search for giant planets down to $\lambda/D$, or $\sim$100 mas, for $L^{\prime}$ and $M_{s}$ bands ($\lambda$=3.4-4.8 $\mu$m). However, pairing the Keck adaptive optics system with a small inner working angle coronagraph does not immediately provide sensitivity to point sources with $\Delta$mag$\gtrsim$3 within 0$\farcs$5 of the host star in our experience. Surpassing this level also requires optimized observing strategies and post-processing (Mawet et al., 2012).

1.1 The limitations of ADI

High contrast imaging surveys with Keck/NIRC2 traditionally make use of angular differential imaging (ADI, Marois et al., 2006) in the search for giant planets, such as the four planets discovered orbiting HR 8799 (Marois et al., 2008, 2010; Konopacky et al., 2016). In this strategy, the telescope’s field rotator is set to vertical angle mode, which keeps the beam fixed in azimuthal orientation with respect to the elevation axis and detector while the field of view rotates over the course of an observation. This assists in differentiating between true astrophysical objects, which revolve in the image, and starlight that leaks through the coronagraph, which typically leaves a residue of speckles in the image.

Although ADI is the most successful observational strategy for detecting giant planets at wide separations to date (see recent review by Chauvin, 2018), it has some innate limitations. In order to clearly identify companions, each target must be observed for long enough to allow enough parallactic angle (PA) rotation such that a companion would move by a substantial fraction of the width of the point spread function (PSF). For instance, at the inner working angle of a vortex coronagraph (0$\farcs$08-0$\farcs$1 for NIRC2 $L^{\prime}$ and $M_{s}$ bands), the minimum desired PA rotation is $\sim$30∘, which moves the PSF of a companion by an arclength of roughly half of its full width at half maximum (FWHM). However, due to practical limitations, observations from ongoing surveys with Keck/NIRC2 vortex coronagraph have a median PA rotation of $\sim$11∘ (Xuan et al., 2018).

To illustrate the problems associated with scheduling ADI observations, Figure 1 shows the PA rotation achieved from Maunakea, Hawaii (latitude 19.8∘ N) in a one hour observation versus the target’s declination and the timing of the observation. An “on-time” observation (thick, blue, dashed line) in this case starts 30 min before the target crosses the meridian. The thin black lines shows how the PA rotation degrades when observation is initiated too early or too late by multiples of $\Delta t$=15 min.

Figure 1 has several important implications:

•

Timing constraints: ADI observations must be carried out during a relatively narrow time window to achieve a PA rotation needed to avoid self-subtraction effects (Marois et al., 2006). For example, in order for a planet to move 30∘, or roughly an arclength of 0.5$\times$FWHM at an angular separation of $\lambda/D$, an observation at the most favorable $\delta$ must be initiated within $\pm$68 min of the optimal time. Increasing the path length requirement to 1.0$\times$FWHM (i.e. a PA rotation of 60∘) reduces the allowable timing error to $\pm$36 min. The timing requirement is even more strict at other $\delta$ values.

•

Limited sky coverage: Although targets with declination $>$-40∘ are observable from Maunakea, achieving a PA rotation of $>$30∘ during a one hour observation is only possible on targets between $\delta$=-8∘ and $\delta$=48∘ (assuming the observations are timed perfectly). For PA rotation $>$60∘, the $\delta$ range becomes 8∘ to 32∘, ruling out more than 90% of the sky.

•

Inefficient surveys of star forming regions: Many of the most attractive targets for high-contrast imaging in the infrared are in nearby star forming regions (distances 120-150 pc), such as the Taurus and $\rho$ Ophiuchi molecular clouds (see e.g. Bowler, 2016). However, these regions only extend over a few square degrees, allowing for one or two ADI sequences on these targets per night. Therefore, an ADI survey of these regions would need to be carried out over many nights.

•

Limited “effective” inner working angle: Unless the star falls within the declination ranges described above and the observations are well-timed, self-subtraction effects at small angular separations (Marois et al., 2006) ultimately limit the sensitivity near to the star and inner working angle of the coronagraph.

•

Erases rotationally symmetric circumstellar disks: Imaging scattered light from close to pole-on protoplanetary, transitional, and debris disks with ADI may not be possible because rotationally symmetric features self-subtract in post-processing and complex dust disk distribution become confounded by artifacts (Milli et al., 2012).

1.2 The limitations of SDI

An alternate speckle subtraction method known as spectral differential imaging (SDI; Sparks & Ford, 2002) is applied on instruments with integral field spectrographs, such as Palomar/P1640 (Hinkley et al., 2011), Gemini/GPI (Macintosh et al., 2014), VLT/SPHERE (Vigan et al., 2016), and Subaru/SCExAO/CHARIS (Jovanovic et al., 2015; Groff et al., 2017). SDI decouples speckles from astrophysical objects based on their wavelength dependence rather than azimuthal rotation and therefore alleviates some of the aforementioned limitations of ADI. Between the shortest and longest wavelengths, speckles move radially by $\alpha\times\Delta\lambda/\lambda$, where $\alpha$ is the angular separation in units of $\lambda/D$. For instance, a radial change in the position of $>\lambda/D$ is only seen at angular separations $>$5 $\lambda/D$ using a typical filter with $\Delta\lambda/\lambda\approx 0.2$. Therefore, SDI also suffers from self-subtraction effects and is not ideal for imaging close-in point sources or pole-on circumstellar disks. Additionally, Keck/NIRC2 does not have an integral field spectrograph for this purpose, making contemporaneous SDI impossible.

1.3 Reference star differential imaging

Reference star differential imaging (RDI; Lafrenière et al., 2009; Soummer et al., 2011; Gerard & Marois, 2016) is an alternative approach that uses images of other stars to build a model of the stellar PSF. RDI is a commonly used observational strategy for Hubble Space Telescope imaging observations of debris disks (Golimowski et al., 2006; Schneider et al., 2009, 2014; Choquet et al., 2016) and vortex coronagraph observations with the Hale telescope at Palomar Observatory (Mawet et al., 2010; Serabyn et al., 2010; Mawet et al., 2011) where ADI is not possible due to its equatorial mount. In addition, RDI will very likely be applied for high-contrast imaging with future space telescopes, including the James Webb Space Telescope (JWST; Green et al., 2005) and the Wide Field Infrared Survey Telescope (WFIRST; Spergel et al., 2015; Bailey et al., 2018) as well as the HabEx (Gaudi et al., 2018) and LUVOIR (The LUVOIR Team, 2018) mission concepts. The latter two may push to very small angular separations using vortex or other small inner working angle coronagraphs (Ruane et al., 2018; Pueyo et al., 2017).

Here, we demonstrate how RDI with Keck/NIRC2 mitigates practical issues associated with ADI (and potentially SDI on other instruments). We demonstrate the benefits of the RDI observing strategy for ground-based imaging studies, especially the detection of point sources at small angular separations and mapping of scattered light from circumstellar disks. In the following sections, we detail the RDI method and strategy (Sec. 2) using three example observing nights that illustrate the benefits of RDI for point source detection (Sec. 3) and for imaging of circumstellar disks (Sec. 4) with Keck/NIRC2. Section 5 discusses our findings and the limitations of RDI, including error sources and potential artifacts. Section 6 summarizes our conclusions.

2 RDI method and strategy

RDI requires the observer to decide which reference frames and algorithms to use in order to model the stellar PSF. Reference frames may be specific reference stars purposely taken on the same night as the primary observation (e.g. Ruane et al., 2017), frames from other stars that were coincidentally observed on the same night in the same observing mode (e.g. Xuan et al., 2018), or from an archive of frames across observing programs (e.g. Choquet et al., 2016). Example PSF reconstruction algorithms include subtracting a scaled version of the reference image (Schneider et al., 2009), principal components analysis (PCA; Soummer et al., 2012), and non-negative matrix factorization (NMF; Ren et al., 2018).

Xuan et al. (2018) showed that using all frames from a given night (excluding the target of interest) as reference frames and applying PCA improves the sensitivity to point sources at small angular separations with respect to ADI for many of the targets in our database of Keck/NIRC2 vortex coronagraph observations. In the following, we will discuss improvements to this strategy by combining this approach with frame pre-selection.

For each observation, we apply basic pre-processing steps (see Appendix A for details) and PCA to estimate and subtract the starlight from the images using the Vortex Image Processing (VIP) software package (Gomez Gonzalez et al., 2017). Following Soummer et al. (2012), a single frame may be written as $X=I+A$, where $I$ is the stellar PSF after the coronagraph (a speckle field) and $A$ is the planet PSF. The stellar PSF is

$I=\sum_{k=1}^{K}\left<X,Z^{(k)}\right>Z^{(k)},$

(1)

where $\{Z^{(k)}\}_{k=1,...,N}$ are a basis set of images derived from a series of $N$ reference frames $\{R^{(k)}\}_{k=1,...,N}$, $\left<X,Z^{(k)}\right>$ is the projection of the frame onto the $k$th basis image, and $K$ represents the number of basis vectors used for the PSF model with $K\leq N$. Here, the basis modes are computed using PCA. In the case of point sources, we empirically determine the number of basis vectors, $K$, that provides the best sensitivity to fake companions injected in post-processing. Whereas in ADI-processed images the value of $K$ can drastically affect the sensitivity to point sources (Pueyo, 2016) and the occurrence of disk artifacts (Milli et al., 2012), RDI results tend to be relatively insensitive to the choice of $K$ (Soummer et al., 2012). Rather, the choice of which reference frames to use to model the stellar PSF has a dominant effect on the quality of the final image.

Table 1: Keck/NIRC2 vortex observations in $L^{\prime}$ band on UT 2016 Apr 13

Note. — aWise $W1$ mag (3.4 $\mu$m; Wright et al., 2010).

Table 2: Best fit photometry and astrometry for the detected companions.

Note. — aBailer-Jones et al. (2018). bLower bound on mass; assumes the BT-Settl model (Allard et al., 2012) and an age of 10 Myr.

We explored a number of different methods for reference frame pre-selection. In each case, we compared the potential reference frames, $X_{i}$, with the temporal median of the science frames, $M$, over a speckle noise dominated region of the image (typically within 5-10 $\lambda/D$). We assigned a score to each reference frame using three metrics: the mean square error (MSE), the Pearson correlation coefficient (PCC), and the structural similarity index metric (SSIM) as defined in Wang et al. (2004). The MSE of the $k$th reference frame is given by

$\mathrm{MSE}^{(k)}=\frac{1}{N_{\mathrm{pix}}}\sum_{i=1}^{N_{\mathrm{pix}}}\left(X_{i}^{(k)}-M_{i}\right)^{2},$

(2)

where $N_{\mathrm{pix}}$ is the number of pixels in the comparison region. The PCC is

$\mathrm{PCC}^{(k)}=\frac{\mathrm{cov}(X^{(k)},M)}{\mathrm{std}(X^{(k)})\mathrm{std}(M)},$

(3)

where $\mathrm{cov}(.)$ and $\mathrm{std}(.)$ represent the covariance and standard deviation. Specifically, the covariance is given by

$\mathrm{cov}(P,Q)=\frac{1}{N_{\mathrm{pix}}-1}\sum_{i=1}^{N_{\mathrm{pix}}}(P_{i}-\bar{P})(Q_{i}-\bar{Q}),$

(4)

where $P$ and $Q$ are the images being compared, with pixel-wise means $\bar{P}$ and $\bar{Q}$. The standard deviation is

$\mathrm{std}(P)=\sqrt{\frac{1}{N_{\mathrm{pix}}-1}\sum_{i=1}^{N_{\mathrm{pix}}}(P_{i}-\bar{P})^{2}}.$

(5)

The final metric we consider is the mean SSIM:

$\mathrm{SSIM}^{(k)}=\frac{1}{N_{\mathrm{pix}}}\sum_{i=1}^{N_{\mathrm{pix}}}L_{i}^{(k)}C_{i}^{(k)}S_{i}^{(k)},$

(6)

where $L$, $C$, and $S$ are the luminance, contrast, and structural terms, which are computed over a FWHM$\times$FWHM window centered on pixel $i$. The luminance term is

$L_{i}=\frac{2\;\bar{X}\bar{M}+c_{1}}{\bar{X}^{2}+\bar{M}^{2}+c_{1}},$

(7)

the contrast term is

$C_{i}=\frac{2\;\mathrm{std}(X)\mathrm{std}(M)+c_{2}}{\mathrm{std}(X)^{2}+\mathrm{std}(M)^{2}+c_{2}},$

(8)

and the structural term is

$S_{i}=\frac{\mathrm{cov}(X,M)+c_{3}}{\mathrm{std}(X)\mathrm{std}(M)+c_{3}}.$

(9)

$c_{1}$, $c_{2}$, and $c_{3}$ are small constants chosen to prevent instability when the denominator is otherwise close to zero. Qualitatively, $L_{i}$ is the relative change in luminance and responds similarly to MSE to differences in the pixel counts. $S_{i}$ is almost identical to the PCC. Thus, SSIM can be thought of as a mixture between the MSE and PCC metrics.

In the following section, we investigate the performance gains achieved by modeling the stellar PSF using the principal components (PCs) of a set of reference frames selected based on their MSE, PCC, and SSIM values, as compared to the median science frame.

3 Point source detection

To illustrate the utility of and trade-offs associated with using RDI for point source detection, we revisit Keck/NIRC2 observations taken during commissioning of the vortex coronagraph mode on UT 2016 Apr 13. Table 1 gives the full list of observations, consisting of 11 targets observed over the course of a full observing night for various scientific purposes. The targets span a large range of elevations, spectral types (Spec. type), and magnitudes, which we list in $V$ and Wise $W1$ bands (3.4 $\mu$m; Wright et al., 2010). The observations also vary in the number of frames ($N_{\mathrm{frames}}$) and total integration times ($t_{\mathrm{int}}$).

Hinkley et al. (2015) previously identified low-mass stellar companions orbiting at $\sim$20 au from two of the observed targets, HIP 79124 and HIP 78233, using aperture masking interferometry. HIP 79124 and HIP 78233 are classified as A0V and F0V stars, which reside in the Upper Scorpious subgroup of the Scorpius-Centaurus (USco) association at a distance of $\sim$120-150 pc (de Zeeuw et al., 1999) and, thus, have an estimated age of 5-20 Myr (Pecaut et al., 2012; Song et al., 2012). HIP 79124 is a triple system with an additional companion, HIP 79124 B, at $\sim$1*′′* from HIP 79124 A (Lafrenière et al., 2014). RDI allows us to directly image their close-in companions, HIP 79124 C and HIP 78233 B, even though their angular separation from the host star is $<$0$\farcs$2, or $2~{}\lambda/D$ in $L^{\prime}$ band. We use these previously reported astrometric measurements along with data presented in this study (Table 2) to check for common proper motion and confirm that HIP 78233 B, HIP 79124 B, and HIP 79124 C are all gravitationally bound to their host star.

Figure 2 shows the MSE, PCC, and SSIM for each possible reference frame from UT 2016 Apr 13, compared to the median frame in the HIP 79124 and HIP 78233 observations. As predicted, the MSE and PCC respond to different attributes, whereas the SSIM shows similar features to both MSE and PCC. In the case of HIP 79124 and HIP 78233, the brighter stars (HIP 66704, DX Leo, HIP 41152) had considerably worse MSE and SSIM. Also, HIP 79124 and HIP 78233 were mutually poor reference stars for one another according to the PCC and the SSIM despite being observed consecutively. This is likely due to variable observing conditions at that point in the night.

Figures 3 and 4 show the ADI and RDI images of HIP 79124 C and HIP 78233 B, respectively. HIP 79124 C was previously imaged by Serabyn et al. (2017) using RDI. However, we have recovered HIP 78233 B for the first time. We show both examples here to demonstrate that, using PCA, the companions are detected at higher significance with RDI compared to our ADI reduction.

In order to determine which reference frames to use, we ordered the reference frames according to the MSE, PCC, and SSIM and varied the number of frames from which we derive the PCs (see Fig. 5). In each case, we compute a forward model of the companion’s PSF (see e.g. Lagrange et al., 2009; Soummer et al., 2012). We subtracted a copy of the off-axis PSF at the location of the planet in each pre-processed science frame, varied its location and brightness, and repeated the PCA reduction until the values were minimized in the final images in a 2$\times$FWHM radius about the companion’s position using a downhill simplex algorithm. The forward model is defined as the difference between the original PCA reduction without injected or subtracted companions and the best-fit residuals.

The signal-to-noise ratio, $S/N$, is defined as the ratio between estimated flux, $F$, and the error in the flux fit, $\sigma_{F}$, as calculated using the method outlined in Appendix B. We find that the MSE and SSIM metrics lead to higher $S/N$ for these objects compared to PCC, with SSIM offering a slight improvement. The optimal number of frames according to SSIM is 237 (out of 336) and 202 (out of 323) for HIP 79124 C and HIP 78233 B, respectively. We also optimized $K$ in this fashion for each number of reference frames considered, but found that $K$ has much less of an influence on the $S/N$. At the optimal number of reference frames, we used 11 and 15 PCs for HIP 79124 C and HIP 78233 B, respectively. Figures 3 and 4 show the result using the optimal combination of parameters. The $S/N$ has a steep drop when using more than $\sim$230 frames ordered by MSE and SSIM, which occurs when the reduction includes frames from the brightest stars (HIP 66704, DX Leo, HIP 41152). This is expected since we found previously in Xuan et al. (2018) that the contrast in Keck/NIRC2 images is strongly correlated with stellar magnitude.

Table 2 reports the astrometry and photometry values derived from the highest $S/N$ case. There are a few important differences between the values in Fig. 5 and the final fits. First, we used the flux and position we retrieved using all reference frames to compute the fitting error in Fig. 5. In the final fits, we used the actual flux and position retrieved at the optimal number of reference frames and PCs for the injected, fake companion. Second, we applied a correction to the retrieved flux values to account for the throughput of the coronagraph, which was 81% and 68% for HIP 79124 C and HIP 78233 B, respectively, according to our optical model. In the end, we arrived at lower $S/N$ values for each companion than the estimates in Fig. 5. The final $S/N$ was better with RDI than ADI for both stars. Specifically, for HIP 79124 C, the $S/N$ in RDI was 36 whereas the best we could achieve in ADI was 7.9 with $K=2$ and a PA rotation of 17∘. Similarly, for HIP 78233 B, the $S/N$ in RDI was 8.6, but the best we achieved in ADI was 5.1 with $K=6$ despite having a PA rotation of 32∘.

Comparing to Hinkley et al. (2015), the change in astrometry since 2010 is consistent with orbital motion, leading to the first confirmation that HIP 78233 B is bound. We updated the distances to HIP 79124 and HIP 78233 based on Gaia DR2 (Bailer-Jones et al., 2018), which gave 132$\pm$1 pc and 256$\pm$4 pc for the respective stars, rather than the distances of 123 pc and 145 pc used in Hinkley et al. (2015). Our measured flux ratio for HIP 79124 C, $\Delta L^{\prime}=4.01\pm 0.03$, is inconsistent with both Hinkley et al. (2011) and Serabyn et al. (2017) by a small margin, who report $\Delta L^{\prime}=4.3\pm 0.1$ and $\Delta L^{\prime}=4.2\pm 0.1$, respectively. The increased flux estimate and updated stellar distance predicts a higher mass for both companions. Assuming the BT-Settl evolutionary model (Allard et al., 2012) and an age of 10 Myr, we estimate minimum masses of 199$\pm$4 $M_{\mathrm{Jup}}$ and 198$\pm$19 $M_{\mathrm{Jup}}$, which are well above the hydrogen-burning limit.

For each point source companion, we use previously reported astrometry along with new data presented in this study to constrain the orbits of each companion around their host stars. Although two data points do not provide a well-constrained orbit, the long baseline provides some constraint on the actual semimajor axis and eccentricity. We use the “Orbits of the Impatient” (OFTI; Blunt et al. 2017) algorithm as implemented by the orbitize+ software package111https://orbitize.readthedocs.io/. OFTI is a Bayesian rejection-sampling method that excels at finding constraints on long period orbits with sparse data points. We determine the median semimajor axis and 68% confidence intervals to be $37^{+31}_{-12}$ au for HIP 78233 B, $157^{+82}_{-37}$ au for HIP 79124 B and $26^{+20}_{-8}$ au for HIP 79124 C. We also found that the data favors lower eccentricity orbits over high eccentricity orbits. Details of these orbit fits and the asymmetrical posteriors are discussed in Appendix C. The posteriors are also available online as supplementary data for use in other studies.

Table 3: Keck/NIRC2 vortex observations in $L^{\prime}$ band on UT 2016 Oct 24.

Name Alt. Name RA Dec Type $N_{\mathrm{frames}}$ $t_{\mathrm{int}}$ (min) PA rot. $V$ $W1$

TYC 1845-2048-1 HD 284115 05 02 41.3 +24 07 19.5 K2 25 16.7 23.1∘ 8.13 3.98$\pm$0.10

MWC 758 HD 36112 05 30 27.5 +25 19 57.1 A8V 80 53.3 173∘ 8.27 4.60$\pm$0.08

TYC 1867-221-1 HD 249769 05 58 47.2 +25 15 28.8 K7 50 33.3 171∘ 8.83 4.14$\pm$0.09

In the two examples above, we have demonstrated that RDI offers significantly improved $S/N$ for companions at angular separations of $\lesssim 2~{}\lambda/D$ and that frame pre-selection is critical for maximizing the performance achieved with RDI. However, previous versions of the data reduction pipeline used by our team for Keck/NIRC2 vortex coronagraph observations (Xuan et al., 2018) used all available reference frames from a given night by default. Although there is often an $S/N$ improvement at small separations using all frames over ADI, we are working towards using SSIM or similar metrics to perform the pre-selection of frames in a robust way and reprocessing archival NIRC2 data to reveal new point sources at close separations.

4 Circumstellar disk imaging

In addition to improving the sensitivity to point sources at small angular separations, RDI is also advantageous for disk imaging. In the following two examples, we use Keck/NIRC2 observations of circumstellar disks. In each case, reference stars were specifically chosen to be close in elevation when observed and to have similar brightness in $R$ and $L^{\prime}$ bands, where the wavefront sensor and science camera are sensitive.

4.1 MWC 758

MWC 758 (HD 36112) is a young (3.5$\pm$2.0 Myr; Meeus et al., 2012) stellar object in the Taurus star-forming region (distance of 160 pc; Bailer-Jones et al., 2018) surrounded by a prominent protoplanetary disk with at least two massive spiral arms seen in scattered light observations (Grady et al., 2013; Benisty et al., 2015; Reggiani et al., 2018; Ren et al., 2018). ADI images with Keck/NIRC2 from Reggiani et al. (2018) showed evidence of a potential third spiral arm.

Here, we reprocessed one of the epochs presented in Reggiani et al. (2018) using reference stars observed directly before and after the previously published ADI sequence. Table 3 gives the list of observations we used in the reduction. For ADI, we used all of the science frames and $K=6$. For RDI, we used similar methods as the previous section keeping the best 70% (50 out of 75) of reference frames according to the SSIM metric and used the projection of 25 out of 50 PCs in the PSF model. However, since the reference stars were a good match to the science targets and conditions were stable during these observations, the variation in MSE, PCC, and SSIM was similar to that of one of the well-correlated observations in Fig. 2. Therefore, the RDI image of the disk is visibly unchanged, barring some background variation, whether we use 10% more or less frames or change the number of PCs by $\sim$20%, while the ADI image is very sensitive to number of PCs.

Figure 6 shows the ADI and RDI images of MWC 758. In addition to the famous North and South spirals (labeled ‘1’ and ‘2’ in Fig. 6b), we recovered the proposed third spiral arm feature (labeled ‘3’) with both ADI and RDI for the first time. The ADI image is dominated by self-subtraction effects, while the RDI image is likely more representative of the true morphology of scattered light features. Although we have confirmed that the proposed third spiral arm was not an ADI artifact, it is unclear from the RDI image whether the feature is truly a third spiral arm or a continuation of the north arm wrapping almost 360∘ about the star. Nevertheless, this observation demonstrates that RDI is beneficial for mapping scattered light distributions in circumstellar disks without the sharpening of azimuthal features seen in the ADI images.

4.2 J16042165-2130284

We present an observation of the young stellar object 2MASS J16042165-2130284 (from here on ‘J1604’). J1604 is a K2 star at a distance of 150 pc (Bailer-Jones et al., 2018) and a likely member of USco (Carpenter et al., 2014) suggesting an age of 5-20 Myr (Pecaut et al., 2012; Song et al., 2012). J1604 is surrounded by a nearly pole-on (7∘ inclination; Mathews et al., 2012) pre-transitional disk (Espaillat et al., 2007), the most massive known ($\sim 0.1~{}M_{\mathrm{Jup}}$) in USco (Mathews et al., 2012). Several previous works identified a ring-like disk in scattered light (Mayama et al., 2012; Pinilla et al., 2015; Canovas et al., 2017; Pinilla et al., 2018). This disk has received considerable attention owing to its variability (Pinilla et al., 2015, 2018) and evidence of planet-induced dust filtration (Rice et al., 2006; Canovas et al., 2017).

The observing list from UT 2017 May 10 (see Table 4) consists of seven observations, including J1604 and a reference star, 2MASS J19121875-2137074 (from here on ‘J1912’), specifically selected prior to the observing night because its similar brightness and elevation. The remaining five targets were observed as part of an ongoing survey of M stars (PI: Mawet). The first three M star targets were significantly worse reference stars according to their SSIM values because they were much fainter ($W1>9$) than J1604 ($W1=7.61$). The other two M stars had SSIM values close to that of J1912.

Figure 7 shows the ADI and RDI images of J1604. Since the disk is nearly rotationally symmetric, the ADI image is consistent with noise for all values of $K$. However, the ring of scattered light is clearly visible in the RDI reduction. For the RDI image in Fig. 7b, we used the best 88 reference frames, according to their SSIM values, out of the 150 available and projected 44 PCs to build the stellar PSF model. Again, the morphology of the disk did not change much when changing the number of reference frames or PCs by $\sim 20\%$ in this case. After removing poorly matching reference frames according to the SSIM, the projection coefficients for the first few PCs dominate.

Polarized intensity maps of J1604 from Pinilla et al. (2018) show azimuthal dips attributed to shadowing by an inner disk. We see similar features labeled ‘1’, ‘2’, and ‘3’ in Fig. 7b. These dips change with time indicating that the inner regions are highly dynamic, which may be evidence of a close-in massive companion. In the southern part of the disk, we estimate the inner and outer edges of the bright scattering feature appear at 440 and 500 mas, respectively. Compared to $J$-band (1.2 $\mu$m) observations taken with VLT/SPHERE by Pinilla et al. (2018) who measured a scattered light peak at 430 mas, the scattering surface appears slightly further from the star in $L^{\prime}$ band (3.8 $\mu$m), which may be evidence of spatial segregation of dust particle sizes and/or lower opacities at longer wavelengths. In future work, we will fit a forward model of the disk using radiative transfer and further investigate the wavelength dependence of the scattered light surface and shadowing effects in the disk (Wallack et al., in prep.), all of which would not be possible using ADI.

Table 4: Keck/NIRC2 vortex observations in $L^{\prime}$ band on UT 2017 May 10.

Name RA Dec Type $N_{\mathrm{frames}}$ $t_{\mathrm{int}}$ (min) PA rot. $V$ $W1$

2MASS J11110358-3134591 11 11 03.6 -31 34 59.1 M 20 15.0 6.99∘ 14.4 9.35$\pm$0.02

2MASS J11431742+1123126 11 43 17.4 +11 23 12.6 M 20 15.0 10.1∘ 12.3 9.12$\pm$0.03

2MASS J13412668-4341522 13 41 26.7 -43 41 52.2 M3.5 20 15.0 6.18∘ 14 9.70$\pm$0.02

2MASS J16042165-2130284 16 04 21.7 -21 30 28.5 K2 73 54.8 36.2∘ 11.9 7.61$\pm$0.03

2MASS J18580415-2953045 18 58 04.2 -29 53 04.5 M0V 20 15.0 7.25∘ 11.8 7.86$\pm$0.02

2MASS J19121875-2137074 19 12 18.8 -21 37 07.4 M 50 37.5 19.8∘ 11.3 7.42$\pm$0.03

2MASS J20013718-3313139 20 01 37.2 -33 13 13.9 M1 20 15.0 6.95∘ 12.3 8.14$\pm$0.02

5 Discussion

The above sections exemplify that (1) combining the small inner working angle vortex coronagraph and RDI is a sensitive technique for discovering faint companions $\sim$100 mas from the star and imaging scattering light from circumstellar disks with Keck/NIRC2 and (2) optimizing the reference frame selection process leads to significant improvement in the quality of the stellar PSF subtraction. In this section, we describe the noise characteristics, artifacts, and limitations unique to this observing mode.

5.1 Treatment of false positives

Given the typical number of reference frames available per night and the limited diversity they provide, subtraction of the stellar speckle field is imperfect, and the residual speckles in RDI images tend to obey a highly non-Gaussian noise distribution (Goodman, 1975; Aime & Soummer, 2004), especially within a few $\lambda/D$ of the star. However, the standard detection methods used in high-contrast imaging, including this work, estimate the standard deviation of the speckle noise, $\sigma_{n}$, as a function of angular separation from the star and set a fixed detection threshold (typically $5\sigma_{n}$), which leads to a substantial increase in the number of false positives compared to what is expected from normally distributed noise. In addition, since the spatial scale of speckles is roughly the FWHM of the off-axis PSF, $\sigma_{n}$ is estimated using only a few independent samples at small angular separations (Mawet et al., 2014). It is therefore especially challenging to differentiate between bright speckles and true companions in a single observation and to set robust upper limits on the brightness of unseen companions (see e.g. Ruffio et al., 2018) within a few $\lambda/D$ of the star. More work is needed to accurately model the residual speckle noise distribution given a set of science and reference frames.

5.2 Typical artifacts

In addition to the non-Gaussian properties of the stellar speckle noise, there are other artifacts that may lead to false positives. For instance, several dust spots on the vortex focal plane mask in Keck/NIRC2 resemble point sources in RDI images because their brightness is not correlated with that of the star and background and therefore they are not well subtracted in post-processing. This requires the observer rule out detections corresponding to the location of the dust. We find that having a large amount of PA rotation helps blur these effects in the de-rotation step.

The RDI approach assumes that stellar speckle noise dominates in both the science and reference frames. Generally speaking, RDI tends not to provide improved sensitivity at small angular separations on targets where spatially resolved circumstellar material dominates the speckle noise. In fact, projecting reference frames of point sources onto the science frames in an attempt to model the stellar PSF modifies the appearance of the circumstellar disk and generates additional speckle noise, which may obscure the true scattered light features and generate point-like false positives. Spatially resolved sources also have a mismatch in the radial size of the main stellar residuals that leaves behind a ring-like residual at the outer edge of the central lobe of the stellar PSF. A similar effect may occur for objects with vastly different colors, though we have not come across a clear example where chromatic effects dominate.

Lastly, there are speckles that appear in NIRC2 images at $\sim 7~{}\lambda/D$ from the star whose brightness depends on the telescope elevation. We attribute these speckles to segment-to-segment phasing errors that depend on the direction of the gravity vector with respect to the primary mirror and the azimuthal position depends on the position of the rotator. For this reason, it may be beneficial to use reference stars at similar declination to the science target and similar instrument settings, including the rotator angle.

5.3 Sources of night-to-night variance

The methods outlined above may be generalized to include reference frames from multiple observing nights. However, although the stellar PSF is relatively stable throughout an observing night, it tends to vary considerably from night-to-night for a number of reasons.

The calibration of static wavefront error within NIRC2 uses a phase retrieval algorithm to reconstruct the wavefront from defocused images of the internal source. The solution changes on a nightly basis. Additionally, wavefront errors due to the primary mirror must be sensed and corrected by the Shack-Hartmann wavefront sensor in the Keck AO system (Wizinowich et al., 2000), but segment piston errors (or the “terrace” modes) are not seen by the wavefront sensor and are therefore not corrected. This issue may be mitigated by the pyramid wavefront sensor under development at Keck Observatory (Bond et al., 2018).

For these reasons, reference frames are generally not well correlated from night-to-night. Nevertheless, there may be correlations with instrument settings (e.g. the rotator angle and wavefront sensor parameters) or observing conditions that can assist in frame pre-selection, which will be investigated in future work.

6 Conclusion

Using the examples of HIP 79124 C and HIP 78233 B, we have demonstrated that RDI offers improved detection of point sources at small angular separation, by up to a factor of 5 with respect to ADI, further unlocking the potential of the vortex coronagraph on Keck/NIRC2. Furthermore, the observations of MWC 758 and 2MASS J16042165-2130284 show that RDI is beneficial for imaging of circumstellar disks because it preserves the morphology, which allows observers to directly interpret the imaged scattered light distribution. We find that frame pre-selection using image comparison metrics, such as the SSIM, significantly improves the performance of RDI, particularly for point source detection, astrometry, and photometry. RDI may help other current and future ground-based instruments achieve better performance at small angular separations, especially when used in conjunction with small inner working angle coronagraphs, and perform more efficient surveys by relaxing PA rotation requirements. RDI may also be the primary strategy for high-contrast imaging with future space missions, including JWST, WFIRST, HabEx, and LUVOIR.

G. Ruane is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1602444. The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration (NASA). The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

Appendix A Observing and image processing details

In this section, we outline the details of our observations and the image processing steps used in our team’s pipeline (Xuan et al., 2018). All observations were carried out using the Keck II telescope, natural guide star adaptive optics, the NIRC2 instrument, and the vortex coronagraph mode in $L^{\prime}$ band (3.7 $\mu m$). The angular resolution was 0.08*′′* and the plate scale was 0.01*′′* per pixel (Service et al., 2016). The field rotator was set to vertical angle mode, such that the telescope pupil tracks the elevation axis, to enable ADI. In addition to the science frames (see Tables 1, 3, and 4), we took several images of the off-axis PSF with a discrete integration time (DIT) of 0.008 seconds and 100 coadds as well as images of the sky background with integration times matching that of the science and off-axis PSF frames. The sky, off-axis PSF, and background frames were taken approximately every 30 min during the observing sequence. The alignment of the star and the center of the vortex focal plane mask was maintained by the QACITS tip-tilt control algorithm (Huby et al., 2015, 2017).

In the pre-processing stage, bad pixels identified in the dark frames and sky flats were replaced by the median of neighboring values in all frames. The sky flat was the median of 10 images of a blank patch of sky with the coronagraph focal plane mask removed (DIT of 0.75 s, 10 coadds each). We then subtracted sky background frames, taken with the coronagraph in place, from each frame individually using a scale factor to account for background variability. The frames were centered based on the position of the optical vortex core in the median of the science frames and each of individual science and reference frames were registered using the peak of the cross correlation with the median science frame. Our pipeline crops the raw 1024$\times$1024 pixel frames to 587$\times$587, which is the largest allowable square frame centered on the vortex focal plane mask.

We applied PCA (Soummer et al., 2012) to estimate and subtract the stellar contribution from each image using the VIP software package (Gomez Gonzalez et al., 2017). In the case of ADI, the PCs of a subset of pixels in the science frames make up the stellar PSF model. Although our standard pipeline uses a set of default annular reduction regions, we generally modify these regions to bound the companion or disk of interest on a case-by-case basis. For RDI, we derive the PCs and the stellar PSF model only from the reference frames. Each frame is de-rotated such that North is in the vertical direction and the ADI/RDI image is taken to be the temporal median of the de-rotated frames. We repeat these steps for all possible number of PCs. Later, an optimal number of PCs may be defined by maximizing the detection significance of a real or injected companion, depending on the scientific goal. For instance, an observer may wish to use a different number of PCs for point source detection and disk imaging.

Appendix B Error bars for astrometry and photometry of companions

To estimate the errors in the position and flux estimates, with the best-fit companion subtracted from the data, we re-injected the PSF into the pre-processed science frames at the same separation and brightness, but varied the azimuthal angle and retrieved the photometry and astrometry using the same method for each injected companion tracing a full circle about the star. The step size in azimuthal angle was $360^{\circ}/(2\pi r)$, where $r$ is the radial position in units of the FWHM, to ensure that each fit was performed at the location of an quasi-independent speckle sample. We took the standard deviation of the measured flux and position of the injected companions as the error on each parameter. The signal-to-noise ratio, $S/N$, is defined as the ratio between estimated flux, $F$, and the uncertainty in the flux, $\sigma_{F}$.

Appendix C Orbit fitting of point source companions with orbitize and OFTI

In this section, we provide more details on the orbit fit and posteriors on the point source companions from the systems discussed in Section 3. Table 5 summarizes the astrometry used in our fit. We use stellar distances from Bailer-Jones et al. (2018) of 132$\pm$1 pc and 256$\pm$4 pc for HIP 78233 and HIP 79124, respectively. We use total system mass estimates from Hinkley et al. (2015), corresponding to $1.7\pm 0.1$ $M_{\mathrm{Sun}}$ for HIP 78233 and $2.48\pm 0.45$ $M_{\mathrm{Sun}}$ for HIP 79124. We choose uninformed priors for all parameters. Namely, for the semimajor axis, we choose a prior probability distribution that is uniformly linear in log space (Jeffreys Prior); for inclination, we choose a prior probability distribution corresponding to $p(x)\propto sin(x)$; and for all other parameters, we choose a linearly uniform prior probability distribution. We fit the companions as point masses and for the HIP 79124 triple system, we treat each companion independently as their large mutual separation would make interactions between the B and C companions very small. The inclination convention is that $0^{\circ}\leq i<90^{\circ}$ orbits are counter-clockwise and $90^{\circ}\leq i<180^{\circ}$ are clockwise.

Table 5: Input astrometry for orbit fitting via orbitize

UT date Separation (mas) Position Angle (∘) Reference

HIP 78233 B astrometry

2010 Apr 26 $133\pm 3$ $20\pm 1$ Hinkley et al. (2015)

2016 Apr 13 $141\pm 64$ $6\pm 1$ this work

HIP 79124 B astrometry

2008 May 25 $990\pm 1$ $98.11\pm 0.05$ Lafrenière et al. (2014)

2016 Apr 13 $971.5\pm 1$ $100.61\pm 0.03$ this work

HIP 79124 C astrometry

2010 Apr 5 $177\pm 3$ $242\pm 1$ Hinkley et al. (2015)

2016 Apr 13 $192\pm 2$ $246.5\pm 0.8$ this work

Using orbitize’s OFTI algorithm, we sample 100,000 orbits for each of the three point source companions and compute posterior probabilities. Figure 8 shows the marginalized posterior probabilities for orbital semimajor axis, eccentricity, and inclination. These figures also include the prior probability to demonstrate how well the data constrains the orbits. The marginalized posterior on semimajor axis is the best constraint on the companion’s orbit size when accounting for all available astrometry. For HIP 78233 B and HIP 79124 C, while the eccentricity and inclination posteriors are not well constrained and are not very different from the prior. The eccentricity posterior distribution does show that lower values are more likely than higher values and the inclination posterior distribution disfavors a completely edge-on system (inclination of 90∘). For HIP 79124 B, the majority of the eccentricity distribution is at values below 0.5 and the inclination distribution strongly favors an inclination less than 90∘, indicating that the object is orbiting counter-clockwise from our viewing angle. These posteriors are available as part of this article’s online data for use in further calculations as a table of orbit samples.