A Chromaticity Analysis and PSF Subtraction Techniques for SCExAO/CHARIS data

TL;DR

This paper analyzes the performance of the SCExAO/CHARIS system, highlighting chromaticity issues, and introduces a PSF subtraction pipeline optimized for broadband mode, with implications for future exoplanet imaging.

Contribution

The paper develops a new PSF subtraction pipeline with throughput correction for CHARIS broadband data and compares its effectiveness to narrowband methods.

Findings

Chromaticity in SCExAO/CHARIS is worse than temporal stability.

The new pipeline achieves comparable contrast with broadband and narrowband spectral differential imaging.

Future broadband observations can benefit from the developed data processing framework.

Abstract

We present an analysis of instrument performance using new observations taken with the Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS) instrument and the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) system. In a correlation analysis of our datasets (which use the broadband mode covering J through K band in a single spectrum), we find that chromaticity in the SCExAO/CHARIS system is generally worse than temporal stability. We also develop a point spread function (PSF) subtraction pipeline optimized for the CHARIS broadband mode, including a forward modelling-based exoplanet algorithmic throughput correction scheme. We then present contrast curves using this newly developed pipeline. An analogous subtraction of the same datasets using only the H band slices yields the same final contrasts as the full JHK sequences; this result is consistent with our chromaticity analysis, illustrating that PSF subtraction using spectral differential imaging (SDI) in this broadband mode is generally not more effective than SDI in the individual J, H, or K bands. In the future, the data processing framework and analysis developed in this paper will be important to consider for additional SCExAO/CHARIS broadband observations and other ExAO instruments which plan to implement a similar integral field spectrograph broadband mode.

Tables

Table 1. Table 1: Table of our target observing parameters.

UT date	name	$m_V$	$t_{\text{exp}}$	$n_{\text{frames}}$	${\Delta }_{\text{PA}}$
(d/m/yy)	(HD)		(s)		(deg)
9/4/17	70313	5.51	15	107	25
9/4/17	98673	6.43	31	91	22
9/4/17	109085	4.3	13	93	18
9/4/17	125162	4.18	20	93	21
9/4/17	141378	5.53	31	70	21
7/1/18	23267	6.88	60	50	32
8/1/18	“ ”	“ ”	“ ”	40	24
8/1/18	38206	5.73	31	118	27

Full text

planets and satellites: general, techniques: image processing, techniques: high angular resolution

††journal: AJ

1 Introduction

Ground-based exoplanet imaging is currently uniquely poised to understand how giant planets form and evolve over time. Distinct from other exoplanet detection methods, exoplanet imaging is currently sensitive to wide separations (greater than about 10 AU), young ages (less than about 1 Gyr), and high masses (greater than about one Jupiter mass; e.g., Nielsen et al. 2019). However, the push to closer separations, older ages, and lower masses is constantly improving with new instrumentation and processing algorithms.

Post-processing of high contrast imaging data is currently a crucial component of reaching the necessary contrasts in order to detect new exoplanets by subtracting the correlated speckle noise around their host star. There are three main point spread function (PSF) subtraction observing strategies used to improve the exoplanet signal to noise ratio (S/N): spectral differential imaging (SDI; Rosenthal et al., 1996; Racine et al., 1999; Marois et al., 2000; Sparks & Ford, 2002), angular differential imaging (ADI; Marois et al., 2006a), and the commonly used reference differential imaging RDI algorithm. Although these strategies are often combined with more advanced least-squares subtraction algorithms (e.g., Lafrenière et al., 2007; Soummer et al., 2012), all three methods remain fundamentally limited by PSF chromaticity (SDI and RDI) and/or temporal stability (ADI and RDI). Thus, understanding and mitigating the level of chromatic and temporal decorrelation in high contrast imaging datasets is a crucial step towards improving contrast to detect fainter and/or lower mass exoplanets.

Although significant efforts have been made to understand the relative impact of temporal stability vs. chromaticity from a numerical perspective (e.g. Marois et al., 2008; Vigan et al., 2010), a similar comprehensive on-sky analysis with second generation high contrast imaging instruments that are equipped with an integral field spectrograph (IFS) has not yet been completed. Hinkley et al. (2007), and Milli et al. (2016) have carried out temporal analyses on the quasi-static speckle lifetime, but with no corresponding chromaticity analysis. Following from the work of Marois et al. (2010), temporal correlation analyses have informed some new ADI+SDI PSF subtraction pipelines to select reference images ranked by correlation (e.g. Currie et al., 2012; Wang et al., 2015). Marois et al. (2014) presented an analyses of time- vs. wavelength-dependent correlation in individual target images. In this paper we will extend this analysis to multiple target images; later we will show that such correlation can vary significantly depending on observing conditions and on location in the coronagraphic image.

