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 Fig. 2(a) shows the lower-left corner of a 17-read ramp taken in CHARIS’s low-resolution mode without correcting for this decay. Given CHARIS’s baffling and low level of dark current, the mean count rate in a dark frame is much lower than the read noise; Fig. 1 shows a noisy decay to zero. For the time series in Fig. 1, the 2-D readout of each channel has been mapped back to one dimension assuming a pixel rate of 100 kHz and an $8\text{-}\mu \mathrm{s}$ downtime among reading rows of pixels on the detector. This exponential decay has a time constant of $\sim 21\mathrm{ms}$ and becomes negligible by 100 to 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.^41,37

Fig. 1

(a) The exponential decay of the reference voltage in the first read, shown via the difference between the first and second reads in a dark frame. Four channels are shown in green, violet, blue, and orange, with the reference pixels in burgundy. The solid black lines are exponential fits with a shared decay constant of about 21 ms, and (b) shows the residuals. The residuals deviate slightly from zero due to low-frequency noise on the detector. Removing the exponential decay of the reference voltage allows us to use the first read in all CHARIS ramps.

Fig. 2

The lower-left corner of a 17-read ramp taken in CHARIS’s broadband, low-resolution mode (a) before and (b) after removing the exponential decay of the reference voltage in the first read. Fitting the reference voltage requires a negligible fraction of the information contained in the first read. The read may then be used normally to increase the integration time on-source and reduce read noise. The microspectra are visible as short vertical lines with gaps corresponding to the atmospheric absorption bands between the $J$- and $H$-bands and between the $H$- and $K$-bands.

