Data Reduction Pipeline for the CHARIS Integral-Field Spectrograph I: Detector Readout Calibration and Data Cube Extraction

TL;DR

This paper introduces a comprehensive data reduction pipeline for the CHARIS integral-field spectrograph, detailing calibration, extraction algorithms, and noise modeling to produce high-quality, calibrated data cubes.

Contribution

It presents a novel pipeline with advanced PSF modeling, noise removal, and multiple extraction methods, improving data quality for high-contrast spectroscopy.

Findings

Achieved ~5% residuals with $ ext{chi}^2$ extraction

Enhanced noise removal by modeling correlated read noise

Produced deconvolved, uncertainty-included data cubes

Abstract

We present the data reduction pipeline for CHARIS, a high-contrast integral-field spectrograph for the Subaru Telescope. The pipeline constructs a ramp from the raw reads using the measured nonlinear pixel response, and reconstructs the data cube using one of three extraction algorithms: aperture photometry, optimal extraction, or $\chi^2$ fitting. We measure and apply both a detector flatfield and a lenslet flatfield and reconstruct the wavelength- and position-dependent lenslet point-spread function (PSF) from images taken with a tunable laser. We use these measured PSFs to implement a $\chi^2$-based extraction of the data cube, with typical residuals of ~5% due to imperfect models of the undersampled lenslet PSFs. The full two-dimensional residual of the $\chi^2$ extraction allows us to model and remove correlated read noise, dramatically improving CHARIS' performance. The $\chi^2$ extraction produces a data cube that has been deconvolved with the line-spread function, and never performs any interpolations of either the data or the individual lenslet spectra. The extracted data cube also includes uncertainties for each spatial and spectral measurement. CHARIS' software is parallelized, written in Python and Cython, and freely available on github with a separate documentation page. Astrometric and spectrophotometric calibrations of the data cubes and PSF subtraction will be treated in a forthcoming paper.

Tables

Table 1. Table 1 Basic CHARIS Parameters

Parameter	Value
Detector	$2048\times 2048$ Hawaii2-RG
# of Lenslets	$135\times 135$
Lenslet Size	16.4 mas
Field-of-View	$2{.}^{\prime \prime }2\times 2{.}^{\prime \prime }2$
Wavelength Coverage	1.15 - 2.38 $\mu$m
Microspectrum Length	$\sim$30 pixels
$R=\lambda /\delta \lambda$ (2 pixels)	$\sim$20 (low-res), $\sim$75 (high-res)
Available Modes	$J$, $H$, or $K$ at $R\sim 75$
	$J+H+K$ or ND$\dagger$$\dagger$Filter is ND3 (${10}^{-3}$ transmission) from 1.15 to 2.4 $\mu$m, opaque at other wavelengths. at $R\sim 18$

Table 2. Table 2 Extraction Software Parameters

Parameter	Default	Allowed Values	Description
Ramp Parameters
$\mathrm{𝚛𝚎𝚊𝚍}\_\mathrm{𝟶}$	1	$[1,N_{\mathrm{reads}}-1]$	First read in the ramp to use.
$\mathrm{𝚛𝚎𝚊𝚍}\_𝚏$	$\mathrm{𝙽𝚘𝚗𝚎}$	$[\mathrm{𝚛𝚎𝚊𝚍}\_\mathrm{𝟶}+1,N_{\mathrm{reads}}]$ or $\mathrm{𝙽𝚘𝚗𝚎}$	Last read to use. If $\mathrm{𝙽𝚘𝚗𝚎}$, $\mathrm{𝚛𝚎𝚊𝚍}\_𝚏$ will be the final read.
$\mathrm{𝚐𝚊𝚒𝚗}$	2	$[0,\infty )$	Detector gain, $e^{-}$/count, used to compute photon noise.
$\mathrm{𝚗𝚘𝚒𝚜𝚎𝚏𝚊𝚌}$	0	$[0,\mathrm{\hspace{0.17em}1})$	Additional fractional error in count rate to account for systematic errors in
			PSFlet models. Suggested values: 0.05 for ${\chi }^2$ extraction, 0 otherwise.
$\mathrm{𝚜𝚊𝚟𝚎𝚛𝚊𝚖𝚙}$	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Save the ramp as its own image?
Calibration Parameters
$\mathrm{𝚌𝚊𝚕𝚒𝚋𝚍𝚒𝚛}$	…	string file path	Directory containing the files created by $\mathrm{𝚋𝚞𝚒𝚕𝚍𝚌𝚊𝚕}$
$\mathrm{𝚋𝚐𝚜𝚞𝚋}$	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Subtract a background/thermal ramp? If True, buildcal must have been
			called with background ramp(s). Should be False for ND frames.