The main focus of this paper is thus to present a detailed analysis of chromaticity and temporal stability using our observations from the Subaru adaptive optics system (AO188; Minowa et al. 2010), Subaru Coronagraphic Extreme Adaptive Optics (SCExAO; Jovanovic et al. 2015b), and Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS; Groff et al. 2016, 2017). We also present new modifications to the ADI+SDI-based template locally optimized combination of images (TLOCI; Marois et al., 2014) PSF subtraction pipeline for these datasets. Our SCExAO/CHARIS datasets use a Lyot coronagraph in the broadband JHK mode, covering the full 1.15 - 2.39 $\mu$m range at R$\sim$18 (Groff et al., 2016); this new low resolution broadband mode, currently unique to CHARIS amongst current high contrast imaging instruments, is thought to provide a potential gain in achievable contrast through SDI PSF subtraction for an achromatic system (Marois et al., 2008, 2014), and so our performance analysis and PSF subtractions methods will focus on impacts of using this broadband mode.

The structure of this paper is as follows: in §2 we present details of our target observations, in §3 we present the architecture of our PSF subtraction pipeline (§3.2) and then use the same reference image selection and correlation calculation procedures to carry out a detailed chromaticity analysis (§3.3), in §4 we present final contrast curves and discuss instrument performance implications for the future, and we conclude in §5. Additional correlation analysis is presented in appendix A.

2 Observations

We observed a handful of young, nearby two temperature-component debris disks, the extrasolar analogues of the Asteroid and Kuiper belts in our Solar System (Wyatt, 2008). We identified targets in Kennedy & Wyatt (2014) that had not yet been observed by the Gemini Planet Imager (GPI; Macintosh et al. 2014) or Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE; Beuzit et al. 2008) and/or are inaccessible to either for being too far north. We also observed two additional targets in 2018A through ancillary time (Hodapp, 2018). Our observed targets are listed in Table 1, including dates of observations, exposure times, and field of view (FOV) rotation.

Observations were made in the standard pupil tracking mode for ADI. The coronagraphic setup in this broadband mode utilized the standard Lyot coronagraph, applying the 217 mas diameter focal plane mask (FPM), for which the inner working angle (IWA; separation at which an off-axis companion would have 50% optical throughput relative to the center of the star) is 113 mas; at the time, no other coronagraph options offered sufficient performance over the full broadband mode of CHARIS. The deformable mirror waffle spots (Marois et al., 2006b; Sivaramakrishnan & Oppenheimer, 2006; Jovanovic et al., 2015a), or “satellite spots,” are applied with a non-standard amplitude for the April 2017 data, and so instead an on-sky non-coronagraphic data sequence (using an ND filter but with the same Lyot stop) is taken to measure the star-to-spot ratio directly (see §3.1). The January 2018 data uses the standard 50 nm amplitude spots and corresponding pre-computed star-to-spot ratio as in Currie et al. (2018a).

3 Data Analysis

3.1 Data Cube Extraction and Setup

We used the CHARIS data reduction pipeline (DRP; Brandt et al., 2017) to generate data cubes from recorded IFS detector images, which uses a least-squares method to extract the flux in each slice of the data cube (Draper et al., 2014). A background subtraction, flat field correction, and bad pixel mask are also applied in the extraction process as described in Brandt et al. (2017). All images are then prepared for subsequent PSF subtraction (§3.2) and correlation analysis (§3.3). The satellite spots are used for image registration, magnification with wavelength to align speckles, flux-normalization to remove the wavelength dependence of the broadband transmission filter and stellar spectrum, and contrast calibration. We chose to magnify the images with wavelength centered around slice 16 of 22 ($\lambda_{0}=2.0\mu$m). Slice 16 was chosen instead of slice 11 ($\lambda=1.7\mu$m; in the middle of the bandpass) because of wavelength-dependent sampling issues. The CHARIS PSF in J and H band is slightly under-sampled, and we found that cubic spline interpolation, used to rotate and magnify the images for our ADI+SDI PSF subtraction algorithm (§3.2), provided a more accurate centroiding precision from numerical interpolation in K band. The four satellite spots all lie at 15.91 $\lambda_{0}/D$ separation, and so measured satellite spot positions in the registered, magnified images are used to calculate a plate scale for contrast curve measurements (§4). We were not able to use any on-sky astrometric data to validate these derived plate scale values. Images are then high-pass-filtered using a $4.1\times 4.1\;\lambda_{0}/D$ median boxcar filter to attenuate the smooth AO halo while retaining a close-to-unity algorithmic throughput of $\lambda/D$-scale static speckles and any possible planet(s). The $\lambda^{2}$ dependence of the satellite spot intensity is also removed during the flux normalization procedure. The satellite spots are omitted from reference image selection in §3.2 and correlation calculations for analysis in §3.3. Finally, all registered, magnified, flux-normalized, high-pass-filtered images are converted into “raw” contrast units (i.e., the pixel value-to-star ratio throughout the whole image) using measured and pre-calibrated star-to-satellite spot ratios, whose data acquisition is described at the end of §2. The same registration, magnification, flux normalization, and high-pass filtering procedure is applied to the custom April 2017 star-to-spot calibration data before computing the ratios.

