Post-processing CHARIS integral field spectrograph data with PyKLIP

We present the pyKLIP-CHARIS post-processing pipeline, a Python library that reduces high contrast imaging data for the CHARIS integral field spectrograph used with the SCExAO project on the Subaru Telescope. The pipeline is a part of the pyKLIP package, a Python library dedicated to the reduction of direct imaging data of exoplanets, brown dwarfs, and discs. For PSF subtraction, the pyKLIP-CHARIS post-processing pipeline relies on the core algorithms implemented in pyKLIP but uses image registration and calibrations that are unique to CHARIS. We describe the pipeline procedures, calibration results, and capabilities in processing imaging data acquired via the angular differential imaging and spectral differential imaging observing techniques. We showcase its performance on extracting spectra of injected synthetic point sources as well as compare the extracted spectra from real data sets on HD 33632 and HR 8799 to results in the literature. The pipeline is a python-based complement to the SCExAO project supported, widely used (and currently IDL-based) CHARIS data post-processing pipeline (CHARIS DPP) and provides an additional approach to reducing CHARIS data and extracting calibrated planet spectra.


INTRODUCTION
Direct imaging is a powerful method to discover extrasolar planets and brown dwarfs and characterize their atmospheres and orbits (Currie et al. 2023a). Planet and brown dwarf companions typically lie at small, sub-arcsecond angular separations: imaging them thus requires very high spatial resolution instruments. Thermal emission from these companions is also typically far fainter than that of their host stars -e.g. ∼10 −4 -10 −7 in the near-infrared (IR) for young brown dwarfs and superjovian planets (e.g. Baraffe et al. 2003;Spiegel & Burrows 2012). Companion contrasts in reflected light can be orders of magnitude steeper.
Furthermore, imperfections and deformations in the optics produces quasi-static speckles in the diffraction pattern that cannot be suppressed through increasing integration time (Marois et al. 2000(Marois et al. , 2005Masciadri et al. 2005). To tackle these challenges, the field of exoplanet science has seen exciting advancements in instrumentation such as adaptive optics (AO), coronagraphs, integral field spec-★ E-mail: minghan@physics.ucsb.edu trographs (IFS) etc., as well as in observing techniques and postprocessing algorithms.
The limiting factor in planet detection is the quasi-static speckle noise. Two observing techniques are often used to separate the speckles from any potential astrophysical signal: Angular Differential Imaging (ADI)  and Spectral Differential Imaging (SDI) . ADI is a technique to take sequences of exposures without tracking the rotation of the field of view (FOV), while keeping the telescope optics aligned. An astrophysical signal in such a sequence would move through the FOV relative to the speckle pattern, allowing post-processing algorithms to separate it from the background noise. SDI works on the same principle but along the spectral dimension, where the speckle pattern expands as wavelength increases while an astrophysical signal would remain stationary. Some commonly used algorithms developed over the years include Locally Optimized Combination of Images (LOCI; , Karhunen-Loève Image Projection (KLIP; Soummer et al. 2012;Pueyo et al. 2015), Local Decomposition into Low-rank, Sparse, and Gaussian Noise Components (LLSG; Gomez Gonzalez et al. 2016), etc.
The Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS) is an integral field spectrograph in the near-infrared for the Subaru Telescope (Peters et al. 2012;Groff et al. 2014). CHARIS uses a lenslet-based design to sample the image plane (Peters et al. 2012) and simultaneously obtain a spectrum for each spatial element in a two dimensional (2-D) field of view. This type of design was first implemented on the TIGER IFS (Bacon et al. 1995). This type of full-field spectrograph makes it possible to extract the spectrum of an exoplanet/brown dwarf companion anywhere in the FOV, enabling the characterization of the companion's properties such as temperature, chemistry, and gravity (Barman et al. 2011;Hinkley et al. 2013;Konopacky et al. 2013;Currie et al. 2014Currie et al. , 2018Currie et al. , 2023bMcElwain et al. 2007). CHARIS is also paired with the adaptive optics systems AO188 (Minowa et al. 2010) and the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) instrument (Jovanovic et al. 2015a) to achieve diffraction-limited resolution. CHARIS has a low spectral resolution broadband mode that covers J, H and Ks with 22 evenly spaced wavelength bins in logarithmic scale, as well as high-resolution narrowband modes in J, H, or Ks, with 14 wavelength bins in a single narrowband.
Thus far, all CHARIS science publications have used the cube extraction pipeline from Brandt et al. (2017) -the CHARIS Data Reduction Pipeline (CHARIS DRP) -to convert raw CHARIS data into data cubes consisting of 2D images at different wavelength channels. For processing subsequent to this initial step -e.g. image registration, points spread function (PSF) subtraction and spectrophotometric calibration -nearly all CHARIS studies have used the IDL-based CHARIS Data Processing Pipeline (CHARIS DPP) (Currie et al. 2020b) 1 . The CHARIS DPP pipeline utilizes a variety of algorithms -KLIP and A-LOCI used in combination with ADI, SDI, and reference star differential imaging (RDI) -to model and subtract the stellar PSF, yielding the spectra of known planets, brown dwarfs, and disks and producing new discoveries (e.g. Currie et al. 2018;Goebel et al. 2018;Lawson et al. 2020;Steiger et al. 2021;Lawson et al. 2020Lawson et al. , 2021Kuzuhara et al. 2022;Swimmer et al. 2022;Currie et al. 2022;Swimmer et al. 2022;Currie et al. 2023b). A separate pipeline infrastructure provides an independent check on the detection of companions and their spectral properties (e.g. Macintosh et al. 2015;Keppler et al. 2018) and widens CHARIS's user base.
In this paper, we introduce a CHARIS post-processing pipeline under the framework of the pyKLIP package, a widely used open source python library for direct imaging (Wang et al. 2015). pyKLIP's PSF modeling is primarily based on the KLIP algorithm utilizing principal components analysis (PCA) (Soummer et al. 2012), with an additional PCA-based algorithm, expectation maximization PCA, which we add to the existing framework and present in this paper. This pipeline offers an Python-based option for observers to analyze CHARIS data within a widely-used package framework. The scope of this paper will be limited to processing point sources with ADI and SDI.
We provide an overview of the pyKLIP-CHARIS Post-Processing Pipeline in Section 2. We explain the detailed astrometric calibrations and provide the calibrated plate scales and north pointing in Section 3. Sections 4, 5, 6 and 7 describe the post-processing steps and algorithms implemented in the pipeline. We then demonstrate the final calibrated spectral products of the pipeline and offer comparisons with spectra of well known sources in the literature and with injected synthetic sources in Section 8. Finally, we summarize our work in this paper and discuss potential future updates to the pipeline in Section 9. 1 A Python version of this pipeline should be released sometime in 2023.

PIPELINE OVERVIEW
The pyKLIP-CHARIS Post-Processing Pipeline is part of the pyK-LIP package (Wang et al. 2015), a python library for direct imaging of exoplanets and disks. pyKLIP supports multiple instruments and also has the ability to handle generic data products with appropriate formats. While pyKLIP was able to reduce CHARIS data using the generic data mode, it is cumbersome and time-consuming to manually work out the interfacing and calibrations for CHARIS data. The pipeline we present in this paper fully integrates instrument support into the package with improved image registration and instrument calibrations. We also implement an alternative Weighted Low-rank PSF subtraction method based on the expectation maximization PCA (EMPCA) algorithm (Bailey 2012), available for all supported instruments. We comprehensively test and demonstrate the pipeline's performance and show that it produces consistent results with injected sources as well as known real sources in the literature.
The basic pipeline flowchart is shown in Figure 1. The pyKLIP-CHARIS Post-Processing Pipeline takes as inputs the extracted data cubes produced from the CHARIS raw data reduction pipeline (Brandt et al. 2017). These are wavelength-calibrated 3-D data cubes with two spatial dimensions and one spectral dimension. The first step is to measure the satellite spot positions and triangulate the stellar position for each image taken at each wavelength. This is done automatically after the user passes in a list of extracted data cubes. Then, one typically first runs KLIP or EMPCA reductions over the full frame (still optimizing over each small section) with a few sets of trial parameters to get a detection. If any companion is detected, one can then run KLIP with forward modeling on a single section that contains the companion, which then allows one to fit for the astrometry, photometry as well as extract the spectrum using the KLIP-FM forward modeling technique (Pueyo 2016). A detailed documentation of the pyKLIP-CHARIS Post-Processing Pipeline module and a walk-through of an example CHARIS dataset can be found on the pyKLIP documentation website. 2

Plate Scales and Position Angle Corrections
Distortion correction is the mapping between coordinates in a distorted image to the coordinates in a distortion-free image. Measuring the distortion correction in high contrast imaging data usually only involves calibrating the plate scale (potentially different along and axes) and the position angle (PA) zero point of the instrument. The PA zero point is the direction of true north in the telescope's field of view. For a high contrast imaging instrument, visual binaries with well-known separations can be used as distortion-free references to calibrate the plate scale and north pointing. Higher order distortion effects are usually negligible due to the small fields of view of high contrast imaging instruments. For CHARIS, the lenslets disperse 134 × 134 microspectra onto the detector, giving us a field of view of roughly 2.2 ′′ × 2.2 ′′ .
In this section, we use multiple datasets taken at different times to calibrate the plate scale and north pointing for CHARIS, as well as examine the stability of these parameters over time. We first use a globular cluster to measure the plate scales along the x and y axes separately. This is possible because there are multiple stars in the field of view with different separation vectors, allowing us to constrain the plate scale along x and y axes separately. The drawbacks compared to using a bright visual binary include worse spatial resolution in the reference images (HST images) of the star cluster, and relatively low SNR and Strehl ratio of the stars in CHARIS images, both lead to larger uncertainties in the positions of the reference stars used for calibration. The small field of view of CHARIS can only capture a handful of stars in the cluster (see Figure 2 as an example). The advantage of having multiple stars compared to a binary makes up for the aforementioned drawbacks, such that the constraints on the zeroth and first order distortion coefficients have comparable error bars as those from a binary calibration. However, a few stars is not enough to provide meaningful constraints for higher order distortion corrections with more coefficients. Therefore, we limit our star cluster calibration to only first order corrections. We then use a visual binary to measure an additional calibration result at a different epoch.

Calibration Using the M5 Star Cluster
We express the transformation between lenslet coordinates in a CHARIS image and the corresponding distortion free coordinates as follows: where ( , ) are the lenslet coordinates in an uncalibrated CHARIS image, and are the plate scales along the and axis of a CHARIS image, R is the 2D rotation matrix associated with a counter-clockwise rotation of , ( cal , cal ) are the calibrated coordinates in a distortion free coordinate frame, and ( 0 , 0 ) are translation offsets. We formulate the plate scale along two orthogonal directions separately to test if the plate scales along the x and y axes of the detector are indeed consistent.
We use images of the M5 star cluster taken by CHARIS and HST for our calibration under the formalism of Equation (1). We find the most ideal HST images available for this task to be the ones taken by either the Planetary Camera (PC) of WFPC2, or the Wide Field Channel (WFC) installed on the Advanced Camera for Surveys (ACS), due to their having high resolution, high quality images with few saturated stars. We use the calibrated HST images available on the Hubble Legacy Archive that were taken via these two channels and contain the stars imaged by CHARIS. The images we selected come from HST proposal IDs 6607 (PI Ferraro), 8118 (PI Piotto), and 10775 (PI Sarajedini), listed in Table 1. CHARIS observations have 5 different small fields of view over four epochs, while the HST images cover much larger fields of views. The stars imaged by CHARIS appear in one or two of the three HST images, but not all of them. The pairing of CHARIS and HST images used for calibration is also listed in Table 1. The WFPC2 data was reduced using the OPUS pipeline (Rose 1998) and the ACS/WFC data was reduced using the calibration pipeline CALACS 3 . All HST images we selected have a plate scale of 50 mas/pixel and have the y-axes oriented towards north, after the pipeline calibrations. For ACS/WFC, the relative precision of the pixel scale and the orientation are both on the order of ∼ 6 × 10 −6 (van der Marel et al. 2007), which is negligible compared to CHARIS and won't be factored into our error propagation. For WFPC2/PC, the precision of the the orientation is ∼ 0.03° (Holtzman et al. 1995), which is factored into our error propagation.
To extract the positions of the stars in HST images, we follow the iterative effective PSF (ePSF) building and fitting approach described in (Anderson & King 2000, which is implemented in the Photutils python package (Stetson 1987;Bradley et al. 2021). For CHARIS images, there are not enough stars within the field of view to build ePSFs, and the seeing is also variable over time. Therefore, we fit for a 2D Moffat profile to each star in every halo-subtracted exposure. We match the stars that are present in both CHARIS and HST images, then use the positions of the stars to fit for { 0 , 0 , , , } as given in Equation (1). In this case, and are the scaling factors that convert the CHARIS plate scales along each axis to the known HST plate scale. The angle is the counter-clockwise angle a CHARIS image needs to be rotated to align with an HST image whose y-axis is already aligned with celestial north. Thus, − is equivalent to the PA of the celestial north in a CHARIS image. The parallactic angle, , is recorded for every CHARIS image at the time of exposure. It provides an estimate for the angle : est = + 113°, where 113°is the wavelength dispersion angle of the lenslets on the detector. During post-processing, CHARIS images are rotated counter-clockwise by est to align the y-axis with the celestial north. However, while the dispersion angle is known and very stable, may have some offset relative to the true parallactic angle. This offset may even vary slightly from observation to observation due to small changes in optics alignment. By fitting for the angle , we measure the north pointing correction to the default estimate: This north is the PA zero point in a CHARIS image aligned via the default estimate. The CHARIS field of view is small (2.2 ′′ ×2.2 ′′ ) and can capture only a handful of stars of the M5 star cluster. Between 2017 and 2020, CHARIS took several sequences of the M5 cluster, each sequence capturing a small section of the cluster. We select by eye the CHARIS observations with at least 3 stars within the field of view that are not within 5 lenslet of the boundary. And we require that at least one pair of stars are separated by more than half an arcsecond, because small clumps of closely separated stars would produce large uncertainties 3 https://www.stsci.edu/hst/instrumentation/acs/ software-tools/calibration-tools and biases. We use a total of 61 CHARIS data cubes of M5, observed on 4 different dates with 5 different fields of view that are also imaged by either PC or WFC. A CHARIS data cube consists of images at different wavelengths. We collapse each CHARIS data cube into a single image by taking the mean of the data cube along the wavelength dimension, while discarding the wavelength slices with significant atmospheric absorption (two wavelength channels at 1.37 and 1.87 for CHARIS JHK broadband are discarded). We then use the collapsed CHARIS images for calibration analysis.
We use the emcee package (Foreman-Mackey et al. 2013) to perform MCMC fits of the parameters { 0 , 0 , , , } in equation 1, where the log-likelihood function is given by: where x is the calibrated coordinate vector of the th star, calculated from equation1, and x is the DAOPHOT extracted position of the same star in the HST image.

Calibration Using Visual Binaries
The results of the previous section are consistent with a single plate scale across the field of view, as shown in Section 3.3. Therefore, we also carry out the calibration of plate scale and PA under the assumption of a single consistent plate scale, using images of a visual binary with a well-known separation. In this case, the calibration can be expressed very simply as: where cal is the known separation of the binary in angular units. Δ lens is the measured separation from CHARIS images in lenslets. is the plate scale in mas/lenslet. CHARIS is the PA of the binary companion in CHARIS images using the default north pointing set by the data reduction pipeline. cal is the known PA of the binary. And ℎ is the PA of the true north in CHARIS images. We carried out this type of calibration using two visual binary systems: HIP 55507 and HD 1160. We describe the two datasets in the following two subsections and show our results alongside previous binary calibrations in other works (Kuzuhara et al. 2022;Currie et al. 2018Currie et al. , 2020aCurrie et al. , 2022 in Section 3.3.

HIP 55507
We have a recent CHARIS dataset of HIP 55507, a ≈0. ′′ 75 visual binary. This binary was observed by CHARIS on UT 2022 February 27th as an astrometric calibrator for the characterization of the newly discovered T-Dwarf companion, HIP 21152 B (Franson et al. 2023). HIP 55507 was previously imaged multiple times by Keck/NIRC2, on UT 2012 January 07 and 2015 May 29 (PI Justin Crepp) as part of the TRENDS survey (Gonzales et al. 2020) and also on UT 2021 December 21st (PI Kyle Franson). The calibrated NIRC2 images has an uncertainty in the north position angle of 0°.02 and an uncertainty in astrometry of 0.5 mas (Service et al. 2016). The Keck/NIRC2 astrometry of this companion was precisely measured by Franson et al. (2023) to be 759 ± 5 mas. The non-negligible orbital motion between Keck/NIRC2 and CHARIS observations are accounted for via an orbital fit carried out by the orbit fitting code, Orbits from Radial Velocity, Absolute, and/or Relative Astrometry (orvara), (Brandt et al. 2021). The separation of this binary in CHARIS images is measured by fitting a 2D Airy disk to the companion, HIP 55507 B, in aligned, background subtracted, and exposure-averaged CHARIS images at different wavelengths, and then the results are averaged over wavelength. A full PSF subtraction is not necessary for this calibrator due to the high SNR and diffraction limited PSF of the binary companion. We anchor the CHARIS astrometry to the Keck/NIRC2 astrometry and measure the plate scale and truth north using equations 4 and 5. The results are presented in Section 3.3.

HD 1160
HD 1160 is a system with two wide (sub)stellar companions (Nielsen et al. 2012), one of which (HD 1160 B) is visible within the CHARIS field of view. This system has been used to anchor the astrometric calibration in previous CHARIS papers by comparing CHARIS astrometry of HD 1160 B to that obtained (near-)contemporaneously from the precisely-calibrated Keck/NIRC2 (Kuzuhara et al. 2022;Currie et al. 2018Currie et al. , 2020aCurrie et al. , 2022. Here we add four subsequent calibrations of CHARIS datasets of HD 1160, observed in 2021 October, 2022 January, 2022 September, and 2022 December (PIs: Thayne Currie, Maria Vincent). We pin the astromeric calibration of these data to unsaturated Keck/NIRC2 data obtained in the s and p filters on October 15, 2021 (pinning Oct 2021 CHARIS astrometry), January 15, 2022 (pinning Jan 2022 CHARIS data), and September 8, 2022 (pinning Sept and Dec 2022 CHARIS data) (PIs: Thayne Currie, Taichi Uyama).
The CHARIS data were reduced using the CHARIS DPP; the Keck/NIRC2 data were reduced using a general-use, well-tested code (Currie et al. 2011). In each data set, HD 1160 B is detectable without PSF subtraction techniques, typically at SNR ∼ 10-20 with NIRC2 and SNR ∼ 35-135 with CHARIS. The results are presented in Section 3.3.

Characterizing Precision and Biases
The average centroid error of the itetrative ePSF fitting for stars in the HST images is ∼0.05 pixel, going down to 0.02 pixel for the well isolated bright stars, and up to 0.08 pixel for stars with relatively close neighbors. For CHARIS images, the average centroid error is ∼ 0.7 lenslet. In this section, we further examine the systematic errors and potential biases using mock star cluster data.
For each set of stars within the CHARIS field of view in the real M5 star cluster data, we use the extracted positions of those stars in an HST image as the ground truths for our mock HST observation. To obtain the ground truth positions for the corresponding mock CHARIS observation, we choose a ground truth rotation angle similar to the best fit rotation for the real data, as well as choosing 16 mas/lenslet as the true plate scales for x and y axes, and transform the mock HST true positions into mock CHARIS true positions by reversing equation 1. Then, we generate "measured positions" of these stars by adding random gaussian errors with = 0.1 pixel or lenslet to the ground truth positions for both mock HST and mock CHARIS observations. We then fit for the parameters in equation 1 for a pair of mock observations and compare it to the ground truths to get a residual. To examine the errors and potential biases, we generate 100 pairs of such mock HST and CHARIS observations and look at the average and the scatter of the residuals. Since the precision would depend on the configuration of the calibration stars in the field of view, we perform this investigation separately for each of the 5 fields of view of the real data.

Results
The results for the precision and bias examination of the star cluster calibration method are shown in Figure 3. The data points are the averaged residuals, relative to the ground truth plate scale, over 100 pairs of mock observations for each field of view. The thin black error bars are the RMS scatter of the residuals of the 100 samples in each field of view. The thick blue error bars are the errors on the mean, suggesting that all the residuals are mostly consistent with zero. For every single pair of HST and CHARIS images used in the real data calibration, the RMS scatter (black error bars) are used and propagated to get the final errors of the corresponding field of view for the real data.
The results of the star cluster calibration and binary calibration of the CHARIS plate scale and PA zero point are summarized in Table  2 and Figure 4. The mean values are plotted as the dashed line and the errors on the mean are plotted as the dark shaded regions. The binary calibration results include two measurements (first two in the Binary Calibration section of Table 2) from previous work (Currie et al. , 2020a, with the caveat that these two results are averaged over several epochs. However, our results from individual epochs suggest that there is not significant variation over time.
The M5 star cluster calibration measures the plate scales with errors on the order of ∼0.1 mas. They show that the plate scales along the and axes are consistent with each other to about 1 . The plate scales are also consistent with the plate scales measured from binary calibrations, and seem to be stable over time.
The overall average plate scale of all calibrations as well as all calibrations derived from binary companions is 16.15±0.02 mas/lenslet, which is the same as the 16.15 ± 0.05 mas/lenslet value adopted in previous papers and in the current version of the CHARIS DPP. Measurements from binary calibrations are dominated by uncertainties in the stellar centroid position and the intrinsic SNR of the Keck data to which CHARIS astrometry is pinned. Future unsaturated, high SNR observations of HD 1160 B and other binaries could improve the CHARIS plate scale precision.
Most PA zero point measurements agree within error bars with the default value of −2.2 ± 0.27°, which has been used to analyze CHARIS datasets so far (Currie et al. , 2023b. The average position angle offset from all measurements is −2.03 ± 0.08°. From just the binary calibrators, the average offset is −2.11 ± 0.13°. There are two epochs of M5 calibration data that produced PA zero points that are ∼ 1 above the total average. While 1 is still within reasonable statistical uncertainty, it is worth noting that these two epochs of CHARIS data of the M5 cluster share the exact same field of view, and this correlation led us to suspect there might be some systematic bias. We took a closer look at these two epochs by carrying out a Jackknife resampling of the stars in the field of view used for the calibration. Among the Jackknife samples, the exclusion of the same star in both epochs caused the PA zero point to change more significantly than the exclusion of any other star, which brought both epochs into better agreement with the other measurements. A potential cause is that the star in CHARIS images partially overlaps the diffraction spike created by the support spider of the Subaru Telescope, which could induce some bias in the measured centroid of the star. However, this is not sufficient to conclude that there is real bias induced by the star. Therefore, we present both the results with the exclusion of this star (red data points) and the results without excluding any star (blue data points) in the bottom panel of Figure 4.
Overall, the calibrations suggest that the plate scale and north pointing are stable over time. But it is possible that there could be some small fluctuation of the north pointing over time on the order of the size of the error bars, because the CHARIS + SCExAO instrument is mounted on rubber mounting and can be craned in and out of position between observations (Jovanovic et al. 2015b;Kuzuhara et al. 2022). Therefore, we encourage observers to take contemporary data of a binary calibrator with precisely known astrometry in addition to science data, and use the calibrator to measure the PA zero point for their observation. If contemporary calibration data is not available, then we recommend adopting the average results of all binary calibrations compiled in Table 2 or the average of all binary calibrators. Using this calibration instead of the default one, we note that the revised astrometry for recently-discovered companions shows inconsequential differences with published values: e.g. a difference of ≈0.2 mas in the east and north positions for HIP 99770 b from October 2021 compared to measurement uncertainties of 4 mas (Currie et al. 2023b).

IMAGE REGISTRATION
The first step in processing coronagraphic data is to accurately determine the position of the centroid of the bright primary star in the images. This allows us to re-register all images through interpolation to align them relative to each other in order to model the diffraction pattern of the primary star over a sequence of images.
In high contrast coronagraphic imaging, the bright primary star at the center of the image is occulted by a coronagraph to block out most of the star light. To enable precision astrometry for CHARIS, a diffractive grid generated by sine waves on the deformable mirror (DM) of SCExAO in the pupil plane projects the primary star onto the focal plane with a fixed contrast (Jovanovic et al. 2015a,b;Sahoo et al. 2020), before the star is occulted. These projections of the primary star, which we refer to as "satellite spots", are much dimmer copies of the primary star's PSF that can be imaged within the dynamic range of the camera. The pattern of the diffractive grid can be varied to control the placement of these satellite spots. The period of the grid controls the separation of the satellite spots from the central star in units of / . More periods within the pupil size projects the spots further away from the center. The orientation of the sine waves controls the orientation of the satellite spots, and the sine waves can be placed along one direction or two orthogonal directions on the DM to control the number of spots (2 or 4) (Jovanovic et al. 2015a). Fitting for the positions of these satellite spots allows us to triangulate the centroid of the occulted primary star and then align all the images for PSF/speckle modeling. An early version of the pyKLIP-CHARIS did this with a simple maximum likelihood fit of a 2D PSF model (Airy Function or Gaussian) to every satellite spot, given user specified initial guesses. We refer to this as the "local centroid algorithm" henceforth. This method does reasonably well when initialized close to the correct spot parameters, but is inconsistent at finding the global best-fit and can have large errors for many individual satellite spots. The final centroid of the primary star in each frame is determined by the average position of the spots determined separately at each wavelength; it is thus prone to outliers. We present an alternative global centroid algorithm that fits the centroid of all images in a dataset simultaneously, taking into account the constraint that the spots lie on the  (2018) and Currie et al. (2022). They averaged over three and four epochs of CHARIS data, respectively. b These two PA zero points correspond to the red points in  vertices of a square, as well as the cross-correlation between images, thus improving the stability and robustness of the image centroids.

Global Centroid Algorithm
We adopt a centroiding approach based on both matching the image slices within a data cube at different wavelengths, and matching different data cubes with one another. Both of these components of our method use a variant of the cross-correlation approach commonly used for relative image registration. The scaling of the diffraction pattern with wavelength enables this approach to also constrain the absolute centroid, since diffractive scaling expands the PSF about the star's absolute location on the lenslet array. We adopt an iterative approach that ultimately fits all of the data cubes in an imaging sequence together. We begin by using the central parts of the image to estimate a wavelength-dependent centroid for each cube. We then transition to using the satellite spots to refine the absolute center position, before performing a global fit to the centroid locations in all data cubes. Our first step is to compute an approximate centroid for each image and each wavelength using the inner regions of the images. We adopt the central wavelength channel (typically among the highest signal-to-noise channels) as our reference, and spatially scale all other channels to this reference wavelength. Each image's diffraction pattern will then have the same size in lenslets but, depending on the center used for the spatial scaling, it will have a translation relative to the image at the reference wavelength. The value of the translation may be used to derive the centroid of the images. These are not suitable as final centroids, because the inner regions of the image are strongly affected by the coronagraph and/or saturation (and can be chromatic).
With an initial guess for the centroid as described above, we next introduce a model for the satellite spot locations. We assume them to be at a fixed separation and at four angles from the image centroid (a reference angle plus a multiple of 90 degrees). We further assume that the centroid location itself can be wavelength-dependent, and allow a quadratic shift with wavelength. This could arise, for example, from uncorrected atmospheric dispersion. We refine our estimate of the image centroid at each wavelength, including a possible linear shift with wavelength, use the satellite spots. We optimize for the spot locations that maximize the sum of the interpolated lenslet values at these locations across all wavelength slices. The center of the four spot locations at each wavelength is the image centroid.
Our next step is to repeat the analysis with the satellite spots, but rather than identifying only the location of peak intensity, we use cross-correlation to match the spot locations in different wavelength channels. This is similar to the cross-correlation approach that we used initially. However, rather than the entire central region (whose shape can be strongly chromatic because of coronagraphic effects), we limit the cross-correlation to a region within four lenslets of each satellite spot. This stage of our approach uses the actual PSFs of the satellite spots to achieve better precision without having to rely on an analytic model of these spots. By maximizing the cross-correlation in these limited spatial regions, we derive a new set of centroids at each wavelength, including a possible linear wavelength-dependent drift.
Our final step combines our centroids with a relative image registration approach. For relative registration, we derive the best-fit offset between all pairs of images using cross-correlation. We then perform a least-squares fit to all of the best-fit relative offsets between all pairs of images, obtaining one position per image. This is only defined up to a constant shift, i.e., these centroids are all relative. We therefore constrain the mean position to be the same before and after performing this correction.

Centroid Methods Comparison
To test the performance of the global centroid algorithm and compare it to that of the local centroid algorithm, we examine the stability of the position of a bright companion in a CHARIS dataset under the two centroid algorithms. For this, we use CHARIS observations of the HD 1160 system taken on UT 2018 December 23rd (PI: Thayne Currie). The dataset has 41 exposures taken over roughly half an hour. Each exposure has an integration time of 31 seconds. The parallactic motion over this period is around 20 degrees. The companion, HD 1160 B, has negligible orbital motion during the sequence of exposures. The companion PSFs in these images are diffraction limited, where the first minimum and the maximum of the Airy Disk are visible even without any post-processing. This allows us to locate the companion position without propagating the data through a postprocessing routine, which may introduce additional uncertainties or artifacts that would make it difficult to disentangle the affect of the centroid algorithm alone. We run both centroid algorithms on the dataset. Then we align and derotate the images, and mean-collapse the images in the time dimension. We obtain two combined data cubes this way with high SNRs at all wavelengths, one for each centroid algorithm. We estimate and subtract the background and then fit a 2D Gaussian to the companion at each wavelength of the combine data cubes using MCMC. From this, we obtain best fit positions and widths of the companion PSF at all wavelengths under each centroid algorithm. In principle, the better centroid algorithm would produce companion positions that are more stable over wavelength. The results are shown in Figure. 5. We find that the companion position is much more stable under the global centroid algorithm, while it drifts by over half a lenslet under the local centroid algorithm. This indicates that using global information of all the images to constrain the satellite spots simultaneously is more robust and achieves better alignment. To confirm that the drift present in the local centroid algorithm is in fact due to differences in the measured centroids and not due to unknown errors and biases in the PSF fitting of the companion, we also plot the predicted positions of the companion for the local centroid algorithm by translating the centroid differences found by the two centroid algorithms into the differences in the companion positions. These predictions are shown as the dashed red lines Figure. 5. The agreement between the prediction and the measurements indicates that the drift present in local centroid algorithm is most likely a result of the differences in the two centroid algorithms, rather than unknown errors or biases.

KLIP and Forward Modeling
Once the images are registered and aligned, the next step is to model the diffracted star light that made it past the coronagraph and subtract it from the images. pyKLIP does this using the Karhunen-Loève Image Projection (KLIP) algorithm (Soummer et al. 2012;Pueyo et al. 2015). The KLIP algorithm constructs principal components of the stellar PSF using reference images. For CHARIS, these are typically images of the same star taken at other times and wavelengths to employ angular differential imaging (ADI) and spectral differential imaging (SDI), respectively Sparks & Ford 2002). Images of other stars can also be used to build the principal components for reference differential imaging (RDI), but that ap-proach is not featured in this work. Subtracting the PSF model built using these principal components creates two types of biases: oversubtraction and self-subtraction (Pueyo 2016). Over-subtractions are corruptions that do not depend on the astrophysical signal. This happens when the algorithm fits the companion with speckle noise and over-subtracts flux from the companion. Self-subtractions are corruptions that depend on the real companion signal in question, and can often be interpreted as subtracting part of the companion signal from itself. These biases can be estimated and corrected for using forward modeling and fake source injection techniques. The products produced by subtracting the PSF model from the data will be referred to as "reduced" data/cubes/images henceforth.
KLIP is run on small sectors of the image to optimize the stellar PSF modeling on local patches. These sectors are formed by dividing the image into concentric annuli, and dividing each annulus into several equal-sized sectors. The smaller each optimization sector, the more aggressively the data will be fitted. KLIP also uses a template selection method to model each target optimization section (Marois et al. 2014;Pueyo 2016;Wang et al. 2015). The movement parameter specifies the exclusion criterion for this selection. The ADI observing mode takes exposures with the image rotator turned off. As a result, a potential source moves through the field of view corresponding to the parallactic angle. If a hypothetical source at the target optimization section would have moved less than the number of pixels specified by the movement parameter at another exposure, then that exposure is excluded from the templates in order to mitigate the effect of selfsubtraction. A smaller movement parameter value allows images with smaller field rotations relative to the target image to be included for PSF-fitting, and is more aggressive. Self-subtraction is still present in the PSF-subtracted target section because the companion still remains in all selected templates, only now offset by an amount larger than the movement parameter in the templates compared to the target section. This produces negative wings visible in all KLIP reduced images. This is corrected for with forward modeling. Furthermore, the images are typically high-pass filtered with a Fourier-based implementation in order to remove any smooth large-scale features in the image. This way, the KLIP algorithm is mainly working on modeling individual speckles, each of which is the spatial scale of the diffraction limit.
In addition to the size of the optimization sections and the movement criterion, the number of principal components used to construct the model also influences how aggressively the fit to the stellar PSF is. For pyKLIP, we typically refer to the principal components as "KL modes" and we use "KL " to refer to using the first principal components for the stellar PSF model. Using more KL modes fits for the stellar model/speckle pattern more aggressively. This results in better speckle suppression and, naively, a better detection limit. However, more KL modes also capture more of the signal of any potential companion. This could lead to worse throughput and introduce more over-fitting, which need to be corrected for (Soummer et al. 2012;Pueyo et al. 2015). The SNR of the companion are often maximized at an intermediate level of aggressiveness, which have been mostly empirical choices in previous works (Pueyo 2016;Wang et al. 2018). We use the default values in the pyKLIP package as our starting point to reduce our datasets in this work.
Any point source we detect in the images after running KLIP is distorted due to self-subtraction and over-fitting. To account for this, we use the KLIP Forward Modeling (KLIP-FM) framework (Pueyo 2016) to approximate the over-fitting caused by KLIP. A key input for KLIP-FM is the instrumental PSF, which we can extract from the four satellite spots in the image (see Section 4.1). This enables us to perform accurate PSF fitting to extract information such as the brightness and spectrum of a point source in the data. We will use this in Section 6 to extract the spectrum of a point source in the data.

Weighted Low-rank Approximation
KLIP and related methods based on Principal Component Analysis (PCA) construct a set of basis images. These basis images minimize the squared residual across all pixels and all images used to construct the basis set. The set of basis images can be efficiently computed using the singular value decomposition (SVD). The popularity of PCA is partially due to the simplicity of computing the principal components.
The squared residuals minimized by PCA equally weight all pixels of all images. This makes it impossible to mask a bad pixel or a planet's location by, for example, setting that weight to zero. It also means that pixels with higher counts drive the fit, as the same relative residual translates to a much larger absolute residual. A more flexible approach to overcome these shortcomings is to minimize the weighted squared residuals where the weights , ≥ 0 may be different for each pixel in each image. In Equation (6), the ′ , are linear combinations of the basis images. For example, if the errors on each pixel were known and Gaussian, the weights would be the inverse variances. Unfortunately, the optimization problem defined by Equation (6) may not be solved by SVD, or indeed, by any closed-form approach (Srebro & Jaakkola 2003).
We use the term Weighted Low-rank for this approach. For the general weighted low-rank approximation, iterative algorithms are the only options (Srebro & Jaakkola 2003). Bailey (2012) has implemented one such approach, the expectation-maximization algorithm, or EMPCA. Bailey (2012) applied the algorithm to construct a set of basis spectra for quasars; we turn to the problem of basis images for high-contrast datasets. We adapt the software of Bailey (2012) for use in the pyKLIP-CHARIS Post-Processing Pipeline.
The weighted low-rank approximation makes it straightforward to mask a planet's location in any dataset by setting the weights of the appropriate pixels to zero. Hot pixels and cosmic rays may be handled the same way, without any need for interpolation. Different weights due to photon and background noise are also straightforward to include.
By setting the weights to zero at pixels where a planet appears in any individual image, we may guarantee that none of the planet's flux will contribute to the basis images used for PSF subtraction. In other words, it is possible to completely eliminate self-subtraction. This is not possible with KLIP, as the pixels containing planet light will be in different positions at different times due to parallactic motion. Masking them requires an approach that can handle missing data at any position in any part(s) of a dataset. We can then use these same weights to fit the basis images to the science frames. This step prevents overfitting to enable an unbiased extraction of the planet's flux or spectrum. An example of this is shown in Section 8.4, for an injected synthetic source.
A complete elimination of self-subtraction and oversubtraction requires masking the location of a planet out to several pixels beyond its peak intensity. This reduces the fidelity of the approximation of the stellar PSF, and results in a tradeoff between removing subtraction artifacts and maximizing PSF suppression. The mask must be recomputed for each possible planet location, which also makes it computationally costly to apply to a full dataset.
Weighted Low-rank offers an alternative, complementary data reduction approach to extract spectra with fewer artifacts and smaller required corrections than PCA-based approaches. In fact, PCA is a special case of Weighted Low-rank (with uniform weights), so there is guaranteed to be some region within the vast parameter space of joint pixel weights for which Weighted Low-rank outperforms PCA. Our implementation of Weighted Low-rank in pyKLIP enables users to apply it to any dataset that pyKLIP supports, not just exclusive to CHARIS. Section 8.4 compares the performance of Weighted Low-rank to that of KLIP for a dataset of the planet system around HR 8799. We restrict this comparison to a single case of weights that mask the location of a companion. A full investigation of the possible space of pixel weights is beyond the scope of this paper.

SPECTRAL EXTRACTION
Once a point source is detected in a dataset, the process to extract its spectrum is two-fold, but both use the KLIP-FM features implemented in pyKLIP. First, we measure the position of the source in the wavelength-collapsed broadband image using KLIP-FM and the PSF fitting procedure described in Wang et al. (2016). This generally enables us to measure the position of the point source well below sub-pixel accuracy, even after accounting for uncertainties due to correlated speckle noise in the image. The second step follows the algorithm described in Greenbaum et al. (2018): we fix the position of the point source we want to forward model, and use KLIP-FM again to fit the brightness of the point source in each individual spectral channel. This procedure accounts for how the planet PSF at a given wavelength distorts itself at all other wavelengths due to SDI. Documentation for how to do this full procedure on CHARIS data is available at the pyKLIP docs. 4 To measure the uncertainties on the extracted spectrum, we inject point sources at the same separation as the real companion but at many other position angles and then measure the spread of the recovered spectra for these injected sources. The uncertainties measured this way are those of the PSF subtraction and spectral extraction. These are then propagated with systematic calibration uncertainties from other aspects such as those of the satellite spot contrasts and those of the observed stellar magnitudes. Furthermore, since KLIP-FM corrects for the bias of self-and over-subtraction based on a perturbation-based linear expansion (Pueyo 2016), it may still have significant biases for aggressive reductions for which non-linear effects become more prominent, especially for the brightest targets. Therefore, we also use the injected sources and the recovered spectra to measure the non-linear biases and correct for them.

SPECTROPHOTOMETRIC CALIBRATIONS
After extracting the companion contrast spectrum relative to the satellite spot PSF as described in the previous section, we need to carry out spectrophotometric calibrations to convert the contrast spectrum into physical flux density units.
To do this, we first convert the companion contrasts relative to the satellite spots into contrasts relative to the host star using the satellite spots to host star contrasts: where object, denotes the flux density of the object at wavelength . Note that each individual flux density in this equation is not yet known, only the fractions, i.e. the contrasts, have been calculated. Currie et al. (2020b) used SCExAO's internal source over a narrow bandpass centered on 1.55 m to measure the contrasts of the satellite spots relative to the internal source. They measured a contrast of spot, star, = 2.72 × 10 −3 ± 1.3 × 10 −4 at 1.55 at a 25 nm modulation amplitude of the diffractive grid on the SCExAO deformable mirror. The intensity of the satellite spots, and hence the contrast, should scale as ∝ 2 −2 , where is the modulation amplitude of the diffractive grid, and is the wavelength. Currie et al. (2020b) found excellent agreement with the ∝ −2 scaling in the narrowband of 1.1 − 1.8 (agree to within ∼ 1% or less), as well as excellent agreement with the ∝ 2 scaling. However, they note a difference on the order of ∼ 8% between separate tests which could plausibly be due to slight changes in the placement of SCExAO's focal-plane masks between the times of the two tests. In addition, Wang et al. (2022) also found correlated fluctuations in the satellite spot fluxes over the course of an observing sequence, as well as a 3% difference in the average satellite spot contrast compared to Currie et al. (2020b). As a consequence, high precision spectrophotometric studies are not yet feasible with SCExAO/CHARIS. To minimize the risk of such discrepancies, we recommend the user look into the satellite spot contrast measured on a date closest to the observation date of the dataset the user is analyzing. The satellite spot contrast measured by Currie et al. (2020b) quoted above is recorded as a python class constant in the CHARIS module in the pyKLIP package. It can be trivially assigned to other values if the user finds better suited calibration tests for the satellite spot contrasts. Details of accessing and handling these constants can be found in the CHARIS documentation online.
We then convert the companion's contrast spectrum relative to the host star to a companion spectrum in physical flux density units using information about the host star. Given the star's spectral type, and optionally its temperature and surface gravity, we obtain a stellar model spectrum of the host star. We scale this stellar model to physical flux density units using the observed magnitudes of the star and properly re-sample them at the CHARIS wavelength bins to obtain star, . The companion spectrum in flux density units is then target, = star, target, star, For the stellar model, the Pickles library (Pickles 1998) is widely used. However, for some spectral types, the library of Pickles lacks direct measurements in the near-IR and extrapolates from shorter wavelengths that could lead to miscalibrated companion spectrum . Therefore, we follow the same choice made in Currie et al. 2018 and use the Castelli and Kurucz stellar models (Castelli & Kurucz 2003) as the default option for our pipeline, implemented via the package pysynphot (STScI Development Team 2013). The Pickles library remains as a back-up option for the calibration as well.

FULL REDUCTION DEMONSTRATIONS
In this section, we run the full pipeline reduction on various injected spectra as well as on CHARIS data with known real companions to demonstrate the performance of the pipeline. As mentioned in section 2, CHARIS has broadband (J, H, and Ks) data and narrowband (J, H, or Ks) data with slightly higher resolution. The reduction process is identical and we will only demonstrate the results for broadband data in this paper. We also limit the scope of this paper to point sources.

Contrast Curve
In this section, we show the calibrated contrast curves for two CHARIS observing sequences taken on Sept 1st 2018 and Aug 31st 2020. We also demonstrate the pipeline's performance on known, injected spectra in this section.
We first measure three broadband contrast curves to see how the detection limit depends on the angular separation and integration time for CHARIS. We use two CHARIS datasets to generate these contrast curves: one dataset of HR 8799 (PI Jason Wang) taken on Sept 1st 2018 that covers a parallactic angle rotation of ∼ 165° (Wang et al. 2022); one dataset of HD 33632 (PI Thayne Currie) taken on Aug 31st 2020 that covers a parallactic angle rotation of ∼ 18°. The apparent magnitudes of the host stars and observing conditions are similar for both systems. As a result, the AO corrections are good for both datasets and are of similar quality visually. We also take a subset of the HR 8799 dataset that covers about the same amount of field rotation as the HD 33632 dataset and measure the contrast curve on this subset to examine the consistency of the detection limit for different observation dates with different conditions.
To measure the contrast curve, we use the reduced, wavelengthcollapsed image along with tools implemented in the pyKLIP package. First, pyKLIP assumes azimuthally symmetric noise and compute the 5 noise level at each separation, smoothing out high frequency noise using Gaussian cross correlation and account for small number statistics. This gives us an estimate of the detection limit at each separation, without accounting for the algorithm throughput. Based on this estimate, we inject fake sources a few times brighter than the detection limit, run KLIP reductions and measure the fluxes of these sources to estimate algorithm throughput. We inject fake sources at three different separations: 20 pixels, 40 pixels and 60 pixels. At each separation, we inject the same source at four different PAs. We average the throughput results over PA to obtain throughput estimates at the three different separations. Then, we calibrate the previous contrast curve by correcting for the throughput as a function of separation. The calibrated contrast curves are shown in Figure  6, which provides detection limit estimates as a function of angular separation within the field of view of CHARIS. We see that two datasets with similar observing conditions produce similar contrast curves, while a longer observing sequence with more field rotation improves the detection limit as expected.

Performance on Injected Spectra
We now demonstrate the pipeline's performance on extracting companion spectra using injected sources. For the injected synthetic sources, we use substellar template spectra of three different spectral types-M7, L7 and T1-to represent the typical range of brown dwarfs being observed (Burgasser et al. 2004(Burgasser et al. , 2010. In addition, we also test the pipeline on a flat spectrum to represent sources with few spectral features such as white dwarfs. We inject fake spectra with three different brightnesses: similar brightness to that of the speckles, ∼ 5 times fainter than the speckles, and ∼ 10 times fainter than the speckles. We choose these brightnesses to cover a full range of observations from bright companions to fainter ones just above the detection limit. We show that the pipeline can successfully retrieve spectra in all these scenarios with the appropriate selection of reduction parameters. We inject all spectra at a separation of 45 pixels or ≈725 mas.  Figure 7 shows the images at one wavelength slice of the injected L7 template as an example. From left to right columns are images before reduction, the PSF-subtracted image, and the forward model of the region around the injected template. The effect of oversubtraction is clearly visible as the negative wings around the source; the forward model is shown to account for this effect. From top to bottom panels are the different brightnesses: similar brightness to that of the speckles, ∼ 5 times fainter than the speckles, and ∼ 10 times fainter than the speckles, respectively. Figure 8 shows the four types of extracted spectra: three substellar templates and a flat spectrum, calibrated with both forward modeling and bias correction described in sections 5.1 and 6, respectively. The top to bottom panels are again the three different brightnesses. For the top panel, we use a ADI only reduction and a small KL mode number because the injected source is very bright. The extracted spectra all agree well with the true injected spectra. The middle panel shows the same figure as the top but now for templates that are ∼ 5 times fainter than the speckles. For this reduction, we use a slightly more aggressive reduction by including SDI along with ADI in the PSF modeling. The accuracy of the extracted spectra compared to the injection is visibly a bit worse than the brighter case, but the error estimates also increased correspondingly. There are still good agreements in terms of reduced 2 for this set of extraction. The bottom panel shows the same figure for templates that are ∼ 10 times fainter than the speckles. This is a few times brighter than the measured 5 detection limit. A slightly more aggressive set of reduction parameters is used: 10 KL modes and ADI+SDI. The extracted spectrum still has good agreement with the injection, though the recovered flat spectrum now has more significant deviations from the true spectrum, indicating that a featureless spectrum near the detection limit is difficult to recover. Overall, these tests on injected spectra show that the pipeline can reliably recover the spectra for bright (lower contrast) companions as well as faint (higher contrast) ones that are near the detection limit.

Performance on CHARIS Data
In this section, we use the pipeline to reduce two real CHARIS observations of two different systems with previously characterized companions, HR 8799 and HD 33632, and compare our calibrated spectra with the published results to show that the pipeline successfully re-

Raw Data
Reduced covers spectra that agree with the literature for these well-studied sources.

HD 33632
We first show the pipeline extracted spectrum for a substellar companion to a nearby Sun-like star, HD 33632 Aa. The companion, HD 33632 Ab, was detected in direct imaging by Currie et al. (2020a) using SCExAO/CHARIS. The data cubes were extracted using the CHARIS DRP (Brandt et al. 2017) and then PSF-subtracted and calibrated using the IDL-based CHARIS DPP, which utilized the A-LOCI algorithm in this case (Currie et al. 2012). We use the same CHARIS dataset Currie et al. (2020a) used in their discovery paper, but now we use our pipeline to perform the analysis. The contrast for this companion to its host star in the CHARIS broadband (JHK) is around 5 × 10 −5 in a wavelength collapsed image, about an order of magnitude above the detection limit according to the contrast curve shown in Figure 6. Therefore, a non-aggressive reduction is ideal for this companion in order to obtain a good SNR with as little over-fitting as possible. We choose a few small KL mode numbers and run KLIP using ADI only. We divide every image into 9 annuli and 4 sections for optimization. We then fit for the companion location using the forward model. Finally, we extract and calibrate the spectrum as described in section 6 and 7. In Figure 9, Figure 9. KLIP-FM extracted spectra of HD 33632 Ab at three similar KL modes compared to the calibrated spectra published in Currie et al. 2020a. using the fake source injection method described in section 6. The pyKLIP-CHARIS Post-Processing Pipeline is able to reproduce the published spectrum with good agreement (well within 1 across the entire CHARIS broadband). This agreement is particularly reassuring because the two pipelines use different image registration and PSF-subtraction algorithms, as well as independently implemented spectrophotometric calibrations. This cross-check adds confidence to analyses that depend on the extracted shape of the companion spectrum.

HR 8799
The second dataset we use to demonstrate the pipeline performance is a CHARIS broadband dataset of the HR 8799 system (Wang et al. 2022). HR 8799 is a 1.5 ⊙ "F0" main sequence star (Gray & Kaye 1999;Gray et al. 2003;Baines et al. 2012) with four imaged planets (Marois et al. 2008. This is one of the most famous and well studied directly imaged planetary systems. Three of its four imaged companions, HR 8799 c, d, and e, are within the CHARIS field of view. We extract the spectra for these three planets using the same procedures as for HD 33632 and compare the extracted spectra with the Gemini Planet Imager (GPI) spectra published in (Greenbaum et al. 2018). The CHARIS broadband contrasts of planet c, d and e to the host star are roughly 1.7 × 10 −5 , 1.7 × 10 −5 and 1.2 × 10 −5 respectively. They are more extreme than the contrast for the HD 33632 system. For the case of HR 8799 e, in particular, the contrast is approaching the detection limit of CHARIS due to the small separation (∼25 pixels). Therefore, we use both ADI and SDI to model the PSF, and also increase the KLmodes used for each of these planets accordingly. Figure 10 shows the comparison of our extracted spectra for HR 8799 c, d and e with the GPI spectra. The GPI spectra were also extracted using pyKLIP, but with separate image registration and calibrations unique to GPI (Greenbaum et al. 2018). The GPI spectra cover the H and K bands, with two separate filters K1 and K2 in K band Macintosh et al. (2008). There is a small overlap between K1 and K2 which can be seen in the Figure. For the outer two planets, HR 8799 c and HR 8799 d, our spectra agree with the GPI spectra well with reduced 2 close to 1. For the inner most planet, HR 8799 e, the fluxes in H band again show good agreement between GPI's results and ours. The agreement is worse in K band, particularly for K1 (∼ 1.9 − 2.2 ). We suspect this is an issue with the GPI instrument, as there is also a discrepancy within GPI's own    Greenbaum et al. (2018). Therefore, to compare to an additional reference, we also included GRAVITY's spectrum for HR 8799 e (GRAVITY Collaboration et al. 2019), and we find that CHARIS's spectrum is consistent with GRAVITY's spectrum with a reduced 2 of 0.8. To summarize, we find our spectra for the three HR 8799 planets within the field of view of CHARIS to be consistent with those extracted from GPI's data (with the exception of HR 8799 e in K1 band, but in agreement with GRAVITY's spectrum). This serves as a consistency check for the image registration and calibration steps implemented in the pyKLIP-CHARIS Post-Processing Pipeline, which share the same PSF subtraction and forward modeling implementation as GPI's pyKLIP pipeline.

EMPCA Reduction KLIP Reduction
−5 0 5 10 15 20 Figure 11. PSF-subtracted broadband images (wavelength averaged) of the injected point source that is roughly 5 times fainter than the speckles for EMPCA and KLIP reductions.

EMPCA vs KLIP
In this section we compare the spectra of injected sources extracted using KLIP and EMPCA. EMPCA currently only implemented ADI and does not have SDI or forward modeling capabilities, it is expected to be strictly inferior to KLIP-FM in most situations, as our results will show. However, while implementing forward modeling for EMPCA is a future effort that is beyond the goal of this paper. We show some of its unique advantages when compared with KLIP using ADI-only reductions without forward modeling, and highlight its potential to motivate future effort to develop forward modeling and optimized pixel weights. As mentioned in Section 5.2, EMPCA is an implementation of the Weighted Low-rank Approximation, which is a more general form of PCA that allows individual pixel weights when minimizing the residuals. Therefore, there is guaranteed to be some region within the parameter space of pixel weights that outperforms PCA, which has fixed uniform weights. However, systematically explore and optimize the pixel weights is a task beyond the scope of this paper. We explore several simple choices of pixel weights, , , from Equation 6 and compare them to a KLIP reduction as a demonstration of how the weights can affect the reduction in terms of SNR and spectral fidelity.
We inject a fake source 5 times fainter than the background speckle pattern the same way we did in section 8.2. We run ADI reductions with EMPCA and KLIP, using the same numbers of annuli and subsections. This means the images are divided in to the same optimization sections for both algorithms. But the two algorithms utilize the movement parameters differently. For EMPCA, instead of using the movement parameter for template selection as described in section 5.1, an aperture at the guess location (or known location from a positive detection) of the companion with a radius corresponding to the movement parameter is masked out (i.e. the pixel weights, , , for the aperture are set to zeros), and this mask tracks the parallactic angle such that the same physical position in terms of separation and position angle from the central star is masked out for every exposure. Therefore, the companion will be masked out in all exposures in EMPCA, whereas KLIP only ensures the companion in the templates are not overlapping with the companion in the target exposure. This means that EMPCA could have almost no self-subtraction given a sufficiently large mask over the companion, while KLIP would still have negative wings produced by subtracting the offset companions in the templates.
For EMPCA reductions, we explore a number of different pixel  Figure 12. Comparison of EMPCA and KLIP extracted spectra (ADI only, no forward modeling) of the injected point source that is roughly 5 times fainter than the speckles. The extracted spectra are calibrated using a sample of fake source injections as described in Section 6.
weighting schemes. On a full image scale, we only test two simple choices of weights: 1. Uniform weights , which is equivalent to the weighting in KLIP PCA; 2. Inverse variance weights. On top of this, we explore the addition of a local companion mask, a circular region centered at the given separation and PA for which the pixel weights are set to 0. The extracted spectra have been calibrated for biases using the fake source injection method described in Section 6. Figure 11 shows an example of PSF-subtracted broadband images using EMPCA with a 2-pixel radius mask compared to KLIP with 2-pixel movement. And it shows how the mask in EMPCA is able to effectively eliminate self-subtraction. But after applying corrections for over-subtraction and self-subtraction, measured from the injections of many synthetic sources, the accuracy of EMPCA extracted spectra in terms of reduced 2 is not too different from that of KLIP, though still mostly comparable or slightly better, as shown in Figure  12. Because we are not using SDI or forward modeling in order to compare EMPCA in its current state with KLIP on level grounds, we expect the accuracy of the extracted spectra to be worse than that of KLIP-FM. Spatially varying speckles could introduce bias in this case that are also correlated over wavelength. For our random choice of injection position, it resulted in some excess in H band. However, we only aim to demonstrate that there is room to further optimize the PSF models if we generalize the PCA algorithm to allow free pixel weights using the EMPCA implementation.
A drawback of masking out the companion is that it would decrease speckle suppression within that mask. And if this penalty outweighs the reward of having less self-subtraction, then the SNR of the companion could be worse. We test this on an injected source on the two real datasets described in Section 8.3. The the HD 33632 dataset was taken on UTC 2020 Aug 31st (PI: Thayne Currie), with an exposure time of 31 seconds for each exposure. The whole sequence covers a parallactic field rotation of ∼ 18°. The much longer sequence on HR 8799 was taken on UTC 2018 Sept 1st (PI: Jason Wang), with an exposure time of 20.7 seconds, and covers a field rotation of ∼ 160°. These two datasets have different configurations, observing conditions, and total rotation angles, all of which could influence the performances of the two algorithms. Among these factors, the rotation angle is something we can separate from the other factors and control artificially in this analysis. To do this, for the HR8799 dataset, we modified the parallactic angles of the exposures and injected the point source at positions that simulated a slow and a fast rotating sequence, resulting in two datasets that covered different amounts of field rotation, with one dataset covering ∼ 160°( the original location with un-modified parallactic angles), and the other covering a parallactic rotation of ∼ 18°. As a result, we have a total three datasets that we will reduce using the two algorithms for comparison. To measure the SNR for KLIP, we fit for the peak flux of the companion, then we use an annulus around the central star with a radius same as the companion separation to estimate the noise. For EMPCA, we estimate the noise by masking out many empty regions (no companion) in this annulus and estimating the noise in the reduced data using only these masked regions. This ensures that we are treating the signal region and the noise region the same way so that the SNR is measured from a fair estimate of noise. Figure 13 shows the comparison of the SNRs on a injected source that is 5 times fainter than the speckles for these three datasets. The top two panels show that when the field rotations are small, depending on various other factors, EMPCA could produce similar levels of SNR as KLIP, or it may pay more penalty for the increased noise than the reward of a less distorted signal, at least for the weights that we explored. The bottom two panels show that, when the conditions are the same, for a longer observing sequence with a larger field rotation, the reward from having less distortion in EMPCA outweighs the increased noise in the masked region, and yields a better SNR.
Forward modeling for KLIP accounts for over-subtraction and self-subtraction effects well, and produces better results than EM-PCA, evident by comparing Figure 12 with the middle panel of Figure 8 in section 8. Nevertheless, it is encouraging to see the advantages of EMPCA compared to KLIP without forward modeling, especially because we have only explored the simplest of the pixel weight schemes and there is a lot of room for in-depth analysis on more complex the pixel weights that optimizes the reduction. This also provides motivation to implement both SDI forward modeling for EMPCA as well in the future. Even at its current state, EMPCA can still provide benefits over KLIP (with or without forward modeling) in certain situations. One of which is when there are exposures with bad pixels or small artifacts. KLIP has to discard these exposures as it cannot apply individual weights to pixels, while EMPCA can simply mask out these bad pixels by setting their weights to zero.

SUMMARY AND FUTURE ADDITIONS
We have presented in this work the pyKLIP-CHARIS Post-Processing Pipeline under the framework of the open source python library pyKLIP, which performs PSF modeling and subtraction based on the KLIP algorithm, and extracts the astrometry and spectrum of a companion using forward modeling. We have carried out calibrations for the CHARIS plate scale and PA zero point. We conclude that the plate scale is consistent with the initially designed nominal value and is also stable over time. However, the PA zero point varies significantly from observation to observation, and should be re-calibrated using calibration data taken close to the science data observation time, with minimal instrument set up changes. We have also fully integrated CHARIS compatible spectrophotometric calibrations for extracting planet/brown dwarf spectra. We introduced a new global centroiding algorithm for better image registration and alignment of CHARIS data. We also implemented the EMPCA algorithm that could serve as a potential alternative to KLIP in the event of a significant fraction of bad pixels in the data. We have demonstrated that the pipeline can reproduce the spectra of several real companions published in the literature that are extracted with other pipelines and/or from other instruments' data, as well as recover the spectra of injected   Figure 13. Comparisons of SNRs of an injected point source for EMPCA vs KLIP reductions in three datasets. The top panel shows the results for the source injected into the HD 33632 dataset, which covers a field rotation of ∼ 18°from first to last exposure. The middle panel shows the same comparison for the "slow rotating" HR 8799 dataset that also covers the same field rotation. The bottom panel is for the "fast-rotating" HR 8799 dataset that covers a field rotation of ∼ 160°.
sources. We limit the scope of this paper to point source reductions with ADI and SDI, and we hope this Python-based pipeline within the framework of the widely used package pyKLIP would serve as a complimentary alternative to the IDL-based CHARIS DPP (Currie et al. 2020b) and increases CHARIS's user base to produce more science results. Moving forward, there are several functionalities on the horizon that would be useful to be added to the pyKLIP-CHARIS Post-Processing Pipeline in the future. These include the capability to perform polarimetric differential imaging (PDI) for observations of disks, forward modeling and SDI capabilities specifically derived and implemented for EMPCA, and an in-depth exploration of the vast parameter space of pixel weights for EMPCA to optimize EMPCA reductions.