The black lines in Fig. 1 take the form

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 (1-D) 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 $\mathrm{d}{\chi }^2/\mathrm{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.

Figure 2(b) shows the 17-read ramp of (a) 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\text{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’s 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 in Fig. 3(a) are labeled with the pixel number, and an inset in Fig. 3(b) 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.

Fig. 3

Leakage of electrons into neighboring pixels upon saturation. (a) The actual counts (assuming a gain of 2) as the individual pixels shown in the inset saturate. Neighboring pixels absorb the extra electrons upon saturation. (b) In spite of the saturation of pixels 4, 5, 6, 7, and 8 (accounting for nearly 70% of the counts rate in the initial few reads), the total count rate in the $5\times 5\text{pixel}$ box decreases by just $\sim 35\%$ between the first few and the last few 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 have 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 16 120- to 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 with the expected count rates from the initial reads assuming perfect linearity.

Figure 4 shows the density of points in actual versus predicted counts on the detector. The density is independently normalized at each count rate ($x$-coordinate) to limit dynamic range, and 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 $\mathrm{10,000}{e}^{-}$, continuously matched to a cubic fit up to saturation. Our method is distinct from Finger et al.,³⁶ who only use unsaturated reads to fit a ramp but assume the pixel response to remain linear.

Fig. 4

The nonlinear response of CHARIS’s detector as measured (a) in a 1000-read ramp taken before the full baffling was installed and (b) in a series of 120- to 160-read ramps taken with the detector sparsely illuminated by monochromatic light passing through the lenslet array. The color shows the density of measurements normalized to the maximum density at a given extrapolated number of counts ($x$-coordinate); each plot contains data from several ${10}^9$ points. The black-dashed lines show a perfectly linear response, while the blue lines show our fiducial fit to the nonlinear response. The frames differ slightly in the nonlinear response, at the level of $\sim 1\%$ to 2%, as the pixels approach saturation ($\sim \mathrm{100,000}\hspace{0.17em}{e}^{-}$ above reset). Our fit assumes that the pixels have a linear response of up to $\mathrm{10,000}\hspace{0.17em}{e}^{-}$ ($\sim 10\%$ of well capacity).

Different ramps suggest nonlinear responses that differ by $\sim 1\%$ to 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, and an appropriate value (typically $\sim 5\%$) gives a ${\chi }^2$ per pixel of order unity after fitting all of the 2-D microspectra. This 5% uncertainty is significantly larger than the $\sim 1\%$ uncertainty in the nonlinear response.

Our adopted pixel response is linear up to $10000e^{-}$ or $\sim 10\%$ of well capacity. For pixels that remain below this value, we, therefore, use the UTR fit described in the first paragraphs of Sec. 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

In our case, $f[x·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 Eq. (3). However, at a 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 2-D nonlinear minimization in $x$ and $b$ to a 1-D nonlinear minimization in $x$ (the coefficients of the cubic function define the nonlinear response and remain fixed). We begin with a guess from the UTR 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 $\mathrm{d}{\chi }^2/\mathrm{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’s H2RG detector has $1/f$ read noise that is correlated among the readout channels; this behavior is typical of H2RGs.³³ It may be suppressed by fitting out power on various timescales using interspersed reference pixels,³³ using (nearly) unilluminated light-sensitive pixels⁴² or using different weightings of the reference pixels at the detector edges.⁴³ 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\text{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 \text{half}$ the noise suppression of using different templates for the even and odd channels.

The power spectrum of CHARIS’s read noise spikes at a range of frequencies in the 1- to 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 10 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.

Fig. 5

Sample realizations of the read noise on three different dates. The read noise is computed as the difference between the actual number of counts at a single read and the expected number of counts based on the other reads in the ramp. The July 2016, September 2016, and March 2017 images have root-mean-square noise of about $33e^{-}$, $31e^{-}$, and $53e^{-}$, respectively. The strong banding in July 2016 (much fainter but still visible in September) has a modulation at a frequency of 120 Hz. In March of 2017, the read noise shows a dramatic difference in amplitude between the even and odd readout channels. Much of the noise in all of these images is in two patterns, one common to all of the even channels and one to all of the odd channels. The lower-right panel shows the same March 2017 read as the lower-left panel, but with a scaled pattern subtracted from the even readout channels and a different scaled pattern subtracted from the odd channels. The read noise falls by a channel-dependent factor of 1.3 to 5.8 (the variance falls by a factor of $\sim 2$ to 30), when this correlated read noise is removed, and the overall root-mean-square falls from about $53e^{-}$ to about $14e^{-}$. We discuss the removal of correlated read noise from CHARIS ramps in Sec. 5.3.

Figure 5 shows that the read noise on Subaru’s Nasmyth platform is a serious problem, with CDS-equivalent read noise $\gtrsim 100e^{-}$ 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 Fig. 5 demonstrates how much of this read noise is shared among 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 16 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 (Sec. 5.2) and the generous spacing of the spectra on the detector does enable us to achieve effective read noise suppression. We discuss this in more detail in Sec. 5.3.

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 the wavelength with a step size comparable to or slightly larger than CHARIS’s spectral resolution. This yields a series of about 10 to 15 images per observing mode with well-separated lenslet PSFs. Figure 6 shows these spectrally unresolved flatfields from 1.35 to $2.15\mu \mathrm{m}$ for a small subregion of the detector.

Fig. 6

A partial calibration sequence through CHARIS’s broadband mode showing the motion of the lenslet PSFs, or PSFlets, in the dispersion direction (toward the bottom of the page). Only a small region of the detector is shown ($20\times 60\text{pixels}$; a microspectrum is just over 30-pixels long), and the intensity is plotted on a logarithmic scale normalized to the peak of the brightest PSFlet in the region shown. The diagonal spikes are due to diffraction from the square lenslets.

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

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 offset with the polynomial coefficients from the previous step to form the initial guess for this 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 \mathrm{m}$, these are undersampled by CHARIS’s H2RG. We take an approach very similar to Ref. 44 for Hubble and to Refs. 45 and 46 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⁴⁷ to produce our final oversampled PSFlets. We measure these oversampled PSFlets at about 10 wavelengths in each observing mode.

Figure 7 shows our resulting oversampled PSFlets in 25 subregions of the detector for $1.55\text{-}\mu \mathrm{m}$ monochromatic light. We use bilinear interpolation to estimate the PSFlet between the centers of regions and assume that the PSFlets 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 (Sec. 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.

Fig. 7

Oversampled $1.55\text{-}\mu \mathrm{m}$ PSFlets reconstructed over 25 subregions of the detector. The observed PSFlets are resampled by the integer pixel grid of the detector; a $3\times 3\text{pixel}$ grid in the upper-right shows the pixel size. The PSFlets are normalized to unit intensity after resampling. The oversampled PSFlets allow us to construct the pixellated microspectra corresponding to monochromatic and broadband light incident on the lenslet array. The crosses are diffraction spikes from the square lenslets, while the shape variations show the focus changing across the field.

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, in 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 Sec. 4.2 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 that its nonlinear component remains fixed.

Our strategy for these narrowband flatfield images is to uniformly illuminate CHARIS’s lenslet array using the halogen flatfield 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\mathrm{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 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.^48,49

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 Secs. 5.1 and 5.2. A summary of the wavelength calibration steps and required products for each extraction technique is shown in Fig. 8.

Fig. 8

Wavelength calibration steps and required products for the two main reduction methods. A single monochromatic flat for each observing mode is the only data required by the CHARIS user to update all of these products and proceed to cube extraction. These steps are implemented in the script buildcal.

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 that measures the pixel-to-pixel sensitivity differences and a lenslet flat that measures 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 detector was (relatively) uniformly illuminated. This flatfield does include artifacts from the nonuniform illumination. We therefore apply a high-pass filter, preserving pixel-to-pixel variations on the scale of the microspectra but removing the slow variations in illumination across the chip. Individual CHARIS microspectra are $\sim 30\times 6\text{pixels}$ in size, and our high-pass filter is a 2-D Gaussian with $\sigma =10\text{pixels}$ (FWHM 23.5 pixels) and removes power on scales significantly larger than this. We use prebaffling 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.4e^{-}{\mathrm{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 by 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 cross-talk among neighboring lenslets may be safely ignored (see Fig. 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 Sec. 5.

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.

Fig. 9

(a) A subregion of the pixel-based flatfield constructed before effective baffling was installed. We used a high-pass filter to remove the illumination pattern. (b and c) Lenslet flatfield images constructed by comparing the monochromatic spot pattern expected for perfectly uniform throughput and the spot patterns actually observed. The flatfields shown are averaged over wavelength and are very nearly consistent between the low- and high-resolution modes.

5.1.

Optimal Extraction

Optimal extraction⁵⁰ 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.^51–53 We implement optimal extraction for CHARIS assuming Gaussian profiles with a wavelength- and lenslet-dependent width that we measure from our high-resolution PSFlets (Fig. 7). This is equivalent to computing the spectral intensity as the normalization of a 1-D Gaussian of known position, 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\text{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 (Fig. 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’s 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 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’s 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.^22,54,55 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’s 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’s native resolution would cause severe problems with aliasing and covariance among 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’s intrinsic spectral resolution.

We construct a spot diagram for a narrow spectral range using the oversampled lenslet PSFs, such as those shown in Fig. 7. We use 10 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 that 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 it 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 as shown in Fig. 6, but only slightly.

Fig. 10

An illustration of the calibration products for ${\chi }^2$ extraction. The panels labeled $\delta \ln \lambda =0$ use monochromatic spots, while the panels labeled with $\delta \ln \lambda =0.03$ account for the finite bandwidth separating the spectral measurements and include wavelength-dependent atmospheric transmission. The lenslet spots with $\delta \ln \lambda =0.03$ are very slightly (almost imperceptibly) broader than those with $\delta \ln \lambda =0$. When modeling the microspectra of a source with constant $I_{\nu }$, the difference (right panel) between using monochromatic spots (third from right) and spots of finite bandwidth (second from right) is a high-frequency pattern with an amplitude of around 10% of the maximum intensity, peaking at the edges of the bandpass.

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\mathrm{GB}$ per calibration set) but enables the spectra to be located to $\sim 0.01\text{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 Fig. 10 have been computed and saved, ${\chi }^2$ spectral extraction is relatively straightforward. We begin by cutting out a 7-pixel wide, $\sim 35\text{-}\text{pixel}$-long rectangle around each microspectrum. We then fit each 2-D cutout with 20 to 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,^46,45 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 (a) and the other in a bright speckle where it is negligible (b). The residuals from systematic errors in the PSFlet models are typically $\sim 5\%$ of the intensity. As we discuss below in Sec. 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.

Fig. 11

The performance of our ${\chi }^2$ extraction algorithm on low-resolution microspectra in two regions of an image: (a) one far from the star where read noise is significant and (b) one in a bright speckle where it is negligible. We have used a square root stretch to more clearly show residuals. Our model includes the undispersed background and the correlated component of the read noise, computed as described in Sec. 5.3. In the absence of read noise, the residuals are typically $\sim 5\%$ of the input data.

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 2-D microspectra. The fit shown in Fig. 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 2-D 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, as shown in Fig. 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.

Fig. 12

A single $H$-band slice through the same cube extracted three ways: (a) ${\chi }^2$ extraction (Sec. 5.2), (b) optimal extraction (Sec. 5.1), and (c) an extraction using an unweighted 3-pixel-wide aperture. A ${\chi }^2$ extraction intrinsically deconvolves the microspectra and the line-spread function, resulting in a more monochromatic image. Speckles in the optimal extraction slice are radially extended because of both the convolution with the line-spread function and interpolation onto a common wavelength array. The extraction using an unweighted aperture (c) is similar to optimal extraction but with slightly more noise and an inability to handle bad pixels.

A ${\chi }^2$ extraction naturally enables the removal of CHARIS’s spectral cross talk. In the calibration step, we measured the lenslet PSFs (shown in Fig. 7) out to a radius of $\sim 7\text{pixels}$, slightly larger than the $\sim 6\text{pixels}$ horizontally separating the centers of the microspectra. While we only fit the microspectra over a rectangle of 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 cross talk. In practice, the effects of cross talk 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 cross talk 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 2-D 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’s H2RG detector has an excessive amount of read noise correlated among the readout channels, as discussed in Sec. 3.3 and shown in Fig. 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 Fig. 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 16 realizations on the detector (half of the 32 readout channels); this approach generally gives at least $\sim 10\text{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 an ND filter to prevent saturation of the central star, resulting in low SNR over much of the image. The correlated read noise was severe on the date shown, but it is mostly removed by our algorithm.

Fig. 13

A single $H$-band slice through a low-intensity cube extracted with our ${\chi }^2$-based algorithm (a) without and (b) with fitting out the correlated read noise and an undispersed background. The scales are the same on both images. Low-background regions are about 55% as noisy in the lower image as in the upper image. The noisier readout channels are visible as stripes in the slice without read noise suppression.

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’s redundancy. While the detector has $2048\times 2048\text{pixels}$, we extract only $\sim 22$ spectral measurements for each of our $\sim 135\times 135$ illuminated lenslets, using $\sim 10\text{pixels}$ on average to fit each spectral measurement. The problem is sufficiently overconstrained in 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 2-D ramp (i.e., a matched dark). This is common practice for near-infrared IFSs, including OSIRIS⁷ and GPI.²² However, it is only worthwhile for CHARIS if a very high SNR 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 SNR 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 \text{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 ND 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.2e^{-}{\mathrm{s}}^{-1}$. The microspectra in $K$-band are dispersed about four times as much as the low-resolution spectra, resulting in $\sim 1/4$ of the peak intensity.

Fig. 14

(a) Thermal backgrounds in the low-resolution, $K$-, and $H$-band modes. (b) Spectral intensity of the dispersed background in low-resolution and $K$-band modes showing its mean and standard deviation across the lenslet array. The background is composed of both a thermal component incident on the lenslet array and dispersed into microspectra by the prism and a nearly uniform, undispersed component. The dispersed component, which dominates the low-resolution and $K$-band backgrounds, is mostly from $K$-band light and is blocked by the $J$- and $H$-band filters.

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 2-D 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 2-D background. This is especially true of low-resolution data taken with the ND filter. Exposures with this filter are usually longer integrations to achieve good SNRs and, therefore, require long background images to perform a 2-D subtraction without adding noise. Subtracting a background taken with the broadband (rather than ND) 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 SNR and 2-D background, the software does have the ability to subtract an undispersed background lenslet-by-lenslet as part of a ${\chi }^2$ extraction (Sec. 5.2) or combined with read noise suppression (Sec. 5.3). The dispersed background will still be included in the final data cube (though it can easily be removed by standard postprocessing algorithms).

5.5.

Extracted Cube

The script extractcube produces a flexible image transport system (FITS) file with the data cube (the intensity at each lenslet and at each wavelength) as the second Header/Data Unit (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 of only a header with key information about the observation and the data reduction. The last HDU, HDU 3, contains the original header of the file containing 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

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 postprocessing.