3.2 PSF Subtraction

Our ADI+SDI PSF subtraction algorithms build on the least-squares-based locally optimized combination of images (LOCI; Lafrenière et al., 2007) algorithm, further developed by Marois et al. (2010) and Marois et al. (2014) as Template LOCI (TLOCI). A conceptual illustration of the procedure used in our TLOCI pipeline to subtract a single region of a single target image is shown in Fig. 1.

In addition to the process of generating least-squares coefficients to subtract a region of the target image, Fig. 1 also illustrates the principle of forward modelling (FM), where the same coefficients are also used to model the algorithmic exoplanet throughput during the subtraction. As shown, a selection criterion, or aggressiveness, specifies that reference images be selected beyond a radial and/or azimuthal gap in wavelength and/or time, respectively, from the target image in order to minimize algorithmic exoplanet throughput losses. The highest aggressiveness of 1 (using dimensionless units calculated from relative aperture photometry between the target and reference image models) would allow using any reference image in the full target sequence, including the target image. In contrast, the lowest aggressiveness of 0 would not allow any reference images in the target sequence to be selected. Aggressiveness is also illustrated in Fig. 2. In this figure, a range of times and wavelengths around the target image are not selectable because, if used, a planet present in the target image would be partially subtracted by the same planet in the reference images. This un-selectable region of Fig. 2, illustrated by a white color, would be either larger or smaller for a lower or higher value of user-defined aggressiveness, respectively. The size of this un-selectable region also depends on the separation in the image, assumed exoplanet spectrum, wavelength, and parallactic angle relative to the target image. The target image is highlighted in the center of this region with a correlation of one. Then, reference images are chosen from the remaining sample based on correlation, also illustrated in Fig. 2. Generally, the most correlated reference images are near the target image time or wavelength, although this depends on the temporal stability and chromaticity of the sequence. A more correlated set of reference images (i.e., those with a more similar speckle pattern to the target image) provides a better basis for the least-squares algorithm and will yield a deeper final contrast without changing any other free parameters (Marois et al., 2014). We will examine the chromatic and temporal dependence of these correlations further in §3.3. Additional details are described below about the setup and new modifications of our TLOCI pipeline.

To help define, identify, and remove uncorrelated images from each dataset, we will introduce a new free parameter into our PSF subtraction pipeline called a correlation cut. For a single target image, the value of a correlation cut will define the minimum level of correlation that is rendered acceptable in a reference image used to subtract the target image speckle pattern. A given target image/region will be discarded from the subtracted sequence if there are not at least $N_{\text{ref}}$  reference images above the correlation cut value, all of which must pass a selection criterion, where $N_{\text{ref}}$  and aggressiveness are also free parameters. For example, if correlation cut = 0.5 and $N_{\text{ref}}$ = 10 but there are only five reference images that are greater than 50% correlated to the target image using a moderate aggressiveness, the target image/region will be discarded from the sequence. The target image would be discarded in this case even if, e.g., one of those five reference images is 99% correlated to the target image. Thus, it is clearly important to use a physically motivated correlation cut so that a maximal number of target images in an observing sequence can be subtracted while still rejecting a maximal number of uncorrelated reference images. Marois et al. (2010) and subsequent papers have shown that using correlated reference images is crucial to prevent unnecessary noise propagation in the covariance matrix inversion step of a least-squares/principal component analysis subtraction algorithm. The correlation cut values are dimensionless by definition and have a minimum and maximum of 0 and 1, respectively. Analysis in appendix A inform the decision to use the following correlation cut values for the four different radial sections of an image, progressing from innermost to outermost separation: 0.8, 0.74, 0.67, and 0.59. Although the specific angular separations of these regions varies for each target based on our method of measuring the pixel scale using the satellite spots (§3.1), the average separation and standard deviation across all targets is: 104$\pm$2 mas $<\rho<$ 286$\pm$4 mas, 286$\pm$4 mas $<\rho<$ 468$\pm$7 mas, 468$\pm$7 mas $<\rho<$ 649$\pm$10 mas, and $\rho>$ 649$\pm$10 mas. A graphical illustration of these zones is illustrated later in §3.2. Note that we chose to fix the radial separation of subtraction regions based on pixel separation, which is why the derived angular separation varies slightly from target to target. In pixels, the radial separations for each annulus for all targets are $7<\rho<19$, $19<\rho<32$, $32<\rho<44$, and $\rho>44$.