$\mathrm{𝚖𝚊𝚜𝚔}$	$\mathrm{𝚃𝚛𝚞𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Mask bad pixels?
$\mathrm{𝚏𝚕𝚊𝚝𝚏𝚒𝚎𝚕𝚍}$	$\mathrm{𝚃𝚛𝚞𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Apply the pixel and lenslet flatfields?
$\mathrm{𝚏𝚒𝚝𝚜𝚑𝚒𝚏𝚝}$	$\mathrm{𝚃𝚛𝚞𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Fit for a subpixel shift of the microspectra?
			If $\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝚋𝚞𝚒𝚕𝚍𝚌𝚊𝚕}$ must have computed oversampled PSFlets.
Extraction Parameters
$𝚁$	30	$[1,\infty )$	Spectral resolution $R={\left(\delta \ln \lambda \right)}^{-1}=\lambda /\delta \lambda$ of the extracted cube
			Recommended: 30 for low res, 100 for high res
$\mathrm{𝚖𝚎𝚝𝚑𝚘𝚍}$	$\mathrm{𝚕𝚜𝚝𝚜𝚚}$	$\mathrm{𝚕𝚜𝚝𝚜𝚚}$, $\mathrm{𝚘𝚙𝚝𝚎𝚡𝚝}$,	Method used to extract cube: ${\chi }^2$ extraction, optimal extraction,
		$\mathrm{𝚊𝚙𝚙𝚑𝚘𝚝𝟹}$, $\mathrm{𝚊𝚙𝚙𝚑𝚘𝚝𝟻}$	or aperture photometry across 3 or 5 pixels.
$\mathrm{𝚛𝚎𝚏𝚒𝚗𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Iterate once to remove crosstalk? Only used if method is $\mathrm{𝚕𝚜𝚝𝚜𝚚}$.
$\mathrm{𝚜𝚞𝚙𝚙𝚛𝚎𝚜𝚜𝚛𝚗}$	$\mathrm{𝚃𝚛𝚞𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Use the method of Section 5.3 to remove read noise?
$\mathrm{𝚖𝚒𝚗𝚙𝚌𝚝}$	70	$(0,\mathrm{\hspace{0.17em}100})$	Minimum percentage of pixels used to estimate read noise (only
			used if $\mathrm{𝚜𝚞𝚙𝚛𝚎𝚜𝚜𝚛𝚗}$ is $\mathrm{𝚃𝚛𝚞𝚎}$).
$\mathrm{𝚏𝚒𝚝𝚋𝚔𝚐𝚗𝚍}$	$\mathrm{𝚃𝚛𝚞𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Fit and remove an undispersed background lenslet-by-lenslet?
$\mathrm{𝚜𝚖𝚘𝚘𝚝𝚑𝚊𝚗𝚍𝚖𝚊𝚜𝚔}$	$\mathrm{𝚃𝚛𝚞𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Interpolate over bad lenslets to display a smooth cube?
			The inverse variance of the bad measurements remains zero.
$\mathrm{𝚜𝚊𝚟𝚎𝚛𝚎𝚜𝚒𝚍}$	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$	$\mathrm{𝚃𝚛𝚞𝚎}$, $\mathrm{𝙵𝚊𝚕𝚜𝚎}$	Save the 2D residual from the ramp? Only used if method is $\mathrm{𝚕𝚜𝚝𝚜𝚚}$.
$\mathrm{𝚖𝚊𝚡𝚌𝚙𝚞𝚜}$	$\mathrm{𝙽𝚘𝚗𝚎}$	$[1-N_{\mathrm{cpus}},N_{\mathrm{cpus}}]$, or $\mathrm{𝙽𝚘𝚗𝚎}$	Maximum number of threads to use. If $\mathrm{𝙽𝚘𝚗𝚎}$, use all threads.
			If negative, it is the number of threads to reserve.

Table 3. Table 3 Suggested Parameter Settings

Parameter	Fast & Rough	Slower & Better
	Reduction	Reduction
read_0	1	1
read_f	None	None
gain	2	2
noisefac	0	0.05
bgsub	$\dagger$$\dagger$See Sections 5.4 and 6 for a discussion. Never set both bgsub and fitbkgnd to True.False	$\dagger$$\dagger$See Sections 5.4 and 6 for a discussion. Never set both bgsub and fitbkgnd to True.False
mask	True	True
flatfield	True	True
fitshift	False	True
R	$\dagger\dagger$$\dagger\dagger$30 in low-res mode, 100 in $J$, $H$, or $K$.30 or 100	$\dagger\dagger$$\dagger\dagger$30 in low-res mode, 100 in $J$, $H$, or $K$.30 or 100
method	optext	lstsq
refine	…	True
suppressrn	False	True
minpct	…	70
fitbkgnd	$\dagger$$\dagger$See Sections 5.4 and 6 for a discussion. Never set both bgsub and fitbkgnd to True.False	$\dagger$$\dagger$See Sections 5.4 and 6 for a discussion. Never set both bgsub and fitbkgnd to True.True
smoothandmask	True	True

Equations

$$R=\lambda /\delta \lambda ={\left(2\frac{\mathrm{d}\ln \lambda }{\mathrm{d}x}\right)}^{-1},$$

$${\chi }_i^2=\sum _{\mathrm{reads}j=1}^N\frac{{\left(j\delta t\cdot C_i+b_i-F_{ij}\right)}^2}{{\sigma }^2},$$

$$C_i\delta t=\frac{12}{N^3-N}\sum _{j=1}^N\left(j-\frac{N+1}{2}\right)F_{ij},$$

$$\frac{{\mathrm{SNR}}_{\mathrm{UTR}}}{{\mathrm{SNR}}_{\mathrm{CDS}}}=\sqrt{\frac{N(N+1)}{6(N-1)}}.$$

$$y=A_i\exp \left[-t/t_0\right],$$

$${\chi }^2[x,b]=\sum _j\frac{{\left(F_j-f[x\cdot j]+b\right)}^2}{{\sigma }^2},$$

$$\left\{\begin{array}{c}{x}_{ij}[\lambda ]\\ y_{ij}[\lambda ]\end{array}\right\}=\sum _{m=0}^3\sum _{n=0}^{3-m}\left\{\begin{array}{c}{a}_{nm}[\lambda ]\\ b_{nm}[\lambda ]\end{array}\right\}{i}^m{j}^n.$$

$${\lambda }_i=LAM\_MIN\times \exp \left[i\cdot \mathtt{DLOGLAM}\right]$$

Full text

1111affiliationtext: Astrobiology Center of NINS, Osawa, Mitaka, Tokyo, Japan

1. Introduction

Due to the advent of large format detectors, integral-field spectrographs (IFSs) have become an increasingly popular class of astronomical instrumentation. IFSs are hybrids of traditional imaging cameras and slit spectrographs: they obtain a spectrum from each spatial element in a two-dimensional field-of-view for an $(x,y,\lambda)$ data cube. The first realization of an IFS used a bundle of fibers to create a pseudo-slit (Vanderriest, 1980), while TIGER (Bacon et al., 1995) was the first IFS to use a lenslet array. Modern IFSs generally use either fiber bundles (MaNGA, Drory et al., 2015) or image slicers (SINFONI, MUSE, NIRSpec Eisenhauer et al., 2003; Bacon et al., 2010; Bagnasco et al., 2007) to rearrange the field-of-view into a long pseudo-slit, or lenslet arrays (OSIRIS, Larkin et al., 2006) to reimage spatial elements into small spots suitable for dispersion.

IFSs have become especially popular tools for high-contrast imaging. Diffraction speckles have a different chromatic behavior from astrophysical sources; an IFS data cube can exploit this to achieve higher contrasts (Sparks & Ford, 2002; Marois et al., 2014). An IFS also naturally enables the extraction of a planet or brown dwarf’s spectrum, providing a probe of the object’s temperature, chemistry and gravity (McElwain et al., 2007; Barman et al., 2011; Konopacky et al., 2013; Hinkley et al., 2013; Currie et al., 2014). IFSs combined with second-generation adaptive optics systems are now operational on Gemini South (GPI, Macintosh et al., 2008), the VLT (SPHERE, Claudi et al., 2008), and Palomar (Project 1640, Hinkley et al., 2011). These new high-contrast instruments have recently discovered and characterized the low-mass companion to 51 Eri (Macintosh et al., 2015). Future NASA mission studies such as Exo-C (Stapelfeldt et al., 2014) and the WFIRST Coronagraph Instrument have baselined high contrast IFSs as their science cameras (Spergel et al., 2015; McElwain et al., 2016).

Data reduction and processing for IFSs has long presented problems. The reduction pipeline for GPI (Perrin et al., 2014, 2016) is an ongoing, years-long effort partially built on legacy software from OSIRIS (whose pipeline also remains, to some degree, a work in progress). This is the result of many complexities inherent in IFS data. There are now two flatfields (one for the detector and one for the illumination of the fibers, lenslets, or sliced image plane). For a lenslet-based IFS, the point-spread function (PSF) of the input optics and of the lenslets are both important. The finite size of the lenslet PSFs means that neighboring spectra partially overlap one another; this must be corrected or accounted for during the extraction. An IFS requires the wavelength and spectrophotometric calibrations of a spectrograph as well as the astrometric calibration of an imager.

This paper presents the data reduction pipeline for the CHARIS IFS on the Subaru telescope. CHARIS, the Coronagraphic High Angular Resolution Imaging Spectrograph, is a lenslet-based near-infrared IFS located on the Nasmyth platform behind the adaptive optics systems AO188 (Minowa et al., 2010) and SCExAO (Jovanovic et al., 2015b). Section 2 summarizes the design and properties of the IFS, while the rest of the paper presents the software that extracts the data cube. Section 3 discusses the construction of the pixel-by-pixel count rates from a sequence of raw reads, Section 4 discusses our calibration procedure and associated data products, and Section 5 presents our algorithms for extracting the data cube. Section 6 summarizes the software’s parameters and settings and Section 7 shows the software’s performance. We discuss and conclude with Section 8.

This software package constructs the data cube and its inverse variance from a sequence of CHARIS reads. We defer the necessary steps of image registration and spectrophotometric and astrometric calibration to a separate software package that is currently under development. These steps depend on the SCExAO system in front of CHARIS and the calibrations changed when SCExAO was rebuilt in July of 2016; they also rely on a system of induced satellite spots that still must be manually controlled by the SCExAO team (Jovanovic et al., 2015a). These elements of the software, in addition to algorithms for angular and spectral differential imaging, will operate only on the data cubes produced by this pipeline and will be presented in a follow-up paper.

2. The CHARIS Integral-Field Spectrograph

CHARIS is a new high-contrast IFS for the Subaru Telescope. Its scientific, conceptual, optical, and mechanical designs are summarized by McElwain et al. (2012); Peters et al. (2012), Galvin et al. (2014), and Peters-Limbach et al. (2013), respectively. Groff et al. (2016) summarize laboratory testing performed after CHARIS was built but before it was transported to the summit. This section briefly reviews the basic parameters and observing modes of CHARIS; we refer the reader to these other papers and to a forthcoming instrument paper for details.

CHARIS is a lenslet-based diffraction-limited spectrograph operating in the near-infrared. Table 1 summarizes its basic as-built properties. CHARIS uses one of two prisms behind the lenslet array to disperse the light from each lenslet into a $\sim$30-pixel-long microspectrum. The detector image consists of about $135\times 135$ of these microspectra, each containing the light incident on a single $16.4$ mas square lenslet, for a $\sim$2$.\!\!^{\prime\prime}2\times 2.\!\!^{\prime\prime}2$ field-of-view.

CHARIS offers five observing modes, three with a high-resolution prism and two with a low resolution prism. We measure its as-built spectral resolution $R$ using the definition

[TABLE]

with $x$ in pixels, so that the dispersion is a wavelength shift per two pixels (slightly larger than the full width at half maximum, or FWHM, of the lenslet PSFs). With this definition, $R$ with the high-resolution prism varies from $\sim$85 at the short end of $J$ to $\sim$65 in the middle of $H$ to $\sim$85 at the long end of $K$. The low-resolution prism has an $R$ that varies from about 18 to 22 across the $J$, $H$, and $K$ bands. The high-resolution mode uses either a $J$, $H$, or $K$ band filter, while the low-resolution prism may be used either with a broadband or a neutral density (ND) filter. The broadband filter has nearly unit transmission from 1.15 to 2.38 $\mu$m and sharp cutoffs toward both shorter and longer wavelengths. The ND filter is ND3 (10*-3* transmission) between 1.15 and 2.38 $\mu$m and opaque at other wavelengths; it is inside the CHARIS dewar to reduce the thermal $K$-band background. The ND filter is intended to allow observers to obtain unsaturated images of bright stars. CHARIS’ location at Subaru’s Nasmyth platform gives it better stability than GPI or Project 1640, while its extremely flat dispersion across the $J$, $H$, and $K$ bands is unique among high-contrast IFSs. CHARIS’ use of pinholes on the lenslet array and its relatively generous spacing of microspectra result in very low spectral cross-talk (Groff et al., 2016) and enable us to model and remove correlated read noise (Section 5.3).

CHARIS is controlled by a Linux-based software and is integrated into the software environment of the observatory. It has only three moving parts to be controlled during observations: a shutter, the five-slot filter wheel ($J$, $H$, $K$, broadband, and ND), and the three-position prism slider (low-resolution, high-resolution, and empty). The only other command is to reset and read out the detector. The rest of this paper presents the software for reconstructing the individual lenslet microspectra, i.e. the data cube, from these reads.

3. From the Reads to a Ramp

CHARIS’ detector is a Hawaii2-RG (H2RG), a HgCdTe CMOS device in which each pixel has its own amplifier. The pixel is read out by measuring the voltage across it relative to a reference voltage in the system. In CHARIS’ configuration, 32 readout channels each read pixels at a rate of 100 kHz. Including some down time as the readout proceeds to the next row of pixels, it takes 1.47 seconds to read out the full $2048\times 2048$ pixel array. Resetting the detector is done pixel-by-pixel, and also takes 1.47 seconds. H2RG detectors have some generic shortcomings, including persistence after exposure to a bright source (Smith et al., 2008) and $1/f$ read noise strongly correlated between readout channels (Moseley et al., 2010).

CHARIS always saves every one of its reads; this sequence of raw reads is then saved to disk. The first step of the data extraction is to fit for the count rate at each pixel from the raw reads. In the limit of a linear pixel response and the dominance of read noise over photon noise, the best-fit count rate $C_{i}$ in pixel $i$ may be derived from a $\chi^{2}$ fit to the sequence of reads $j$:

[TABLE]

where $\sigma^{2}$ is the variance from the read noise, $b_{i}$ is the pixel’s reset value, $j\delta t$ is the time from reset at which pixel $i$ was read out for the $j^{\rm th}$ time, and $F_{ij}$ is the number of counts in read $j$ of pixel $i$. Minimizing $\chi^{2}$ with respect to $C_{i}$ and $b_{i}$, we find that the best-fit count rate $C_{i}$ is given by a linear combination of the reads $j$, and the readout is called up-the-ramp (UTR, Fowler & Gatley, 1990):

[TABLE]

where $N$ is the total number of reads. Up-the-ramp, or a variant using only some of the available reads, is now commonly used to read out infrared arrays both on the ground and in space (Rauscher et al., 2007; Finger et al., 2008; Perrin et al., 2014; Dressel, 2017).

In the limit of uncorrelated read-noise, up-the-ramp improves on the signal-to-noise ratio of correlated double sampling (CDS), the normalized difference of the first and last reads, by a factor

[TABLE]

If there are only two reads, i.e. $N=2$, up-the-ramp and CDS are equivalent, and the ratio in Equation (4) is unity. It is never less than one, and increases asymptotically as $\sqrt{N/6}$, or $\sqrt{t_{\rm exp}/(9\,{\rm s})}$ assuming 1.5 seconds/read. When photon noise dominates, up-the-ramp is asymptotically noisier than correlated double sampling by about 23% (Fixsen et al., 2000). This may be corrected by dynamically choosing different weights for each read at each pixel (e.g. Robberto, 2014), but such an approach is not easily compatible with our nonlinear fit (Section 3.2). In the high signal-to-noise regime we are limited by the fidelity of our models of the microspectra rather than by photon noise.

We compute a variance on a given pixel’s count rate from the variance in the count rates of the reference pixels in its channel. We then add photon noise assuming our configured gain of 2 $e^{-}$/count to be correct. We have verified this gain by using a frame-to-frame scatter in count rate as a crude measure of shot noise. In addition, we allow the user to enforce an additional error equal to a fixed fraction (with a suggested value of a few percent) of the count rate seen by each pixel. This additional error accounts for our imperfect model of the microspectra and produces reduced $\chi^{2}$ values close to unity when fitting these models to the ramps; we add it in quadrature to the other errors. The variances are saved as inverse variances so that bad pixels may be masked by giving them an inverse variance of zero. Cosmic ray hits are rare in short CHARIS exposures, and we do not implement cosmic ray rejection within our ramps (Fixsen et al., 2000; Offenberg et al., 2001). Instead, we mask pixels in which a single read exceeds that pixel’s mean count rate by a factor of five and its mean read noise by a factor of ten. Less than 0.03% of pixels are masked in this way in a typical two minute exposure.

The CHARIS software generally uses the up-the-ramp coefficients in Equation (3). The following subsections discuss the removal of an artifact in the first read of a CHARIS ramp and our handling of nonlinearity and saturation. We then briefly discuss the read noise properties of our ramps. In Section 5.3, we will take advantage of the relatively low dimensionality of much of the read noise, fitting it out when extracting data cubes.

3.1. Correcting Artifacts in the First Read

The first read in a CHARIS ramp is contaminated by an exponential decay of the reference voltage. Figure 1 shows the difference between the first and second reads of a dark frame, while the top panel of Figure 2 shows the lower-left corner of a 17-read ramp taken in CHARIS’ low-resolution mode without correcting for this decay. Given CHARIS’ baffling and low level of dark current, the mean count rate in a dark frame is much lower than the read noise; Figure 1 shows a noisy decay to zero. For the time series in Figure 1, the two-dimensional readout of each channel has been mapped back to one dimension assuming a pixel rate of 100 kHz and an 8 $\mu$s downtime between reading rows of pixels on the detector. This exponential decay has a time constant of $\sim$21 ms and becomes negligible by 100–200 ms ($\sim$10% of the 1.47 s full-frame readout time). Similar artifacts have previously been noted in H2RGs and in related detectors (Bacon et al., 2004; Rauscher et al., 2007).

The black lines in Figure 1 take the form

[TABLE]

where the time constant $t_{0}=21.1$ ms is the same for all of the lines and the $A_{i}$ are fitted separately to each readout channel $i$ (and to the reference pixels). The lower panel shows the residuals, indicating a good fit with some remaining low-frequency read noise. Removing the first-read artifact requires fitting 34 parameters in all: one time decay constant $t_{0}$, 32 amplitudes for the 32 readout channels, and an additional amplitude for the reference pixels. We perform the fit as follows.

In a CHARIS ramp, the difference between the counts in the first and subsequent reads is the sum of the read noise, the exponential decay of the reference voltage, and the photon rates scaled by the gain. We first estimate the photon count rates by fitting a ramp only to the second and subsequent reads, ignoring the contaminated first read. We then use the fitted pixel-by-pixel count rates and reset values to obtain a modeled number of counts at the first read. The difference between the actual and modeled counts in the first read is the sum of read noise and the reference voltage artifact that we wish to remove. We must then perform a 34-parameter nonlinear fit. Luckily, the fit is only nonlinear in a single parameter, the decay constant $t_{0}$. Once the decay constant is fixed, the 34-parameter nonlinear optimization becomes 33 decoupled one-parameter linear optimizations. By solving these linear problems for each value of $t_{0}$ we may reduce the problem to a one-dimensional nonlinear optimization for $t_{0}$. We first guess the value of $t_{0}$ from the known behavior of the CHARIS detector, and then iteratively fit parabolas to converge to the best-fit decay constant. This is equivalent to Newton-Raphson iteration on ${\rm d}\chi^{2}/{\rm d}t_{0}$. Once we have the decay constant, the other 33 parameters may be obtained by straightforward linear optimizations. The entire process for a typical CHARIS ramp takes a few hundred milliseconds on a laptop computer.

The lower panel of Figure 2 shows the 17-read ramp of the top panel after removing the exponential decay of the reference voltage. There are no longer any artifacts visible, and the first read may now be used to increase the integration time on-source and reduce the read noise. The first read consists of $2048\times 2048$ pixels, many orders of magnitude larger than the 34 parameters describing the reference voltage artifact: fitting out this artifact uses a negligible amount of the information available in the first read.

3.2. Nonlinearity and Saturation

As CHARIS’ H2RG approaches saturation, its response becomes nonlinear. The nonlinearity sets in gradually before the pixel’s response drops sharply to zero. The H2RG also suffers from bleeding into adjacent pixels. After a pixel saturates, its four nearest neighbors see an immediate and compensating increase in their count rates, while its four next-nearest neighbors (the nearest along a diagonal) see a smaller, but still significant, increase. Pixels that do not adjoin the saturated pixel show no significant increase in their count rates until an adjacent pixel saturates.

Figure 3 shows the detector’s response to saturation in a cutout of a 120-read ramp in which the most strongly illuminated pixels saturate in $\sim$30 reads. The lines on the top panel are labeled with the pixel number; an inset in the lower panel provides the key. The increase in count rates in the peripheral pixels closely approximates the lost counts from the saturated pixel(s). The total instantaneous count rate in a $5\times 5$ box decreases by just $\sim$35% between the first few and the last few reads, despite the fact that saturated pixels account for $\sim$70% of the photons in the initial reads.

The H2RG’s saturation behavior means that the saturation of one pixel immediately corrupts the count rates of its eight nearest neighbors, but can be neglected for pixels that are farther away. We fit for each pixel’s count rate using only the reads for which neither the pixel in question, nor any of its immediate neighbors, has saturated.

We measure our H2RG’s nonlinearity using a series of long exposures: a 1000-read ramp taken with the detector almost uniformly illuminated, and sixteen 120–160-read ramps with the detector sparsely illuminated by monochromatic light passing through the lenslets. In both cases, we fit for each pixel’s count rate using only the first $\sim$5% of the reads. We then compare the actual counts at subsequent reads to the expected count rates from the initial reads assuming perfect linearity.

Figure 4 shows the density of points in actual vs. predicted counts on the detector. The density is independently normalized at each count rate ($x$-coordinate) to limit dynamic range; there are several billion points in each figure. A perfectly linear detector should be nearly symmetric about the dashed line $y=x$ (biases from read and photon noise are $<$1% at these count rates). Our H2RG falls below the line, indicating a loss of sensitivity as the pixels approach well capacity. The solid blue line shows our adopted nonlinear response: linear up to 10000 $e^{-}$, continuously matched to a cubic fit up to saturation. Our method is distinct from Finger et al. (2008), who only use unsaturated reads to fit a ramp but assume the pixel response to remain linear.

Different ramps suggest nonlinear responses that differ by $\sim$1–2% near saturation. We make no attempt to explain or account for this in the data, but we do allow the user to add a fixed fraction of the count rate as an uncertainty. This error term also accounts for imperfect modeling of the microspectra; an appropriate value (typically $\sim$5%) gives a $\chi^{2}$ per pixel of order unity after fitting all of the two-dimensional microspectra. This 5% uncertainty is significantly larger than the $\sim$1% uncertainty in the nonlinear response.

Our adopted pixel response is linear up to 10000 $e^{-}$, or $\sim$10% of well capacity. For pixels that remain below this value we therefore use the up-the-ramp fit described in the first paragraphs of Section 3. We individually fit pixels that exceed this threshold. With a fixed nonlinear pixel response, the nonlinear fit has two free parameters: the initial count rate and the reset value. The best-fit count rate $x$ minimizes

[TABLE]

where the sum is over the reads $j$, $b$ is the reset value, and we neglect photon noise. For a linear detector, $f[x\cdot j]=x\cdot j$ with $x$ in units of counts/read and the best-fit $x$ is given by the usual up-the-ramp weights (Equation (3)). We exclude from the fit all pixels with a predicted count rate more than $>$105 $e^{-}$ above reset or with $>$60,000 raw counts ($1.2\times 10^{5}$ $e^{-}$ at a gain of 2; values above 65535 raw counts cannot be represented by unsigned 16-bit integers).

In our case, $f[x\cdot j]$ is a nonlinear function of $x$, and the count rate may no longer be obtained as the weighted sum of the reads with the coefficients from Equation (3). However, at fixed count rate, finding the best-fit reset value $b$ remains a linear problem, and the best $\chi^{2}$ at this count rate is trivial to compute. We use this fact to reduce the two-dimensional nonlinear minimization in $x$ and $b$ to a one-dimensional nonlinear minimization in $x$ (the coefficients of the cubic function define the nonlinear response and remain fixed). We begin with a guess from the up-the-ramp count rate and locally fit a parabola to $\chi^{2}[x]$; the vertex of the parabola is our next guess for the count rate. As with our procedure to fit for the initial exponential decay of the reference voltage, this is equivalent to Newton-Raphson iteration on ${\rm d}\chi^{2}/{\rm d}x$. The algorithm converges in only a few steps and requires a negligible amount of computation for a typical ramp.

3.3. Read Noise Properties

CHARIS’ H2RG detector has $1/f$ read noise that is correlated among the readout channels; this behavior is typical of H2RGs (Moseley et al., 2010). It may be suppressed by measuring by fitting out power on various timescales using interspersed reference pixels (Moseley et al., 2010), using (nearly) unilluminated light-sensitive pixels (Brandt et al., 2013), or using different weightings of the reference pixels at the detector edges (Kubik et al., 2014). In CHARIS, there is an additional component of read noise that is largely shared by alternating readout channels. We read out our detector using 32 readout channels, each $2048\times 64$ pixels in size. The even and odd channels each have their own shared component of read noise: removing a single scaled template from all channels achieves only $\sim$half the noise suppression of using different templates for the even and odd channels.

The power spectrum of CHARIS’ read noise spikes at a range of frequencies in the 1–10 kHz range, most of which do not match any known frequencies of the system. The noise is also highly variable with time, both in amplitude and in its power spectrum. In lab tests during CHARIS’s assembly, the shared read noise was comparable in power to the independent read noise in each channel. On Subaru’s Nasmyth platform, the correlated component of the read noise has gotten much worse: its variance can be more than ten times that of the independent read noise depending on the date and the readout channel. We are currently investigating this source of noise and attempting to fix it in hardware.

Figure 5 shows three representative realizations of CHARIS read noise on three widely separated dates. In all cases, we have used a long ramp to predict the count rate at the second read of the ramp, and then used the difference between the actual and expected counts in the second read as a measure of the read noise. The noise in a CDS image (a difference of two reads) would be $\sqrt{2}$ times the noise shown.

Figure 5 shows that the read noise on Subaru’s Nasmyth platform is a serious problem, with CDS-equivalent read noise $\gtrsim$100 $e^{-}$ in some channels at some times, and that the read noise properties vary with time as CHARIS and neighboring instruments are moved and reconfigured. The lower-right panel of Figure 5 demonstrates how much of this read noise is shared between alternating readout channels. In this panel, we have constructed two read noise templates: one for the even channels, and a separate one for the odd channels. We then rescaled the appropriate template to each channel and removed it. This procedure reduced the read noise by a channel-dependent factor of 1.3 to 5.8. If there were no correlated read noise, this procedure would only be expected to reduce the variance by $\sim$1/16, or the noise by $\sim$3% (since sixteen channels were used to construct each template).

Real CHARIS ramps are packed with microspectra, preventing us from simply removing the correlated read noise as described above. However, the combination of our $\chi^{2}$ method of fitting the microspectra (Section 5.2) and the generous spacing of the spectra on the detector do enable us to achieve effective read noise suppression. We discuss this in more detail in Section 5.3.

4. Calibrations for Cube Extraction using Monochromatic Flatfields

A CHARIS image consists of microspectra arrayed in a grid on the detector (see Figure 2); these must be extracted into a data cube. Each microspectrum is the integral of the monochromatic lenslet spots over the spectrum of light that that lenslet sees, convolved with the pixel response function. Extracting the source spectrum requires knowing the pixel locations where each wavelength of light falls for each lenslet; this is the wavelength solution. Optimal extraction and $\chi^{2}$ extraction, the two main algorithms included in the CHARIS data reduction software, also require knowledge of the monochromatic lenslet spots (the lenslet point-spread function, or lenslet PSF). In the remainder of this paper we will use the terms lenslet PSF and PSFlet interchangeably.

In order to both derive the wavelength solution and measure the wavelength-dependent PSFlets, we inject a supercontinuum source with a narrowband tunable filter into an integrating sphere to uniformly illuminate CHARIS’ lenslet array with monochromatic light. The tunable filter has a width of 5 nm, for a spectral resolution of $\sim$300 at 1.6 $\mu$m, well in excess of CHARIS’ resolving power in its high spectral resolution mode. Such a calibration strategy would need to be revised for an instrument with higher spectral resolution, though a lamp with well-spaced emission lines might serve as an effective substitute for our tunable filter. We have also integrated three narrowband filters, one each in the $J$, $H$, and $K$ bands, into SCExAO’s optics. With one of these filters in place, we may uniformly illuminate CHARIS’s lenslet array with any infrared-bright lamp or even the twilight sky.

We break the CHARIS calibration procedure into two steps. We first use our supercontinuum source and tunable filter to gradually step through wavelength. This enables us to measure the locations of the lenslet PSFs on the detector and to derive the wavelength solution. We also use these images to extract the position-dependent lenslet PSFs. This first step, a full sequence of calibration images, is rarely (if ever) repeated. As of publication, the pipeline includes calibration products derived from a July 2016 calibration sequence. A new calibration sequence offers negligible improvements even one year and several cooling cycles later. The second step is to use a narrowband flat to derive a correction to the wavelength solution once per night. This section describes the process in detail, and summarizes the final calibration files that it produces.

The typical user of CHARIS will not need to derive the full wavelength solution or the monochromatic lenslet PSFs; we have calculated these and distribute them with the source code. The user will, however, need to build the appropriate calibration files for a given observing mode from a single narrowband flat: thermal cycles of the instrument induce small shifts of the microspectra on the detector. We include the script buildcal for this purpose. It takes as input the raw reads from a single narrowband flat and writes all of the necessary calibration files to the directory from which it was run.

4.1. Full Calibration Sequences

A full calibration sequence consists of connecting our supercontinuum source to a tunable filter, coupling it to an integrating sphere, and injecting the light into CHARIS to illuminate the lenslet array. We then step through wavelength with a step size comparable to or slightly larger than CHARIS’s spectral resolution. This yields a series of about 10–15 images per observing mode with well-separated lenslet PSFs. Figure 6 shows these spectrally unresolved flatfields from 1.35 $\mu$m to 2.15 $\mu$m for a small subregion of the detector.

We use these calibration sequences to measure the PSFlet locations for each lenslet as a function of wavelength. We assume a cubic polynomial mapping between integer lenslet coordinates $(i,j)$ and floating point pixel coordinates $(x_{ij}[\lambda],y_{ij}[\lambda])$, with the polynomial coefficients being functions of wavelength:

[TABLE]

We first lightly smooth our monochromatic images with a narrow Gaussian; we use the known lenslet pitch and rotation as our initial guesses for the linear coefficients of the polynomial. We then maximize the sum of the interpolated intensities at the lenslet spot locations by adjusting the coefficients. Once we have derived these coefficients for one wavelength, we estimate the offset in the dispersion direction for the next wavelength using a grid search. We combine this new offset with the polynomial coefficients from the previous step to form the initial guess for this new optimization. We proceed to derive the cubic polynomial transformation from lenslet to detector coordinates for all wavelength steps.

For any wavelength of interest, we can now compute all of the coefficients of the lenslet-detector transformation polynomial (and hence the full wavelength solution) by fitting a cubic polynomial as a function of $\log\lambda$ to each coefficient. We have derived these full wavelength solutions for each CHARIS observing mode using calibration sequences taken in July 2016 and distribute them as part of the CHARIS software package.

Our next step is to reconstruct the lenslet PSFs. For wavelengths $\lesssim$2 $\mu$m, these are undersampled by CHARIS’s H2RG. We take an approach very similar to Anderson & King (2000) for Hubble and to Ingraham et al. (2014) and Draper et al. (2014) for GPI. By deriving the wavelength solution, we already have the location of each lenslet PSF’s centroid and can place it on an oversampled grid. The lenslet PSFs are not spaced by an integer number of pixels nor by a ratio of small integers; as a result, they populate an oversampled PSFlet reasonably densely. We iteratively construct a slightly smoothed PSFlet and fit for each PSFlet’s normalization to refine the smoothed template. Finally, we deconvolve with the smoothing kernel using the Richardson-Lucy algorithm (Lucy, 1974) to produce our final oversampled PSFlets. We measure these oversampled PSFlets at about ten wavelengths in each observing mode.

Figure 7 shows our resulting oversampled PSFlets in 25 subregions of the detector for 1.55 $\mu$m monochromatic light. We use bilinear interpolation to estimate the PSFlet between the centers of regions, and assume the PSFlets to remain constant from the centers of the outer regions to the edges of the detector. The convolution and deconvolution are not completely equivalent because some subpixel offsets are sampled more than others; this could lead to a systematic underestimation of the PSFlet width. Random errors in the PSFlet centers would push in the other direction, leading to a systematic overestimation of PSFlet width. An inspection of the residuals after our $\chi^{2}$ fitting of the microspectra (Section 5.2) indicates that these errors are negligible. Our use of only 25 images to represent PSFlet variation over the detector (motivated by the need to average over a large number of lenslets) is probably a bigger source of error.

As of publication, the pipeline includes one set of calibration files derived from calibration sequences taken in July 2016. We have been unable to conclusively measure any differences in the dispersion, the nonlinear part of the wavelength solution, or in the PSFlet shapes over several calibration sequences and cooling cycles. Wavelength solutions based on more recent calibration sequences have not improved the extracted data cubes so long as a contemporaneous narrowband flat is available (see the following section for details). We will add updated calibration files to the pipeline if and when they prove necessary.

4.2. Narrowband Flats

We use the long calibration sequences described above to measure the lenslet PSFs as a function of both wavelength and position, and to compute the wavelength solution. This wavelength solution changes as CHARIS thermally cycles and as it is craned onto and off of its bench. Small shifts due to motion of the optics or lenslet array change the position and orientation of the wavelength solution, but do not have a measurable impact on its nonlinear component. We use a single narrowband flat for each night to compute these shifts and adjust the wavelength solution accordingly, assuming its nonlinear component to remain fixed.

Our strategy for these narrowband flatfield images is to uniformly illuminate CHARIS’s lenslet array using the halogen flat-field lamp in front of AO188, the system optically upstream of SCExAO. One of three narrowband filters within SCExAO creates a flatfield image that CHARIS cannot resolve spectroscopically. These filters, at 10 nm width, are somewhat broader than our tunable source ($\sim$5 nm); the 1200 nm $J$-band light ($R\sim 120$) is marginally resolved by CHARIS in its $J$-band mode ($R\sim 80$). All filters are completely unresolved in CHARIS’s broadband mode.

We process a narrowband flat using the same procedure as for the full calibration sequence above, except that we use the existing wavelength solution as a starting guess for the spot locations. We perform a grid search to find the approximate offset and then iterate using Powell’s method within the minimization function of scipy.optimize. We keep the offset and linear terms of this new solution, and replace the higher-order terms with their values in the main wavelength solution. These refinements to the wavelength solution are all that we extract from the narrowband flats; they are the only aspects of CHARIS calibration data that are unstable from run to run. Because CHARIS sits on Subaru’s Nasmyth platform and because we designed the imaging relay to be thermally stable (with aluminum mirrors and an aluminum bench), the location of the PSFlets is very stable within a night. Project 1640 and GPI are both mounted at the Cassegrain focuses of their respective telescopes; their PSFlets shift due to flexure as the telescope elevation angle changes (Zimmerman et al., 2011; Wolff et al., 2014).

We include the script buildcal with the CHARIS software for the purpose of building all of the necessary calibration files. The user runs buildcal with the raw file for the narrowband flat as a command-line argument. The software will then compute the calibrations and write all of the files to the directory from which buildcal was run. Some of these files are specific to optimal extraction and $\chi^{2}$ extraction, the two main spectral extraction techniques used by the pipeline. These files are described together with the methods themselves in Sections 5.1 and 5.2. A summary of the wavelength calibration steps and required products for each extraction technique is shown in Figure 8.

4.3. Masking and Flatfielding

Flatfielding an ordinary image divides by a single sensitivity for each pixel; that sensitivity is a combination of illumination by the optics and the quantum efficiency of the pixel. In an IFS like CHARIS, these two components of the flatfield separate into a detector flat, measuring the pixel-to-pixel sensitivity differences, and a lenslet flat, measuring differences in the illumination and transmission of the various lenslets.

We construct the pixel flat from early ramps taken before CHARIS was effectively baffled, when the lenslet array was (relatively) uniformly illuminated. This flatfield does include artifacts from the nonuniform illumination. We therefore apply a high-pass filter, preserving lenslet-to-lenslet variations, but removing the slow variations in illumination across the chip. Individual CHARIS microspectra are $\sim$3$0\times 6$ pixels in size; our high-pass filter is a two-dimensional Gaussian with $\sigma=10$ pixels (full width at half maximum 23.5 pixels) and removes power on scales significantly larger than this.

We use pre-baffling images, together with later images where the background count rate is very low, to identify bad pixels. These are either hot, with a very high dark current, or they respond to light much more weakly than their neighbors. Later images effectively identify the hot and warm pixels, while our method of constructing the pixel flat identifies pixels that are not light sensitive. We flag a pixel as “bad” if it has a dark current 15 sigma above the read noise of its neighbors in a series of long exposure dark frames (i.e. above $\sim$0.4 $e^{-}$ s*-1*), or if its sensitivity in our high-pass-filtered pixel flat is below 80%. We flag 0.6% of pixels as bad in this way; they are clustered in groups large and small across the detector.

Our second step is to construct a lenslet flat. For this step, we return to our full calibration sequences from which we derived a wavelength solution and measured the high-resolution lenslet PSFs. Having extracted these, we reconstruct the monochromatic spot pattern expected for a uniformly sensitive detector and a uniformly illuminated lenslet array. We then scale this ideal spot pattern to the pixel flat described above and determine the best fit lenslet-by-lenslet amplitude of the PSFlets. Because the spot pattern is monochromatic, the spots are very well separated and crosstalk between neighboring lenslets may be safely ignored (see Figure 6).

We thus obtain a lenslet flat for each wavelength in our sequence. The flats are consistent with one another within the observing mode, as is expected for illumination of the field by nearly all reflective optics. This also indicates that the filters, which lie immediately behind the lenslet array, have transmission curves that are nearly spatially uniform. We median-combine the flats to create our lenslet flat for each observing mode. We have also verified that these lenslet flats are consistent with those that we would derive from our narrowband flatfield images, described in the following subsection.

Figure 9 shows the lenslet flat and a subregion of the pixel flat. The pixel flat is the same for all observing modes (it was measured using white light while the detector was poorly baffled), while there is, in principle, a different lenslet flat for each observing mode. In practice, the lenslet flat is nearly the same in all modes due the almost exclusive use of reflective optics throughout CHARIS and the AO systems. The pixel sensitivity generally varies by just a few percent from one pixel to the next. The lenslet flat, apart from a few poorly illuminated (or effectively opaque) lenslets, has typical variations of $\sim$20% across the array.

5. Cube Extraction

Once we have constructed all of the calibration materials, extracting a data cube is relatively straightforward. CHARIS includes two main techniques to perform the extraction: optimal extraction (Horne, 1986), and $\chi^{2}$ extraction. We describe our implementation of each below. The CHARIS pipeline also implements a simple aperture photometry-based extraction. However, this algorithm is unable to account for bad or noisy pixels and has no advantage, either in performance or run time, over optimal extraction.

5.1. Optimal Extraction

Optimal extraction (Horne, 1986) computes the spectral intensity at each wavelength (i.e. each row of pixels perpendicular to the dispersion direction) using both the measured shape of the line-spread function and the pixel-specific measurement errors to weight the pixels. This extraction method has long been the standard approach for spectrographs (Baranne et al., 1996; Cushing et al., 2004; Bolton et al., 2012). We implement optimal extraction for CHARIS assuming Gaussian profiles with a wavelength- and lenslet-dependent width that we measure from our high-resolution PSFlets (Figure 7). This is equivalent to computing the spectral intensity as the normalization of a one-dimensional Gaussian of known position and width and unit area.

We compute the position of each lenslet’s microspectrum as part of the calibration process, but the wavelengths that correspond to integer pixels along the spectrum differ. Each microspectrum is $\sim$30 pixels long in the dispersion direction. We obtain $\sim$30 ($\lambda,\sigma,x,y$) quadruples for each lenslet, where the positions are integers in the dispersion direction $y$ and floating point numbers in the perpendicular direction $x$ (which runs along the center of the microspectrum). The wavelength-dependent width $\sigma$ of the microspectrum is given by computing the second moment of the flux along a line passing through the center of the corresponding PSFlet (Figure 7). Because the PSFlets are not circularly symmetric, the actual profile will include variable contributions from the diffraction spikes at other wavelengths and will differ depending on a lenslet’s spectrum. Optimal extraction takes the ($\lambda,\sigma,x,y$) quadruples for each lenslet (which are computed and saved as part of the calibration step) and uses a weighted sum to calculate a corresponding spectral intensity.

Optimal extraction returns the spectral intensity at the wavelengths corresponding to a given microspectrum’s sampling on the detector. Each microspectrum has its own exact spectral resolution and subpixel sampling, and is therefore defined on its own native wavelength array. It is possible to extract the data without interpolating onto a common wavelength array, but it is difficult to visualize and manipulate such data. CHARIS’ software therefore returns a cube in which these microspectra are interpolated onto a common wavelength array at the spectral resolution requested by the user. The resulting data cube couples neighboring wavelengths because of this interpolation and also because the extracted microspectra are convolved with the line-spread function (a PSFlet’s extent along the dispersion direction of a microspectrum).

The output of optimal extraction is a pair of data cubes: one for the spectral intensities, and one for their errors. The code could easily be modified to return cubes at CHARIS’ native wavelength sampling (which differs for each lenslet). In that case the software would return three cubes: one for the spectral intensities, one for their errors, and one for the native wavelengths.

The CHARIS software also has the ability to perform a naïve aperture extraction using unit weights perpendicular to the dispersion direction. This approach is currently used by SPHERE and GPI (Perrin et al., 2016; Pavlov et al., 2008; Mesa et al., 2015). For CHARIS, aperture extraction produces cubes with slightly more noise (depending on the aperture) and lacks an ability to treat errors and bad pixels: bad or noisy pixels must be replaced with guesses from their neighbors. We do not attempt to replace bad pixels in this case but simply set them to zero. A more careful implementation of aperture photometry would produce a result no better than that from optimal extraction, and would save a negligible amount of computational cost.

5.2. $\chi^{2}$ Extraction

CHARIS’ software implements a $\chi^{2}$-based extraction, fitting every microspectrum with a linear combination of narrowband spots. Once the narrowband spot templates have been computed, it is relatively straightforward to extract the spectrum, i.e., the coefficients of their best-fit linear combination. Computing the templates themselves, however, is more subtle. The spectra on the detector are the spectra seen by the lenslets, and convolved with both the lenslet PSFs and the pixel response function. An attempt to use monochromatic PSFlets to fit the microspectra will suffer from the fact that no real spectrum can be represented as the sum of delta functions. Attempting to overcome this by extracting a spectrum at a resolution significantly higher than CHARIS’ native resolution would cause severe problems with aliasing and covariance between neighboring wavelengths. Our solution is to fit the microspectrum of each lenslet as a series of top-hat spectra in units of $I_{\nu}$, with the spectral resolution of this sampling chosen to be comparable to, or slightly higher than, CHARIS’ intrinsic spectral resolution.

We construct a spot diagram for a narrow spectral range using the oversampled lenslet PSFs like those shown in Figure 7. We use ten monochromatic spots to construct each narrowband spot, first scaling each monochromatic spot by the atmospheric and filter transmission appropriate to the very narrow range of wavelengths it covers. In this way, our narrowband spectra are what we would expect for an astrophysical source with a perfectly flat $I_{\nu}$, average atmospheric transmission, the filter for the given observing mode, and achromatic optics in the rest of the system. The correction is far from perfect, but does mitigate the effect of the atmosphere and filter on our recovered microspectra.

Our calibration routine produces and saves a series of narrow, but not quite monochromatic, arrays of spots stepping through wavelength and built as described above. Together, these spots cover the full wavelength array of a given observing mode with no gaps, and approximately correct for the chromatic throughput of the atmosphere and filter. Figure 10 shows these narrowband spot diagrams. They are broader than those shown in Figure 6, but only slightly.

Finally, we produce one last set of narrowband spots where we retain an oversampling by a factor of 5 in the direction perpendicular to the dispersion. This allows us to use cross-correlation to find the appropriate offset for an individual CHARIS ramp even if a calibration data set is not available from that night (or if it proves a relatively poor match). This takes a significant amount of disk space ($\sim$2 GB per calibration set), but enables the spectra to be located to $\sim$0.01 pixels and more accurately fit by our extraction routine. We do not oversample in the dispersion direction because we currently have no way of getting a sufficiently accurate wavelength solution for an individual image without a narrowband calibration flat.

Once all of the narrowband templates shown in Figure 10 have been computed and saved, $\chi^{2}$ spectral extraction is relatively straightforward. We begin by cutting out a 7 pixel wide, $\sim$35 pixel long rectangle around each microspectrum. We then fit each two-dimensional cutout with 20-25 narrowband spots, minimizing the squared residuals weighted by the pixels’ inverse variance. We use the singular value decomposition for this purpose, implementing it ourselves in Cython to enable lenslet-by-lenslet parallelization using OpenMP. This approach naturally includes the errors on individual pixels and masks hot pixels, and it returns the spectral covariance matrix for each lenslet. It also avoids any interpolation onto a common wavelength array as was necessary for optimal extraction. A similar approach has been implemented for GPI (Draper et al., 2014; Ingraham et al., 2014), but measurements of GPI’s wavelength-dependent lenslet PSFs were insufficient for the algorithm to perform well.

Figure 11 shows the performance of our $\chi^{2}$ extraction on the microspectra in two regions of a low-resolution image: one far from the star where read noise is important (upper panel), and a region in a bright speckle where it is negligible (lower panel). The residuals from systematic errors in the PSFlet models are typically $\sim$5% of the intensity. As we discuss below in Section 5.3, our models (center panels) include the correlated component of the read noise and the undispersed background. This is visible in the top panel, where the left few columns of microspectra are in a noisy readout channel while the rightmost columns are in a much cleaner channel.

As an optional feature, the software can use the PSFlet template file that is oversampled in the direction perpendicular to the dispersion. It uses cross-correlation over $32\times 32$ subregions of the detector to fit a position-dependent subpixel shift in the locations of the spectra. We adopt a range of prospective shifts spaced by 0.2 pixels, compute the cross-correlation in each case, and then fit a parabola to the three cross-correlation values nearest their minimum to obtain the exact subpixel offset. We then use bilinear interpolation to compute an offset for each lenslet from the $32\times 32$ cross-correlation offsets, and interpolate the oversampled PSFlet templates onto the appropriate pixel sampling. This approach typically produces an offset of no more than a few tenths of a pixel, but noticeably improves the residuals in the two-dimensional microspectra. The fit shown in Figure 11 includes this subpixel offset.

A $\chi^{2}$-based extraction has several features that differentiate it from optimal extraction and aperture photometry on the microspectra. Because we fit the entire two-dimensional spectrum, it is trivial to also fit for and remove an undispersed background that is uniform over the microspectrum. For microspectra that straddle two readout channels, we allow for the undispersed background to have different values in the two channels. This fitting of an undispersed background is an optional setting in the software’s configuration.

An important feature of the $\chi^{2}$ extraction is that it automatically attempts to extract the intrinsic source spectrum: it performs a deconvolution with the instrumental line-spread function. This deconvolution results in a negative covariance between neighboring spectral bins and a superficially noisier cube than that produced by optimal extraction. A small amount of smoothing in the spectral dimension, i.e. a reapplication of the convolution with the line-spread function, removes this extra apparent noise. The effect of the line-spread function is apparent in a single slice through a data cube, shown in Figure 12. With optimal extraction, the speckles are radially extended due to contributions from a range of wavelengths. This arises both from the line-spread function and from the fact that we interpolate all of our microspectra (each of which has a different native wavelength sampling) onto a common wavelength array. The speckles in a $\chi^{2}$-extracted cube show a negligible radial extent beyond that of the monochromatic PSF incident on the lenslet array.

A $\chi^{2}$ extraction naturally enables the removal of CHARIS’ spectral crosstalk. In the calibration step, we measured the lenslet PSFs (shown in Figure 7) out to a radius of $\sim$7 pixels, slightly larger than the $\sim$6 pixels horizontally separating the centers of the microspectra. While we only fit the microspectra over a rectangle 7 pixels wide, our model microspectrum extends out to a box roughly twice as wide. Subtracting all of the microspectra results in a small, negative residual as the broader boxes remove a few photons from the neighboring spectra. We iterate one time on the cube to remove crosstalk. In practice, the effects of crosstalk in CHARIS are extremely minor due to the spacing of our microspectra and our use of pinholes on the back of the lenslet array. The correction from an iteration to remove crosstalk is generally $\ll$1%.

5.3. Residual Intensity and Read Noise Suppression

A particularly important feature of our $\chi^{2}$ extraction is that it produces a two-dimensional model of the entire detector readout. Subtracting this model from the actual ramp produces a residual image that, in the limit of perfect models of the microspectra, is pure noise. Residuals from actual CHARIS ramps do show some systematics in the centers of the microspectra but are noise-dominated over much of the detector, both between microspectra and in the less-illuminated lenslets. CHARIS’ H2RG detector has an excessive amount of read noise correlated between the readout channels, as discussed in Section 3.3 and shown in Figure 5. We use our measurement of a residual intensity to model and remove the correlated component of the read noise, achieving a suppression approaching that shown in the lower-right panel of Figure 5.

We fit for the correlated read noise as two patterns, one shared by the even readout channels and a second shared by the odd channels. We also fit for a scalar coupling between each channel and the appropriate noise pattern. We compute the correlated read noise using a trimmed mean over a user-specified fraction, by default 70%, of the pixels with the smallest fitted intensity relative to the read noise, i.e., the pixels where the residual is most dominated by the read noise. Each element of the correlated read noise pattern has sixteen realizations on the detector (half of the 32 readout channels); this approach generally gives at least $\sim$10 pixels from which to calculate the pattern. We then scale the correlated read noise by our fitted couplings and subtract it from each readout channel.

We find that our approach provides excellent suppression of the correlated component of the read noise while avoiding the addition of systematics back into the data. Figure 13 shows an example of a slice through an extracted cube with and without suppressing the read noise and fitting out an undispersed background. The image was taken through a neutral density filter to prevent saturation of the central star, resulting in low signal-to-noise ratios over much of the image. The correlated read noise was severe on the date shown, but it is mostly removed by our algorithm.

We also enable the suppression of correlated read noise when using optimal extraction or aperture photometry. In this case, we use $\chi^{2}$ extraction but save the read noise and (optionally) the undispersed background rather than the data cube. We then subtract the fitted read noise before performing the requested extraction algorithm. This approach has the same effectiveness at removing correlated read noise as simply using $\chi^{2}$ extraction.

Our approach to read noise suppression works because of CHARIS’ redundancy. While the detector has $2048\times 2048$ pixels, we extract only $\sim$22 spectral measurements for each of our $\sim 135\times 135$ illuminated lenslets, using $\sim$10 pixels on average to fit each spectral measurement. The problem is sufficiently overconstrained that we can also fit two $2048\times 64$ noise patterns without facing significant degeneracies.

5.4. Background Subtraction

CHARIS ramps have a background overwhelmingly composed of light leaks and thermal photons in the $K$-band–the true dark current is negligible. The software does have the ability to subtract a background count rate from a two-dimensional ramp (i.e. a matched dark). This is common practice for near-infrared IFSs including OSIRIS (Larkin et al., 2006) and GPI (Perrin et al., 2016). However, it is only worthwhile for CHARIS if a very high signal-to-noise ratio dark frame is available in the same instrument configuration as the science data. Shot noise and read noise from the background frame will add to each image in an observing sequence, and because the background has a single realization of the noise, it will add coherently to all images. The signal-to-noise ratio in the background must be substantially higher than in each science frame (i.e. the integration time must be longer) to avoid this becoming a problem. Matching exposure times will result in the background subtraction contributing $\sim$half the noise to each science frame, and most of the noise to a stack of science frames (because the same realization of noise is added to each frame).

The complexity of subtracting a background image is largely due to the fact that the background contains both a dispersed and an undispersed thermal component. The dispersed component is overwhelmingly due to longer wavelengths, and is orders of magnitude lower with a filter in place that blocks $K$-band light. Microspectra from the thermal background are clearly visible on the detector in the low-resolution and $K$-band modes, but are nearly invisible in the $J$ and $H$ bands or with the neutral density filter in place. Figure 14 shows these backgrounds over part of the detector in the low-resolution, $K$, and $H$ bands. The dispersed component dominates the background in the low-resolution and $K$-bands, while the $H$-band shows a relatively uniform background of $\sim$0.2 $e^{-}\,{\rm s}^{-1}$. The microspectra in $K$-band are dispersed about four times as much as the low-resolution spectra, resulting in \sim$$\frac{1}{4} of the peak intensity.

Any change in the alignment of the lenslet array will change the locations of the thermal microspectra on the array. The CHARIS extraction script can handle shifts in the location of the microspectra by applying a subpixel offset to the template PSFlets, but it cannot use this to shift the two-dimensional array of thermal microspectra. Any mismatch, even if it is just a fraction of a pixel, will degrade the ability of $\chi^{2}$ extraction to fit the background-subtracted microspectra. Unfortunately, we lack the hardware (e.g. a mechanism to block the lenslet pinholes) to measure only the undispersed background.

It is easy to degrade data by subtracting a poor two-dimensional background. This is especially true of low-resolution data taken with the neutral density filter. Exposures with this filter are usually longer integrations to achieve good signal-to-noise ratios, and therefore require long background images to perform a two-dimensional subtraction without adding noise. Subtracting a background taken with the broadband (rather than neutral density) filter results in a dispersed thermal background orders of magnitude too high and a negative inferred sky background in the $K$-band.

In the absence of a matched, high signal-to-noise ratio, two-dimensional background, the software does have the ability to subtract an undispersed background lenslet-by-lenslet as part of a $\chi^{2}$ extraction (Section 5.2), or combined with read noise suppression (Section 5.3). The dispersed background will still be included in the final data cube (though it can easily be removed by standard post-processing algorithms).

5.5. The Extracted Cube

The script extractcube produces a FITS file with the data cube (the intensity at each lenslet and at each wavelength) as HDU 1 and the inverse variance of the data cube (the diagonal of the covariance matrix if using $\chi^{2}$ extraction) as HDU 2. The first HDU, HDU 0, consists only of a header with key information about the observation and the data reduction. The last HDU, HDU 3, contains the original header of the raw reads.

The cube is defined on a logarithmic wavelength array. It is given by the nonstandard header keywords LAM_MIN, LAM_MAX, and DLOGLAM, with

[TABLE]

for $i=0,\ldots,{\tt NLAM}$; the spectral resolution of the extracted cube is $R=1/{\tt DLOGLAM}$. The wavelength scale is also given by the standard FITS keywords (Greisen et al., 2006), with CTYPE3$\,=\,$’AWAV-LOG’. We use a logarithmic scale because of CHARIS’ nearly wavelength-independent dispersion: each pixel in the dispersion direction is a nearly constant increment in $\log\lambda$. A uniform spacing in wavelength would result in aliasing problems at the short wavelength end of the spectrum.

With $\chi^{2}$ extraction, the data cube is never interpolated except over bad spectral measurements (which in any case have their inverse variance set to zero). With aperture photometry or optimal extraction, each lenslet’s microspectrum has been interpolated to place it on a common wavelength array. In any of these cases the inverse variance does not capture the full spectral covariance. Still, it does capture lenslet-to-lenslet variations in the quality of the cube, and prevents bad measurements from corrupting the images in post-processing.

6. Software Parameters and Settings

In order to reduce CHARIS data, we provide a software package that accomplishes the two required tasks: wavelength calibration (through the script buildcal, see Section 4), and cube reduction (the script extractcube, see Section 5). The wavelength calibration only needs to point to a monochromatic flat taken as close as possible to the actual date the science data were taken, ideally the same night. The user can choose whether or not to compute over-sampled PSFlets, which are required to fit for sub-pixel shifts in the ramp as described in Section 5.2.

The cube extraction script, extractcube, has several parameters that the user can set; these will change the photometry and morphology of the reduced spatio-spectral datacube as described in Section 5 and Figure 12. The set of parameters is presented in Table 2 and can be edited within a .ini text file that is read by the extraction software.

Many of the parameters in Table 2 are intended to be left as fixed. Indeed, the only parameter that a typical user will need to modify from its default value is calibdir (which must point to a local directory). Some should never be changed from their default values, while others may be modified to suit a particular reduction or a user’s particular needs. Several of these parameters, however, allow the user to run a faster reduction, to change the treatment of the background, and to write intermediate data products.

Table 3 lists two sets of parameters, one intended for a quick reduction and the other for a more detailed reduction that removes read noise and performs a $\chi^{2}$ extraction. The parameters read_0, read_f, gain, mask, and flatfield are the same in both reductions; we generally recommend that the user never change these values. Changing either flatfield or mask to False, for example, will save a negligible amount of computational effort. In order to extract a cube with bgsub, a thermal background must have been computed by buildcal, while using fitshift requires oversampling the PSFlets in buildcal. The parameter smoothandmask is intended only for cosmetics: it replaces spectral measurements that are much noisier than their neighbors with values taken from a smoothed cube. To avoid biasing the results, any modified intensities have their corresponding inverse variances set to zero. Parameters marked by ‘…’ are ignored, refine because method is not lstsq and minpct because suppressrn is False.

One of the more problematic aspects of the reduction is background subtraction, discussed in Section 5.4. Subtracting a good match to the two-dimensional thermal background (comprising both the microspectra and undispersed background) can improve the reduction. However, it is easy to degrade the final cube by subtracting a noisy background or one that is a poor match to either the thermal microspectra or the undispersed background. For this reason, we generally recommend that a user set bgsub to False. If the user has a background image matching the observing mode and with at least a factor of a few longer integration time than any of the science frames, then setting bgsub to True may improve the reduction. If bgsub is set to False, then the undispersed background may be fit and removed by setting both fitbkgnd and suppressrn to True, or by setting fitbkgnd to True with method equal to lstsq. The user should never attempt to remove the background twice by setting both bgsub and fitbkgnd equal to True.

7. Software Performance and Scaling

The CHARIS software is written in Python and Cython (Dalcin et al., 2011), and uses the extensive libraries available in the NumPy, SciPy (Van Der Walt et al., 2011), and Astropy (Astropy Collaboration et al., 2013) packages. In addition to the parallelization performed natively in these packages, we have used the multiprocessing module to parallelize the construction of model PSFlets and OpenMP to parallelize many of the Cython routines.

We have tested the performance and scaling of the software on a compute server with two 12-core 2.5 GHz Intel Xeon processors. We use maxproc to set the maximum number of threads to allocate and test the two sets of parameters given in Table 3. We use a 41-read ramp (about a 60-second exposure); a shorter exposure would require correspondingly less time to compute the ramp.

Figure 15 shows the results. The software has a few seconds of overhead from the initial module imports and from astropy computations of important parameters of the observation (such as the parallactic angle). The computation of the ramp requires a large amount of I/O and scales with the number of reads; using only the first ten of the 41 reads reduces the computational cost by about a factor of four. For the fast and rough reduction, OpenMP provides only a slight performance gain and the cube extraction is a minor component of the total computational cost. The slower and better reduction, however, sees large performance gains from parallelizing the $\chi^{2}$ extraction step, scaling reasonably well up to $\sim$10 CPUs. This step requires the cube to be extracted three times: an initial extraction, a second extraction once a first estimate of the correlated noise has been removed (Section 5.3), and finally after removing updated estimates of the correlated read noise and spectral crosstalk. Setting refine and suppressrn to False (which we do not recommend) would reduce the computational cost of the $\chi^{2}$ extraction step by nearly a factor of three. Using optimal extraction but setting suppressrn to True would require at least two $\chi^{2}$ extractions (three if also computing an undispersed background). In this case, the total computational time would be almost identical to the results shown here for the slower and better reduction with $\chi^{2}$ extraction.

8. Discussion and Conclusions

In this paper, we have presented the data reduction pipeline for the CHARIS IFS. The software begins with the raw detector reads, corrects for an artifact in the first read, and then reconstructs the count rate at each pixel using either up-the-ramp weighting or the full nonlinear response. This approach handles saturation gracefully as long as there are at least a few reads before a pixel saturates, and it uses only those reads that are uncontaminated by bleeding from a neighboring saturated pixel. Fitting the full nonlinear response to pixels that exceed $\sim$10% of well capacity adds a negligible amount to the required runtime.

There are two steps after constructing the ramp: first, all of the calibration files must be built from a narrowband flatfield image, and second, the data cube must be extracted. Our calibration step uses a range of data products that we have built from early lab data and from detailed image sequences taken with a supercontinuum source and a tunable filter. These products consist of the detector flatfield, the lenslet flatfield, a bad pixel mask, a reference wavelength solution (a wavelength-dependent mapping from lenslet to pixel coordinates), and reconstructed position- and wavelength-dependent lenslet PSFs. We supply a script buildcal with the software; buildcal takes a narrowband filter from the same night as the data to be reduced and computes a perturbation to the wavelength solution and resampled PSFlets for each wavelength and position to be extracted.

Once the calibration step has been completed, the script extractcube extracts the data cube. This script operates on a sequence of reads and requires a set of input parameters that the user may set using a .ini file. It first constructs the count rate at each pixel from the sequence of reads and computes basic parameters of the observation using keywords in the FITS header. Next, buildcal optionally computes a sub-pixel shift in the positions of the microspectra, and finally extracts the data cube.

We implement three cube extraction algorithms: aperture photometry, optimal extraction, and $\chi^{2}$ extraction. The GPI and SPHERE pipelines currently use aperture photometry as their primary extraction methods (Perrin et al., 2016; Pavlov et al., 2008; Mesa et al., 2015). While a least-squares inversion technique has been implemented for GPI (Draper et al., 2014), it has not yet equaled the performance of simple aperture photometry. Project 1640’s pipeline (Zimmerman et al., 2011) comes closest to our approach, but after reconstructing the lenslet PSFs, it uses a weighted sum with the PSFlet templates rather than a $\chi^{2}$ inversion. This decision increases the coupling between neighboring wavelengths, as it effectively performs a second convolution with the line-spread function.

All of the aperture photometry and optimal extraction techniques require interpolation to place the microspectra onto a common wavelength array. This interpolation couples neighboring wavelengths and degrades the data cube, while aperture photometry also requires bad data to be fixed. Our $\chi^{2}$ extraction algorithm avoids all of these drawbacks, and also performs a deconvolution of the microspectra with the line-spread function, resulting in much more monochromatic slices through the data cube. The $\chi^{2}$ algorithm produces a full two-dimensional residual map that allows us to estimate and remove correlated read noise, which results in a dramatic improvement to CHARIS data. This algorithm also allows us to fit out and remove spectral crosstalk and an undispersed lenslet-by-lenslet background, and it produces a full covariance matrix for each lenslet’s microspectrum (though we presently save only the diagonal elements).

The final data product of the CHARIS pipeline is an extracted $(x,y,\lambda)$ data cube with inverse variances for each spatial and spectral measurement. This cube must subsequently be calibrated, using satellite spots induced by SCExAO (Jovanovic et al., 2015a) or some other method. The pipeline currently saves an estimated transformation from lenslet to sky coordinates; we defer a full astrometric calibration to a future paper. We caution against naïvely interpolating our data cube, as doing so would mix spatial and spectral elements that could have very different uncertainties and will inevitably discard information in the cube.

The CHARIS software is written in Python and Cython, and is freely available on github at https://github.com/PrincetonUniversity/charis-dep. It is parallelized, taking anywhere from a couple of seconds to run on a short ramp with a rough reduction to $\sim$100 seconds on a longer ramp with a detailed reduction using a single processor. With a moderately powerful compute server, the most detailed reductions may be completed in $\sim$20–30 seconds per frame. The CHARIS pipeline has a documentation page that we maintain at http://princetonuniversity.github.io/charis-dep/. The data cubes that it produces are suitable inputs for high-contrast image processing, a topic that we will address in forthcoming work.

This work was performed in part under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. This work was performed in part under a Grant-in-Aid for Scientific Research on Innovative Areas from MEXT of the Japanese government (Number 23103002). This research is based on data collected at the Subaru telescope, which is operated by the National Astronomical Observatories of Japan. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.