Based on the definitions and work in Marois et al. (2010) and Johnson-Groh (2016), for all subsequent analyses in this paper we use an L0 or T6 input spectrum (J02281101+2537380 and J02281101+2537380 from Burgasser et al. 2008 and Burgasser et al. 2004, respectively, subsequently convolved and binned to the appropriate CHARIS spectral resolution), an aggressiveness of 0.5 and 15 reference images per least-squares subtraction. Gerard & Marois (2016) found, on GPI datasets, that optimizing aggressiveness and number of reference images yielded a negligible gain in contrast, and so here we use the same aggressiveness, number of reference images per subtraction, and input spectral types (to define a selection criterion and for subsequent forward modelling). We use a singular value decomposition (SVD) cut-off value for the covariance matrix inversion of $1\times 10^{-4}$, as in Gerard & Marois (2016), where a maximum value of 1 would set all the LOCI coefficients to zero. Although this is a low value compared to other algorithms (e.g. Currie et al., 2018b), note that there are many degeneracies between the number of reference images used, reference image correlation, and SVD cut-off (Marois et al., 2010); our approach here uses a lower SVD cut-off but compensates by selecting more correlated reference images.

In a new addition to our TLOCI pipeline, the central region of each subtraction region is defined by both the parallactic angle and wavelength at any given timestamp, and so the location of a single region will both rotate with time and magnify with wavelength over a single observing sequence, similar to the approach in Wang et al. (2015). These definitions are used so that the subtraction regions line up after de-rotation and de-magnification, which allows the generation of a FM map for each target. The LOCI optimization regions are defined for each subtraction region as the other two azimuthal regions at the same radial separation. This subtraction/optimization region geometry ensures that the least-squares coefficients are never generated from the subtraction zone where a planet would be. S/N maps are calculated using back-rotated images, as in Gerard & Marois (2016). A few example data products of a single target sequence from our PSF-subtraction pipeline are shown in Figure 3, including a subtracted collapsed data cube, corresponding S/N map, and corresponding FM map.

For each target we also run a bootstrapping procedure, illustrated in Figures 4 and 5.

As in Marois et al. (2010) and Marois et al. (2014), we implement a FM algorithmic throughput correction procedure in each subtraction region of the image, including a model for the radial variation of the algorithmic throughput within each region. In addition to the center, four additional FM throughput corrections using the same LOCI coefficients are computed at smaller and larger separations within each subtraction region. A fifth order polynomial is then fit to the radial variation of throughput correction values to generate a continuous function across the region. As a result of this throughput correction procedure, a bootstrapped planet at any location in the image (i.e., not necessarily at the radial or azimuthal center of a subtraction region), is extracted at the input value and does not need any additional throughput correction. This is illustrated in Fig. 5. The simulated planets were injected throughout the image, one per subtraction region, such that none lies at the radial or azimuthal center of the region, as in the left panel of Fig. 4. Almost all recovered fluxes match the input value within their expected scatter, materialized by the 5$\sigma$ contrast value at the location of the planet.

3.3 Chromaticity Analysis

As discussed in §1, in a perfectly achromatic system, the larger “lever arm” of wavelength coverage over the JHK filter of CHARIS should provide a better algorithmic exoplanet throughput than a narrower bandpass, and so if both cases provide the same level of speckle subtraction, the final throughput-corrected contrast should be better in the former case. However, in reality broadband systems are not achromatic, and the level of wavefront decorrelation, or chromaticity (i.e., separate from wavelength magnification), will ultimately determine the limits of speckle subtraction via SDI processing (e.g., Marois et al., 2008). In this section we will measure the chromaticity in our CHARIS datasets, which will then inform the SDI speckle subtraction limits in §4.

A generalization of Figure 2 extends to Figure 6, which shows the median correlation of all of our targets as a function of time and wavelength.

Fig. 6 is generated as follows:

For a single target image at a single time, wavelength, and user-defined separation, spectral type, and aggressiveness, a selection criterion is applied to all of the images in the oberving sequence; images that do not meet this criterion are discarded. 2. 2.

Correlation of the target image with every remaining image from step 1 is computed, normalized to unity using the standard deviation in both the target image and reference images. Steps 1-2 are illustrated in Fig. 2. 3. 3.

In order to create a common zero-point for every subsequent target image (i.e., so the central index of Fig. 6 is zero for every target image), each correlation value from step 2 is saved in a two dimensional array of size $2i_{t}\times 2j_{t}$, where $i_{t}$ and $j_{t}$ are the total number of time stamps and wavelength slices, respectively (i.e., so both the negative and positive indices can be shown in Fig. 6). The $x$ and $y$ indices that fill a quarter of the total area of the two dimensional array for each correlation value from step 2 are $(2i_{t}-i_{\text{target}})$ and $(2j_{t}-j_{\text{target}})$, respectively, where $i_{\text{target}}$ and $j_{\text{target}}$ are the indices of time and wavelength for the target image. This index definition ensures that the target image is in the center, or zero point, of each $2i_{t}\times 2j_{t}$ array. Both the region removed by the selection criterion and the other three quarters of the array remain empty and are filled with “not a number” (NaN) placeholders. 4. 4.

Steps 1 - 3 are repeated over time and wavelength for every target image within a single observing sequence, generating a vector of two dimensional arrays of size $i_{t}j_{t}\times 2i_{t}\times 2j_{t}$. The autocorrelation of each target image is in the center of each $2i_{t}\times 2j_{t}$ array. 5. 5.

A median is taken across the zeroth axis of the vector in step 4 (i.e., the axis of dimension $i_{t}j_{t}$), ignoring NaN values, to generate a two dimensional array of size $2i_{t}\times 2j_{t}$, representing the median correlation between the target image and every other image in the sequence. A slice along the $x$ direction of this array at $y=j_{t}$ (the zero point of the wavelength direction) reveals the median correlation of the target image (at $x=i_{t}$) vs. time, whereas a slice along the $y$ direction of this array at $x=i_{t}$ reveals the median correlation vs. wavelength. 6. 6.

Steps 1-5 are repeated over each user-defined spectral type. 7. 7.

Steps 1-6 are repeated over each user-defined separation. 8. 8.

Steps 1-7 are repeated for each target, computing the median and standard deviation of correlation over time and wavelength for each separation.

The outcome of the above procedure is illustrated in Figure 6, whereas steps 1-7 are illustrated in Fig. 7, for which we used a separate dataset from our sample, the $\kappa$ And data from Currie et al. (2018b). Rather than a single target in our sample, the Currie et al. (2018b) dataset was chosen for Fig. 7 to understand how the conclusions from Fig. 6 change with higher quality AO correction.

The main empirical conclusion from Figures 6 and 7 is that SCExAO/CHARIS datasets are generally more correlated as a function of time than as a function of wavelength, most notably at small and large separations. This is particularly important at small separations, where PSF subtraction algorithms generally underperform compared to the rest of the image. For a given difference in parallactic angle and/or wavelength between a target and reference image, the respective azimuthal and/or radial movement of an exoplanet in between the two images is larger at larger separations. In other words, For a given aggressiveness, more reference images near the target image, both in time and wavelength, have to be discarded at lower separations compared to higher separations. Although a longer wavelength lever arm would help overcome this problem, Figures 6 and 7 illustrate that this would generally add unnecessary noise in the covariance matrix inversion. A more correlated image at the same wavelength and different time will be selected instead. Later in §4 we will test the predictions from Fig. 6 on our PSF subtraction pipeline.

With that said, at small separations in Fig. 6 and at all separations in Fig. 7, reference images are actually more correlated at the closest wavelengths than at the closest times. Although these images will be preferentially selected for SDI processing, this approach will not necessarily reach the deepest (throughput-corrected) final contrasts. Selecting a closer wavelength will reach a better non-throughput-corrected contrast but will also cause more self-subtraction, therefore requiring a larger throughput correction. Thus, the optimal reference selection based on final, throughput-corrected contrast will ultimately depend on optimizing the trade off between correlation and self-subtraction, as in Gerard & Marois (2016). Additionally, although the chromaticity is still worse than the temporal-stability in Fig. 7, note that there is a significant improvement in both temporal stability and chromaticity between Fig. 6 and 7. Although the improvement from temporal stability is most likely from improved AO performance (strehl ratios for the Currie et al. 2018b dataset are around 0.92), the origin of chromaticity improvement is not clear and will require further investigation and analysis, beyond the discussion below.

Although the specific origin of the observed chromaticity is beyond the scope of this paper, we discuss a few possibilities here that could explain these results. In general, Fresnel propagation generates a chromatic evolution of the wavefront from any out-of-pupil plane optics, particularly for optics that are transmissive and/or near the focal plane (e.g., Marois et al., 2008), as well as from atmospheric scintillation (Guyon, 2005), although both of these effects are expected to be very weak at the current level of raw contrasts (but see Madurowicz et al. 2018, who show that in some cases scintillation effects are observed in GPI data). Additional sources of the observed behavior may be algorithmic in nature, arising from either numerical interpolation errors and/or DRP extraction errors. These numerical interpolation errors are wavelength-dependent when aligning data cubes, since, as discussed in §3.1, the full CHARIS bandpass goes from near the Nyquist sampling limit in J band to more oversampled in K band. Although numerical interpolation errors are more likely to disproportionally effect the J band (from under-sampling effects) while DRP extraction errors are more likely to affect the K band (from increased thermal background effects), further investigation is needed to understand the relative impact of either effect. At small separations there are two additional possible origins of the observed chromaticity: the atmospheric dispersion corrector (ADC) and the FPM IWA. The ADC does not sufficiently offset differential atmospheric tip/tilt to center the FPM on the star at the red edge of the K band. For the FPM IWA, because our chromaticity analysis uses magnified data cubes to align speckles as a function of wavelength, the IWA on magnified data cubes varies between about 2 and 4 $\lambda/D$ at the red and blue ends of the bandpass, respectively, which could explain the observed chromaticity in the inner-most radial zone.

4 Contrast Curves

Our 5$\sigma$ contrast curves are shown in Fig. 8. Final contrasts are calculated using the standard deviation at each separation in a 3 pixel-wide annulus on the collapsed, back-rotated images. Raw contrasts are calculated on a median image of the stack of the registered, magnified, flux-normalized, high-pass-filtered images (i.e., the same image preparation used before the correlation analysis in addition to a median stack of the full sequence). Also note that the frames used to calculate raw contrast are not de-rotated to align the parallactic angle, thus keeping quasi-static aberration aligned during the median collapse.

Although other deeper contrasts have been reached with CHARIS on more stable datasets (e.g., Currie et al., 2018b; Rich et al., 2018), the main purpose of this paper is instead to provide a detailed correlation analysis (§3.3) that can be used as a reference for future CHARIS observations (e.g., considering the impact of SDI in broadband mode).

With that said, Fig. 8 a and b does show that T-type subtractions generally reach up to $\sim$50% deeper contrasts than the corresponding L-type subtractions for the same target, illustrating that SDI is improving contrast. Contrasts in a single slice are shown in addition to the final broadband contrasts to control for the effects of noise averaging in comparing an L- vs. T-type reduction. In applying a weighted average using the assumed planet spectral type to collapse the final subtracted cube, noisier frames that lie further away from the bright part(s) of the spectrum are given less weight to the final collapsed cube (Marois et al., 2014), and so looking only at a single wavelength slice at the peak of the H band T-type spectrum removes this affect. Although other wavelengths will not show as much of a performance gain further away from the T-type spectrum peak, a deeper contrast reached by a T-type vs. L-type, at any wavelength, illustrates that SDI is being used to improve contrast. This improvement occurs because, within the 15 reference images that are chosen to subtract a target image at a given time and wavelength, some are at the same time but different wavelength as the target image. These results are consistent with our analysis in §3.3 which showed that chromaticity generally prevents using a reference image for SDI with a large wavelength lever arm, but that reference images closer in wavelength were often more correlated than images closer in time. For a T-type subtraction, the relative lack of long dark shadows in the spectral (i.e., radial) direction of the right panel of Fig. 3 is also consistent with these results; selecting the reference images closest in wavelength, especially for a T-type subtraction, are still often more correlated than the closest images in time.

Along these lines, because the main conclusions of Figures 6 and 7 indicate that reference images at a large wavelength difference from the target image are unlikely to be selected (i.e., separate from the hypothesis of whether or not SDI is improving contrast), we carried out the exact same PSF subtraction procedure as described in §3.2 on all seven targets using only the H band wavelengths within the broadband data cubes (i.e., slices 6 - 13). The idea here is to see what level of contrast improvement is lost by removing the broadband feature of CHARIS while still using SDI in a standard $\sim$20% bandpass. Note that using the same parameters as in §3.2 (e.g., N${}_{\text{ref}}$=15) means that there are less available reference images as a function of wavelength to select for the H band-only sequence. However, because reference images are selected by correlation, our chromaticity analysis in §3.3 and Figure 6 has determined that images at these wider wavelength differences are generally not selected, and so we expect this discrepancy to have a negligible effect on final contrast.

Fig. 9 illustrates the contrast difference on a single slice in the center of H band for all of our targets between the full broadband mode and the H band mode. Note that comparing the collapsed cubes between the broadband and H-only subtractions would introduce a systematic bias; in the final wavelength collapse the broadband dataset would average more slices than the H-only dataset, producing a deeper final contrast for the former even if the individual subtracted cubes are the same for both sequences. In general, all of our April 2017 observations show little to no gain using the broadband dataset, as expected from Fig. 6. Interestingly, the January 2018 observations, which reach the deepest raw and final contrasts in Fig. 8, instead consistently show a $\sim$0.5 - 1 magnitude gain from using the broadband mode. Also when using the broadband mode, in a few cases at separations less than about 200 mas, the L-type subtractions appear to decrease in contrast by up to a magnitude.

Ultimately, the origin of this observed chromaticity that may be improving and/or degrading contrast will need to be better understood before a detailed analysis leading to mitigation strategies can be done. As discussed in §3.3, many different possible sources of chromaticity could reproduce the results in Figures 6, 7, and 9. These sources may be algorithmic and/or optical in nature, and could be quasi-static and/or time-variable. However, independent of the origin, in this paper we illustrate that chromaticity in our broadband datasets is limiting the final achievable contrast, thus motivating the need for further work towards more achromatic ExAO systems.

5 Conclusions

In this paper, we have presented a detailed correlation analysis, data processing architecture, and final contrast curves of our SCExAO/CHARIS broadband datasets to the community. Our main findings are as follows:

A correlation analysis of the SCExAO/CHARIS system across all seven of our targets shows that for small ($\lesssim$300 mas) and large ($\gtrsim$650 mas) separations, chromaticity across the full JHK bandpass is worse than temporal stability across the full observing sequence (Fig. 6). This result suggests that at these separations, the large wavelength “lever arm” of our CHARIS broadband datasets could not be used advantageously to improve SDI post-processing (i.e., slices in J band can generally not be used to subtract slices in K band, and vice versa). 2. 2.

With that said, we find that our T-type PSF subtractions generally reach up to $\sim$50% deeper final contrasts than our L-type subtraction of the same target (Fig. 8), illustrating that more correlated wavelength slices close to the target image wavelength are being used to improve contrast. This is consistent with the findings of point 1 above (i.e., Fig. 6); together, these results illustrate that improved contrast from SDI processing of our targets comes from selecting reference slices closer in wavelength to the target image where chromaticity is less significant. 3. 3.

We present a modified framework for exoplanet algorithmic throughput correction using forward modelling that is designed for the CHARIS broadband mode (§3.2). As a result, simulated planets are recovered in good agreement with their input values even when close to the 5$\sigma$ noise floor (Fig. 5). 4. 4.

PSF-subtracted sequences using only the H-band slices within the broadband dataset generally reach final contrasts that are similar to the final contrasts of the same slices in the broadband sequences (Fig. 9), validating the claim made in points 1 and 2 above.

Beyond the work presented in this paper, as long as ExAO systems remain limited by temporal stability (i.e., preventing subtraction to the photon noise limit using ADI and/or RDI post-processing, thus motivating the additional need for SDI-processing), chromaticity will also limit the deepest contrasts that can be reached. Looking towards the future, coherent differential imaging (CDI; e.g., Guyon, 2004; Bordé & Traub, 2006; Baudoz et al., 2006; Give’On et al., 2007; Serabyn et al., 2011; Sauvage et al., 2012; Jovanovic et al., 2018, and references therein) is a complimentary approach towards reaching deeper contrasts on ExAO systems. CDI can operate in a narrow bandpass (i.e., in contrast with SDI) and is therefore, in principle, not limited by chromaticity. However, detailed characterization of colder/lower mass exoplanet atmospheres will continually push the ExAO community towards reaching deeper broadband contrasts, in which case chromaticity will ultimately become a fundamental limitation to this goal. The future push towards reaching deeper contrasts will therefore also require a push towards implementing more achromatic broadband systems.

Acknowledgements

The data in this paper was obtained through 2017B open use Canadian and US Gemini exchange time (Gerard, 2017). We gratefully acknowledge research support of the Natural Sciences and Engineering Council of Canada through the Postgraduate Scholarships-Doctoral discovery grant (DG), and Technologies for Exo-Planetary Science Collaborative Research and Training Experience programs. The development of SCExAO was supported by the JSPS (Grant-in-Aid for Research #23340051, #26220704 #23103002), the Astrobiology Center of the National Institutes of Natural Sciences, Japan, the Mt Cuba Foundation and the directors contingency fund at Subaru Telescope. CHARIS was built at Princeton University under a Grant-in-Aid for Scientific Research on Innovative Areas from MEXT of the Japanese government (# 23103002). We thank Jason Wang and Jean-Baptiste Ruffio for helpful discussions and suggestions about forward modelling. The authors thank the anonymous referee for his or her comments and suggestions that have significantly improved this manuscript. The authors acknowledge and support the cultural and spiritual importance of the summit of Mauna Kea to the Hawaiian native community. We are grateful to share and respect this land with the community.

Appendix A Additional Correlation Analysis

Additional correlation analysis informed our choice of values for a correlation cut, as discussed in §3.2. Still using the HD 23267 Jan. 7 dataset, frame 17, slice 10 as an example (as in Fig. 2), Fig. 10 illustrates the trends in correlation of selected reference images at either the same time or same wavelength (i.e., i = 17 or j = 10, respectively). The top two panels show the correlation trends of the 15 most correlated images as defined in §3.2, N${}_{\text{ref}}$=15) at large separations as a function of time (a) or wavelength (b) at large separations, using an aggressiveness of 0.5, and using an assumed L0 input spectrum. Uncorrelated frames can be identified as vertical “streaks” in these images (e.g., the water bands at j = 5 and j = 14, as in Fig. 2). The chosen correlation cut value at this separation should reject these uncorrelated frames from the sequence in order to prevent unnecessary noise amplification in the least-squares algorithm (Marois et al., 2010). Because we require that all 15 reference images lie above our defined correlation cut value, the bottom panels of Fig. 10 show a one dimensional slice of the top two panels: the correlation of reference image number 14 (i.e., the 15th most correlated selectable reference, using a zero-index counting system) as a function of time (c) and wavelength (d) for both spectral types and all separations. Similarly, at all separations and spectral types we would like to discard the frames where there are “dips” in these plots in order to optimally reject uncorrelated images and avoid propagating them through the least-squares algorithm. Fig. 10 c and d also suggest that larger separations may be less correlated overall and that a different correlation cut value will be needed for each separation. Note that the low correlation at larger separations may be a result of wavefront decorrelation at those spatial frequencies and/or increased background noise because the speckles at those separations are detected at a lower S/N. At a given separation, using a correlation cut value that is too high will reject too many images at that separation and return only a minimal number of subtracted frames to the final collapsed cube; the “optimal” correlation cut value should balance rejecting too many images with propagating too much noise through the least-squares algorithm.

The extension of Fig. 10 c and d is to calculate the median behavior over a full sequence and for all targets. Thus, for a single target at a given separation, in Fig 11 we compute the median and robust standard deviation for the 15th most correlated reference image over all target images (time and wavelength) and input spectral types. Fig. 11 illustrates that, for most targets, the 15th most correlated image decreases in correlation with separation, motivating the use of different correlation cut values at each separation. The horizontal dashed lines in Fig. 11 illustrate our chosen correlation cut values at each separtion. We set these values to generally lie near the bottom of the corresponding 1$\sigma$ error bars over all targets. Using this definition to define the correlation cut values will balance the rejection of bad frames relative to the median correlation value without discarding too many frames from the sequence. In general, the optimal correlation cut value will vary for each least-squares subtraction and is degenerate with SVD cutoff and number of reference images (Marois et al., 2010); this will ultimately require a more optimized approach (e.g. Gerard & Marois, 2016), but is beyond the scope of this paper.