The Hyper Suprime-Cam Software Pipeline

TL;DR

The paper presents the development of a comprehensive data processing pipeline for the Hyper Suprime-Cam instrument on the Subaru Telescope, incorporating custom algorithms and large-scale processing for high-quality imaging data reduction.

Contribution

It introduces a customized, scalable pipeline for HSC data processing, integrating novel algorithms and adaptations from LSST's system, optimized for both specific and general observational data.

Findings

Pipeline effectively reduces HSC data to produce science-ready images and catalogs.

Incorporates novel algorithms now part of LSST codebase.

Supports large-scale processing for extensive imaging surveys.

Abstract

In this paper, we describe the optical imaging data processing pipeline developed for the Subaru Telescope's Hyper Suprime-Cam (HSC) instrument. The HSC Pipeline builds on the prototype pipeline being developed by the Large Synoptic Survey Telescope's Data Management system, adding customizations for HSC, large-scale processing capabilities, and novel algorithms that have since been reincorporated into the LSST codebase. While designed primarily to reduce HSC Subaru Strategic Program (SSP) data, it is also the recommended pipeline for reducing general-observer HSC data. The HSC pipeline includes high level processing steps that generate coadded images and science-ready catalogs as well as low-level detrending and image characterizations.

Equations

$${\rho }_1(r)$$

$$z-c\ast N(z)>PSF(d)\left(\bar{z}+cN(\bar{z}))\right)$$

$$y_i\equiv s_i+n_i$$

$${\tilde{s}}_p=\sum _j{D}_{pj}{y}_j;$$

$$L(\alpha ,\mathbf{\mu })$$

$$L(\alpha ,\mathbf{\mu })=L(\hat{\alpha }(\mathbf{\mu }),\mathbf{\mu })+\frac{{\partial }^{2}L}{\partial {\alpha }^2}(\alpha -\hat{\alpha }(\mathbf{\mu }){)}^2$$

$$\hat{\alpha }(\mathbf{\mu })=\frac{\sum _{i}z({\mathbf{r}}_i)\varphi ({\mathbf{r}}_i-\mathbf{\mu })}{\sum _j{\left[\varphi (\mathbf{\mu }-{\mathbf{r}}_j)\right]}^2}=\frac{1}{A}\sum _{i}z({\mathbf{r}}_i)\varphi ({\mathbf{r}}_i-\mathbf{\mu })$$

$$A\equiv \sum _j{\left[\varphi (\mathbf{\mu }-{\mathbf{r}}_j)\right]}^2.$$

$${\sigma }_{\alpha }^2(\mathbf{\mu })={\left[\frac{{\partial }^{2}L}{\partial {\alpha }^2}\right]}^{-1}=\frac{{\sigma }^2}{A}$$

$$\nu (\mathbf{\mu })=\frac{\hat{\alpha }(\mathbf{\mu })}{{\sigma }_{\alpha }(\mathbf{\mu })}=\frac{1}{\sigma \sqrt{A}}\sum _{i}z({\mathbf{r}}_i)\varphi ({\mathbf{r}}_i-\mathbf{\mu }).$$

$$T_i(\mathbf{r})=\min \left(z(\mathbf{r}),z(2{\mathbf{p}}_i-\mathbf{r})\right).$$

$${\alpha }_i⟶\underset{{\alpha }_i}{\min }\left[\sum _{\mathbf{r}}{\left(z(\mathbf{r})-\sum _i{\alpha }_i{T}_i(\mathbf{r})\right)}^2\right].$$

$$C_i(\mathbf{r})=\frac{{\alpha }_i{T}_i(\mathbf{r})}{\sum _j{\alpha }_j{T}_j(\mathbf{r})}z(\mathbf{r}).$$

$$f_{\mathrm{coadd}}^A=\sum _i{w}_i{f}_i^A$$

$$f_{\mathrm{coadd}}^B=\sum _i{w}_i{f}_i^B.$$

$$r_i=\frac{f_i^B}{f_i^A}.$$

$$r_{\mathrm{coadd}}=\frac{f_{\mathrm{coadd}}^B}{f_{\mathrm{coadd}}^A}.$$

$$r_{\mathrm{coadd}}=\frac{f_{\mathrm{coadd}}^B}{\sum _i{w}_i{f}_i^A}=\frac{f_{\mathrm{coadd}}^B}{\sum _i{w}_i\frac{f_i^B}{r_i}}.$$

$$r_{\mathrm{coadd}}=\frac{f_{\mathrm{coadd}}^B}{\sum _i{w}_i\frac{f_{\mathrm{coadd}}^B}{r_i}}={\left(\sum _i\frac{w_i}{r_i}\right)}^{-1}.$$

$$r_{\mathrm{coadd}}\approx \sum _i{w}_i{r}_i.$$

$$\mathbf{M}=\frac{\sum _{\mathbf{r}}\mathbf{r}{\mathbf{r}}^{T}z(\mathbf{r})e^{-\frac{1}{2}{\mathbf{r}}^T{\mathbf{C}}^{-1}\mathbf{r}}}{\sum _{\mathbf{r}}z(\mathbf{r})e^{-\frac{1}{2}{\mathbf{r}}^T{\mathbf{C}}^{-1}\mathbf{r}}},$$

$$\mathbf{M}=$$

$$\approx$$

$$=$$

$$\mathbf{Q}={\left({\mathbf{M}}^{-1}-{\mathbf{C}}^{-1}\right)}^{-1}.$$

$$P(\mathbf{z}|\alpha )={\left(\prod _i\sqrt{2\pi }{\sigma }_i\right)}^{-1}{e}^{-\sum _i\frac{{\left(z_i-\alpha {\varphi }_i\right)}^2}{2{\sigma }_i^2}}.$$

$${\alpha }_{\mathrm{ML}}=\frac{\sum _i{\varphi }_i{z}_i/{\sigma }_i^2}{\sum _i{\varphi }_i^2/{\sigma }_i^2}.$$

$${\alpha }_{\mathrm{MF}}=\frac{\sum _i{\varphi }_i{z}_i}{\sum _i{\varphi }_i^2}.$$

$${\sigma }_{\mathrm{ML}}^2={\left(\sum _i{\varphi }_i^2/{\sigma }_i^2\right)}^{-1}$$

$${\sigma }_{\mathrm{MF}}^2=\frac{\sum _i{\varphi }_i^2{\sigma }_i^2}{{\left(\sum _i{\varphi }_i^2\right)}^2}.$$

$$z_i=\alpha ({\varphi }_i+{ϵ}_i)$$

$${\sigma }_i^2=b+\alpha ({\varphi }_i+{ϵ}_i)$$

Full text

1 Introduction

One of the most important developments in modern astronomy is the steadily increasing survey power of optical telescopes, which has been driven in large part by a generation of wide field-of-view CCD imaging cameras. The science enabled by modern wide imaging surveys with such instruments is extremely broad, because large samples of objects are important for both finding rare objects and making precise measurements of noisy quantities.

Constructing a scientifically useful large catalog from a wide-field photometric survey is a major challenge in its own right, however. Searches for rare objects such as high-redshift quasars or strong gravitational lenses are easily swamped by false positives due to processing problems such as spurious detections or unmasked instrumental artifacts. While follow-up observations of these objects may be necessary for both confirmation and to extract information about e.g. reionization and dark matter substructure (respectively), even producing a candidate sample small enough to vet by-eye (let alone propose for follow-up observations) from a large survey requires careful image processing to avoid these false positives. Similarly, studies of statistical signals such as cosmological weak lensing or galaxy clustering are very sensitive to even small systematic errors in galaxy shape measurements and photometry. Because the tolerance for systematic errors is essentially proportional to the expected statistical error, the quality (and generally the sophistication) of processing must improve as survey area and depth increase to avoid wasting some of the survey’s statistical power.

Most modern surveys have responded to this challenge by building a dedicated pipeline responsible for producing a catalog that can be used for most science without further pixel-level processing, developed by a team of scientists focused on developing new algorithms and understanding the features of the data. The Hyper Suprime-Cam (HSC) team has taken this approach with both its Subaru Strategic Program (SSP) Survey (Aihara et al., 2017a) and general-use observers, which use the same pipeline. We have accomplished this by building the HSC Pipeline as a customization of the Large Synoptic Survey Telescope (LSST Science Collaboration et al. 2009) Data Management (DM) software stack (Jurić et al., 2015), which in turn builds on an algorithmic and conceptual foundation inherited from the Sloan Digital Sky Survey (SDSS) Photo Pipeline (Lupton et al., 2001). A large fraction of the HSC pipeline team is also part of LSST DM, and this relationship is highly symbiotic; HSC provides an invaluable early testbed for algorithms and software being developed for LSST, while LSST provides a much broader team of developers and support scientists than is available purely within the smaller HSC collaboration.

Both the HSC Pipeline and the LSST DM stack are open source software, released under the GNU General Public License (Version 3). This paper does not attempt to provide a tutorial for downloading, installing, and using the HSC Pipeline or the LSST DM stack. Readers should instead consult the LSST DM111http://dm.lsst.org or HSC222http://hsc.mtk.nao.ac.jp/ssp/pipedoc websites for that information.

1.1 HSC Overview

Hyper Suprime-Cam is a 1.77 $\mathrm{deg}^{2}$ imaging camera at the prime focus of the Subaru telescope. It includes 116 2k$\times$4k CCDs (104 science sensors, 4 guide sensors, and 8 focus sensors) with 0.168\arcsec pixels. Some of the most important science it enables comes from the primary associated survey, the HSC-SSP. The HSC-SSP is a multi-layer survey comprised of a 1400 $\mathrm{deg}^{2}$ Wide layer, a 27 $\mathrm{deg}^{2}$ Deep layer, and a 3.5 $\mathrm{deg}^{2}$ UltraDeep layer, with depths of approximately (5$\sigma$ point source, AB magnitude) $r\sim 26$, $r\sim 27$, and $r\sim 28$, respectively. All three layers utilize $grizy$ broad-band filters, with four additional narrow-band filters also used in the Deep and UltraDeep layers. The HSC instrument is more fully described in Miyazaki et al. (in prep.) and Komiyama et al. (in prep.), while the SSP is described in Aihara et al. (2017a). The HSC Pipeline is also a major component of the on-site quality analysis system used by observers (Furusawa et al., in prep.).

Approximately 100 $\mathrm{deg}^{2}$ of the Wide survey were observed to full depth in all bands as of November 2015, and are included in the first public data release (PDR1) produced by the software version (4.0) described in this paper. This data release is described in detail in Aihara et al. (2017b), which also includes additional validation of the HSC Pipeline specific to the SSP dataset. HSC-SSP data releases are served via a web portal and database described in Takata et al. (in prep.).

Because the SSP wide observing strategy is depth-first, the survey already goes significantly deeper than other ongoing ground-based optical surveys such as the Dark Energy Survey (DES; Dark Energy Survey Collaboration et al. 2016) and the Kilo Degree Survey (KiDS; de Jong et al. 2013). This puts HSC in a qualitatively different regime for image processing algorithms, simply because a large fraction of objects are blended (i.e. they have significantly overlapping isophotes at the surface brightness limit of the survey). While astronomical image processing pipelines have always had to deal with blending at some level, at shallow depths blends are rare enough that they can simply be omitted from most analyses. At HSC (and LSST) depths, that sort of cut is simply unacceptable: 58% of objects in the HSC Wide survey are blended, in the sense that they were detected as part of an above-threshold region of sky containing multiple significant peaks in surface brightness. Because most of the blended objects are galaxies, algorithms developed for crowded stellar fields (e.g. Stetson 1987) are not applicable. One advantage HSC has in this regard is its spectacular seeing (0.6\arcsec median FWHM in the $i$-band), but this can also make this deblending challenge more severe; neighbors that might have been hopelessly unresolved in worse seeing may now be possible to deblend, but these are in general the hardest blends to handle because we are limited more by the unknown morphologies of the objects than the (known) point-spread function.

Like DES, KiDS, and LSST, many of HSC’s science goals are cosmological; we seek to measure the distribution and evolution of structure in the universe and use these to infer the properties of dark matter and dark energy. One of the most important cosmological probes in photometric surveys is weak gravitational lensing (e.g. Hoekstra & Jain, 2008), which uses the shapes of distant galaxies to measure the foreground matter distribution via the coherent distortions of those shapes. These distortions are extremely subtle compared with many observational effects, such as convolution with the point-spread function (PSF) and distortion due to the optics, and this puts extremely stringent requirements on the image processing pipeline. At present, all image processing for weak lensing is performed by the main HSC Pipeline; there is no independent or follow-up pipeline for weak lensing. The HSC weak lensing team has put considerable additional effort into validating the HSC Pipeline’s data products, however, as well as an extensive simulation effort to calibrate its results. We refer the reader to Mandelbaum et al. (in prep.) for more information on the performance of the HSC Pipeline for weak lensing.

1.2 Terminology

We adopt the term visit from LSST to denote a single exposure (for LSST, a visit is actually a pair of back-to-back exposures of the same pointing, but these are combined in the earliest stages of processing). Raw and processed CCD-level datasets are thus identified by the combination of a visit ID and CCD ID. Most HSC observations (including all SSP observations) involve multiple dithered exposures at each pointing. We process all visits that overlap each area of sky together, primarily by building coadded images.

The HSC/LSST pipeline divides the sky (or some region of interest on the sky) into tracts. Each tract is a rectangular region with a common map projection, which is further subdivided into patches, which share the tract coordinate system but define smaller regions that are a convenient size for processing and storing images. Tracts overlap on the sky, and patches overlap within a tract; we process these overlap regions multiple times, and then resolve duplicates at the catalog level.

While the pipeline itself allows essentially arbitrary tract definitions, SSP productions use a variant of the rings sky tessellation initially developed by the Pan-STARRS team (Waters et al., 2016), which splits the sky into a number of constant-declination rings (120, in our case) and then divides these in right ascension into approximately square overlapping tracts of constant size. These tracts are approximately 1.68 $\mathrm{deg}$ on a side (slightly larger than the HSC field of view) and they overlap by at least 1\arcmin on each side. Each tract is divided into 9$\times$9 patches that are 4k$\times$4k pixels (0.03434 $\mathrm{deg}^{2}$), with 100-pixel overlaps; these coadd pixel are the same size on the sky (0.168\arcsec) as the raw data pixels.

1.3 Organization of this Paper

Section 2 describes the HSC Pipeline’s software architecture (which is largely inherited directly from LSST DM). Readers interested only in the algorithms from a scientific perspective could consider skipping this section. Section 3 describes the HSC Pipeline at the highest level, and is organized according to the actual sequencing of pipeline stages. Lower-level algorithms that are novel, unusual, or particularly complex are described more fully in Section 4, which also includes demonstrations of their performance. Overall assessment of the pipeline’s performance and future plans are discussed in Section 5.

2 Software Architecture

The HSC Pipeline described here is a fork of the original LSST DM software stack that began to slowly diverge from the main LSST codebase in early 2014. Over the course of 2016, the features introduced on the HSC side were fully reintegrated into the LSST stack, including the ability to fully process raw HSC data, and all future HSC-SSP data releases will be produced with a minimally-modified version of the main LSST stack. This paper thus describes the current algorithmic capabilities of the LSST stack in a general sense, as these have been driven largely by HSC development recently, but it explicitly focuses on the HSC Pipeline implementation and the features (and shortcomings) present in the version used to process the recent SSP PDR1 (Aihara et al., 2017b). All SSP internal data releases thus far have also been produced with the HSC Pipeline fork.

Both the LSST Data Management codebase and the HSC Pipeline are written in a combination of C++ and Python. Flexibility and extensibility are core development goals at all levels of the codebase: we have attempted to make virtually every algorithm not just configurable but replaceable (even by external code). This flexibility is important for LSST in part because many algorithmic choices in LSST processing still require significant research and are under active development. As the HSC Pipeline explicitly demonstrates, it also makes it possible to begin to commission LSST’s pipelines on similar precursor datasets well before LSST commissioning begins.

2.1 Tasks: Configurable Python Algorithms

High-level algorithms are written in Python, primarily as subclasses of our Task class. Each Task subclass defines a set of configuration options as an associated Python Config class, but otherwise Tasks are essentially arbitrary Python callables. Tasks may also be nested, with a parent Task delegating work to one or more “subtasks”, and this nesting is reflected in the associated Config classes; a top-level Config instance holds the configuration for all subtask Configs, forming a complete tree for that unit of processing. Our configuration files are simply Python files that set the attributes of the configuration tree.

Our configuration system also allows any subtask in the hierarchy to be retargeted to a different Task with the same signature without changing any code, simply by adding a Python import statement to a configuration file to load the replacement Task and calling a retarget method on the configuration node representing a subtask. Retargeting is our simplest extensibility mechanism, and it is used by the HSC Pipeline to customize the generic detrending algorithms provided by LSST’s Instrument Signature Removal (ISR) Task with a variant specialized for HSC. All other HSC Pipeline tasks are the same Tasks used by default in the LSST DM stack (though in many cases these Tasks were first developed on the HSC Pipeline fork and then ported back to the LSST codebase).

The most frequently used points of customization are actually in lower-level algorithms that are not written as Tasks. These specialized customization points include algorithms for PSF estimation (see Section 4.3) and source measurement (Section 4.9). As with Tasks, a user-supplied algorithm can be added to the pipeline simply by importing it in a Python-parsed configuration file; this adds the algorithm to the list (called a Registry) of available algorithms for a pipeline step, allowing it to be enabled in the same way internal algorithms are selected.

2.2 Defining the Python/C++ Boundary

Low-level algorithms that perform direct pixel-level processing or other computationally intensive operations are written in C++. We also define most of our low-level primitive classes and data structures (e.g. images, PSF models, geometric regions) in C++, making them available to Python code using SWIG333We have nearly completed a switch from SWIG (http://www.swig.org) to pybind11 (http://pybind11.readthedocs.io); the latest version of the LSST DM codebase uses pybind11 instead.. This is an important philosophical difference from other astronomical software projects that mix Python and C/C++ but define all class interfaces in Python (and pass only simple primitive types such as numbers, strings, and pointers to C/C++). This difference has limited our ability to utilize third-party Python libraries (because we frequently need access to the functionality they provide in C++ as well), and our broader Python/C++ boundary layer represents a sometimes complex category of code that the more Python-centric philosophy largely avoids. But in return, our approach allows us to use the full object-oriented power of C++ when writing some of our most complex and challenging algorithms, making that code more readable and more maintainable without sacrificing performance. It also makes our code more reusable: by providing both Python and C++ APIs, we allow both Python-centric projects and those that heavily utilize C++ to build on our software.

2.3 I/O and Provenance

The HSC Pipeline uses LSST DM’s Data Butler concept to organize and execute all input and output. The butler manages the file format and path of every concrete data product consumed or created by the pipeline (including intermediates), allowing algorithmic code to read or write these data products with just a dataset name (e.g. src, for source catalog) and a data ID that describes the unit of data as a set of key-value pairs (e.g. visit=1228, ccd=40). The butler allows different cameras to define their own key-value pairs, customize file paths, and associate camera-specific metadata with data IDs. These customizations are invisible to the algorithmic code in Tasks, allowing nearly all code to be completely camera-independent. The names of butler datasets provide a convenient and precise way to refer to pipeline data products, and we will use these names in Section 3 when declaring the inputs and outputs of each high-level processing stage. Note that Aihara et al. (2017b) refers to these data products via their file names, and the mapping between these and the butler dataset names is in some cases confusing (for regrettable historical reasons). Table 3 shows the dataset and file names for the HSC Pipeline’s main data products.

The butler organizes the data products created or consumed by individual runs of the pipeline into data repositories. All data products in a repository are produced with the same versions and configuration of the pipeline. Repositories can be chained together – the outputs from one stage of the pipeline can be saved in one repository, and the next stage of the pipeline can use this repository for its inputs while writing its outputs to a new repository. Multiple runs of one pipeline stage with e.g. different configurations but common inputs thus map naturally to a tree of data repositories. All repositories containing pipeline outputs are ultimately chained to a root input repository containing the raw data, calibrations, and metadata information. HSC data repositories are currently implemented as directories with symbolic links to the repository directories they are chained to, but this too is customizable and hidden from algorithmic code; the butler automatically traverses the tree of chained repositories when retrieving data products.

3 Pipelines

The HSC Pipeline is currently comprised of four high-level stages:

• •

in CCD Processing (Section 3.1) we detect sources on each CCD image and use them to calibrate and characterize the image;

• •

in Joint Calibration (Section 3.2) we use repeat observations of sources in different visits to further constrain the astrometric and photometric calibrations;

• •

in Image Coaddition (Section 3.3) we combine visit images into coadds;

• •

and in Coadd Processing (Section 3.4) we detect and measure objects on the coadd images.

These stages are focused on making measurements of the deep, static sky; while the LSST codebase includes some support for image subtraction (to detect transient objects) and visit-level forced photometry (to measure light-curves of variable objects), these are not yet run regularly in HSC-SSP data release processing and are not described further here. Figure 1 shows the relationships between these steps and the datasets they produce and consume.

Throughout all pipeline steps, we approach the problem of missing or bad data from the perspective that we should attempt to process everything and use catalog flags and image masks to indicate results that may be unreliable. Our mask images use each bit of the integer pixel values to represent a different “mask plane”, and for most mask planes we set a corresponding flag field in our catalog for any object that contains a pixel with that bit set. The most common mask planes and flags are listed in Table 3.

3.1 CCD Processing

In CCD Processing, we process each CCD in a visit to produce a calibrated image and catalog. Each CCD is currently processed completely independently; while some algorithms (e.g. background and PSF modeling) can be improved by considering all CCDs in a visit together, this was not implemented in the version used for HSC-SSP releases.

This pipeline consists of a long sequence of semi-iterative steps: many operations have (weak) circular dependencies, so some steps are repeated to take into account improved inputs generated by previous steps. Most of these steps are explained in greater detail in Section 4, so we focus here on the order of operations and the relationships between them.

The steps are as follows.

We run Instrument Signature Removal (ISR; Section 4.1), to assemble raw amplifier-level images into CCD images, perform basic detrending (i.e. flat-field, bias, and dark correction), and apply nonlinearity, brighter-fatter (Section 4.2), and crosstalk corrections. 2. 2.

We detect pixels affected by cosmic rays (Section 4.4) and mask them. This algorithm uses a circular Gaussian with 1\arcsec FWHM width as a stand-in for the PSF model we have not yet estimated, and hence these detections are preliminary. 3. 3.

We repair bad pixels and those affected by saturation by using a linear predictive code to interpolate the values of nearby good pixels (Section 4.5). We defer repairing pixels affected by cosmic rays to step 11 in order to make use of the PSF model. 4. 4.

We estimate and subtract an initial model of the sky background (Section 4.6). 5. 5.

We run a maximum-likelihood detection algorithm (Section 4.7) with a $50\sigma$ threshold and the same $1^{\prime\prime}$ FWHM Gaussian as above as the smoothing filter. Because we are interested only in obtaining enough bright stars to use for PSF modeling and astrometric and photometric calibration, the fact that this filter may be significantly different from the true PSF is not a serious concern. 6. 6.

We re-estimate and subtract the background, ignoring pixels marked as belonging to detections in step 5. 7. 7.

We run centroiding, adaptive moment shape algorithms, and aperture photometry algorithms (Section 4.9) on the detected sources from step 5. We do not deblend prior to measurement at this stage, since our PSF modeling algorithm requires unblended sources. 8. 8.

We use a custom implementation of the “Optimistic B” algorithm of Tabur (2007) to match these to an external reference catalog. We explicitly correct source positions to account for a static approximation to the optical distortion before matching. For HSC-SSP data release processing, we use the Pan-STARRS PV1 catalog (Chambers et al., 2016). We use these matches to fit both a per-CCD magnitude zero-point and a gnomonic world coordinate system transform (TAN WCS) with 3rd-order polynomial distortions. We undo the approximate static distortion correction before fitting so the full distortion is included in the fitted model. 9. 9.

We use the source measurements from step 7 (particularly adaptive moments) to filter out galaxies and most blends, yielding a catalog of secure stars for use in PSF estimation. We reserve 20% of these stars as a validation sample. 10. 10.

We construct a PSF model from images of the other 80% of our secure star catalog, using outlier rejection to remove any barely-resolved galaxies or blends that passed previous cuts. Most of this work is delegated to a restructured version of the public PSFEx code (Bertin, 2013). The details of the star selection and PSF measurement algorithms are described more fully in Section 4.3. 11. 11.

We repeat cosmic ray detection (step 2), and re-run our pixel interpolation algorithm (step 3) to repair pixels affected by cosmic rays. 12. 12.

We re-run source measurement (step 7) on all detected sources, this time utilizing our full suite of measurement algorithms (many of which require or perform better with an accurate PSF model). We also repeat aperture flux and adaptive moments measurements because these depend on centroid measurements that in turn utilize the PSF model. 13. 13.

We estimate aperture corrections (Section 4.9.2) for each photometry algorithm by modeling the spatial variation of the ratio of each algorithm’s flux to the flux we use for photometric calibration, by default a 4\arcsecdiameter circular aperture. The aperture corrections are determined using the same sample of stars used to construct the PSF model, and then applied to the measurements for all sources. 14. 14.

We re-run our astrometric and photometric calibration (step 8), refining the results from that step by using the improved measurements from step 12. 15. 15.

We construct a deeper catalog (repeating step 5), this time using a detection threshold of $5\sigma$ and a Gaussian smoothing filter with the same size as the PSF model. 16. 16.

We once again re-estimate and subtract the background (repeating step 4), using the deeper catalog from step 15 to mask out pixels containing detected objects. 17. 17.

We deblend (Section 4.8) all above-threshold regions containing multiple peaks. 18. 18.

We run our full suite of source measurement algorithms a third time, now on sources in the deeper catalog, using the deblended pixels, and apply the aperture corrections measured in step 13.

CCD Processing generates two primary data products, which we refer to as calexp (for “calibrated exposure”) and src (for “source catalog”). The calexp dataset includes the detrended, background-subtracted image, an integer mask image, and an image containing a per-pixel estimate of the variance. Individual bits in the integer mask pixels correspond to different kinds of features, such as saturation or cosmic rays (see Table 3); the association with is dynamic and is documented in the image headers. The calexp dataset also contains several smaller but generally more complex objects representing the PSF model, the astrometric calibration (WCS), the photometric calibration (magnitude zero-point), and the aperture corrections. Background models are currently stored in a separate dataset (calexpBackground).

In addition to measurement outputs and diagnostic flags (from the final measurement stage on the deeper catalog in step 18), the src catalog also holds objects called Footprints that record the exact above-threshold detection region. These are similar to a per-source version of SExtractor’s “segmentation map” (Bertin & Arnouts, 1996). For deblended sources the Footprints are extended to hold the corresponding array of deblended pixel values.

3.2 Joint Calibration

After individually processing all CCDs, we can refine their astrometric and photometric calibrations by requiring consistent positions and fluxes for sources that appear on different parts of the focal plane in different visits. This allows us to use many more stars than are available in the shallower reference catalog for internal calibration, providing constraints on smaller spatial scales as well as stronger constraints overall.

For both the astrometric and photometric fit, we begin with a spatial match between all of the CCD-level src catalogs that overlap a particular tract in a particular band (each tract$+$band combination is fit independently). We also include external reference catalog sources in this match to tie to absolute astrometric and photometric systems, but we do not require that all matches have a reference catalog source. We select only point sources (via their extendedness; see Section 4.9.10) that were not blended.

Our astrometric model includes a 9th-order polynomial distortion that is continuous over the full focal plane composed with distinct translation and rotation transforms for every CCD. We use the same translation and rotation parameters for all visits in the per-tract fit. While these parameters should be a function of time or telescope configuration rather than sky position, their variation is currently not significant compared to other sources of uncertainty and we have sufficient information to constrain them well in a single tract, so we do not yet need to move away from this simpler approach. We do fit a different polynomial distortion for each visit, as this term models both the optical distortion (which is mostly stable, but may depend subtly on the telescope configuration) and distortions due to the atmosphere (which vary from visit to visit). We use standard polynomials instead of an orthogonal polynomial set (e.g. Chebyshev polynomials) to permit the results to be written to disk using the popular TAN-SIP (Shupe et al., 2005) FITS convention.

For the photometric fit, the model is composed of a per-visit zero-point, a per-visit 7th-order Chebyshev polynomial over the focal plane, and a constant per-CCD scaling. This combines a correction derived from stars to the flat field and a correction for nonphotometric conditions, which are degenerate in this model. We also include a correction for deviations from the nominal on-sky area of each pixel derived from the local Jacobian of the astrometric model. We do not attempt to remove nonphotometric data; because we fit a different Chebyshev polynomial scaling for each visit, this term naturally accounts for both non-uniform illumination of the flat-field screen (which is approximately constant across visits) and gray atmospheric extinction (which is not).

We fit the astrometric and photometric models separately, starting with the astrometric fit. We treat the true fluxes and positions of the stars as additional parameters to be fit, and use a dense multi-threaded LU decomposition with partial pivoting to solve the normal equations at each iteration. This scales poorly when the number of visits in a tract is very large (as in e.g. SSP UltraDeep fields), and an alternate implementation using sparse matrices is currently being developed.

The details of these mathematical models and how we solve them are described more fully in Appendix A.

3.3 Image Coaddition

The traditional approach to processing multi-epoch surveys for static-sky science is to build coadds: images from different observations are resampled onto a common grid and averaged together to construct a single deeper image. Many common approaches to coaddition either degrade the data or introduce systematic errors, however, and recent surveys focused on weak lensing have avoided these problems by instead using models fit simultaneously to images from all epochs for their most sensitive measurements (CFHTLens, Miller et al. 2013; DES, Jarvis et al. 2016; and KiDS, Kuijken et al. 2015). The HSC Pipeline has instead focused on building coadds without introducing systematic errors or discarding information. While it remains to be seen whether this approach will work as the HSC-SSP survey area grows and our tolerance for systematic errors decreases, at the level required for first-year science we believe we have achieved this goal.

The coadds generated by the HSC Pipeline are computed as a direct weighted average of the resampled CCD images – we do not match the PSFs of the input images before combining them. This combination is not optimal (see e.g. Zackay & Ofek 2015), but its signal-to-noise ratio (S/N) and effective PSF are equivalent to those of a single long exposure with the same total exposure time (and observing conditions) if the exposure time is used to set the weights of the input images. We actually use the inverse of the mean variance of the input images to set the weights, so our coadd S/N is actually slightly better than the hypothetical long exposure. This direct coaddition algorithm is significantly simpler to implement than the optimal algorithm, and can easily handle missing data. It also approaches optimality in the limit that the PSFs of all input visits are the same. In the HSC-SSP Wide layer, the mean loss in S/N due to suboptimal coaddition is only $1.7\%$ (with a standard deviation of $1.5\%$).

Prior to combining CCD images, we resample them to a common pixel grid using 3rd-order Lanczos interpolation. We resample each bit plane of the integer mask images using a 2$\times$2 top-hat filter to avoid spreading out masked areas too much. While this means some pixels originally affected by cosmic rays, saturation, or image defects (and repaired via interpolation; steps 3 and 11 in Section 3.1) do contribute to the coadd without being masked on the coadd, this ensures that these contributions to a given coadd pixel are always small compared to contributions to the same pixel from both other epochs and non-interpolated pixels from the same epoch.

The final step of coaddition is to set the BRIGHT_OBJECT mask plane in the coadded images from a catalog of geometric regions whose processing may be affected by bright stars. The construction of these regions is described in Coupon et al. (2017).

Coadds for each patch (see Section 1.2) and filter are built independently. Because patches overlap, and the relative weights of visits are determined separately for each patch, the coadd pixel values in the overlap regions are not identical. These differences are not scientifically significant, but they can slightly complicate the resolution of duplicate detections between patches described in Section 3.4.

3.3.1 PSF Coaddition

Even slightly different PSFs in different visits make it essentially impossible to correctly model the effective PSF of the coadd using the coadded images of stars: even small dithers to fill the gaps between CCDs will yield discontinuities in the effective PSF on the coadd. As the number of dithers increases, the regions of the coadd with a continuous effective PSF become smaller, and it becomes increasingly unlikely that each region will contain even one star usable for PSF modeling.

Instead, the HSC pipeline includes an implementation of the approach of Jee & Tyson (2011), in which the existing PSF models for the input images are resampled and combined with the same coordinate transformations and weights that are applied to the image data. Annis et al. (2014) followed a similar approach, but coadded the PSF at predefined positions on a grid and then fit a smoothly-varying PSF model to these (despite the fact that the PSF being modeled may contain discontinuities). We instead evaluate the PSF only at the position of each object to be measured, and we perform this evaluation “on demand”. Along with the image data for each coadd patch, we store the PSF models, WCS transforms, combination weights, and boundaries of all of the CCD images that went into it. When an image of the effective coadd PSF is requested at a point on the coadd (by, e.g., one of the source measurement algorithms described in Section 4.9), we transform that point to the coordinate system of each of the input CCD images that went into the coadd at that point, evaluate their PSF models at that point, resample them to the coadd coordinate system, and finally combine them (with appropriate weights) to generate an image of the effective coadd PSF. We cache the results so repeated evaluations of the PSF model at a single point will not require additional computation.

We typically assume that the model PSF at the center of an object is appropriate for the entirety of the object. Usually the spatial variation is sufficiently slow that this approximation is perfectly valid even for large objects, but it is invalid where the number of inputs to the image change discontinuously. This can happen both at the edges of input images and in regions where one or more pixels were not included in the coadd due to masking or clipping (see Section 3.3.2). Because there is no well-defined PSF for the entire object in these cases, we simply attempt to flag such objects as having been measured with such PSFs by setting the flags_pixels_clipped_any flag. In most cases, these PSF inaccuracies are quite small, and science cases that do not require careful control over the PSF generally should not need to filter out these objects. Unfortunately, a bug in the current version of the pipeline prevents the flag from being set on objects at CCD boundaries or containing inputs rejected due to cosmic rays or known bad pixels. This affects approximately 18% of the objects in the HSC-SSP Wide layer, though this fraction is quite sensitive to the subjective definition of the area covered by an object (we have used the detection Footprint region here; see Section 4.7). This will be fixed in a future release.

3.3.2 Safe Clipping

For PSF model coaddition to be valid, the operation used to combine all input pixels at each point on the coadd image must be strictly linear – robust estimators such as the median or sigma-clipped means cannot be used. Nonlinear estimators do not just prevent PSF coaddition from working, however; they prevent the coadd image from even having a well-defined effective PSF. Any estimator that rejects individual pixel outliers will tend to reject pixels in the cores of stars on the best-seeing exposures, and brighter stars will experience more rejection, giving them a different profile than fainter stars. It should be emphasized that this occurs even in the absence of noise, and even with extremely weak outlier rejection (e.g. clipping at $10\sigma$). All robust estimators start from the ansatz that all input values are drawn from the same underlying distribution, and convolution with different PSFs means that they are not.

To detect image artifacts that were not detected morphologically in individual visit CCD processing, then, we need a more sophisticated approach that works beyond the level of individual pixels. Our starting point is to first make both a direct mean coadd with no outlier rejection and an iterative $1.5\sigma$-clipped mean coadd, and difference them. We then run our usual maximum-likelihood detection algorithm (Section 4.7) on the difference with an extremely low $2\sigma$ threshold. The above-threshold regions (coadd difference regions; CDRs) will include both true artifacts and PSF inconsistencies due to clipping. To avoid rejecting pixels due to PSF differences we compare the CDRs to the sources detected in single-visit CCD Processing (Section 3.1, step 15). For each CDR-visit pair, we compute the “overlap fraction”: the number of pixels shared by the CDR and detected sources from that visit divided by the total number of pixels in the CDR. We consider three cases:

• •

If only one visit has an overlap fraction with a CDR greater than $0.45$, we set the CLIPPED mask bit in the CDR region on that visit; this means the CDR is probably a satellite trail, ghost, or other artifact that appears in only a single image.

• •

If no more than two visits have an overlap fraction with a CDR greater than $0.15$ and at least four visits contributed to the coadd at the position of the CDR, we set the CLIPPED mask bit in the CDR region in those visits. The low overlap fraction in this case is intended to identify artifacts that intersect real astrophysical sources, but it will miss cases where the artifact occupies a relatively small fraction of the above-threshold region.

• •

In all other cases, we ignore the CDR. Essentially all differences due to PSF inconsistencies will fall into this category, as they occur on stars or galaxies that are present in every image.

In all of the above cases, we also discard any CDRs that are more than $50\%$ occupied by pixels already masked as bad pixels, saturated, or cosmic rays in single-visit CCD Processing (Section 3.1, steps 1 and 11).

We then look for cases where the clipped area may not be large enough because a single large artifact (frequently satellite trails or optical ghosts) was detected as multiple CDRs because some pixels dipped below our threshold. We look for single-visit sources that already have more than 100 pixels marked CLIPPED; in these cases we mark the full single-visit detection as CLIPPED as well.

Finally, we build a new coadd using a direct mean that ignores all pixels marked CLIPPED in the visit images. We also set CLIPPED in the coadd mask image in these regions, where the complete set of input images was not used. We also ignore pixels marked as CR or BAD (see Table 3) in CCD Processing.

This “safe clipping” procedure is a complicated, highly heuristic algorithm that has been tuned to work with HSC-SSP data, especially in the Wide layer; the small dithers (typically $<10$\arcmin) and large number of exposures in the SSP Deep and UltraDeep layers (see Aihara et al., 2017a) make it harder for to detect and remove artifacts in those fields. The main failure mode in Wide processing is a failure to remove some satellite trails and optical ghosts that significantly overlap real bright objects; these have CDRs with large overlap fractions on multiple visits. A more general replacement algorithm that more directly utilizes image differences is currently under development.

3.4 Coadd Processing

In Coadd Processing, we detect, deblend, and measure objects on the coadd images. To maximize cross-band consistency in these measurements, we do this in five steps, alternating between processing images of each band independently and processing catalogs from all bands jointly: we detect in each band separately, merge the detections across bands, deblend and measure in each band, identify a reference band for each object, and finally perform forced measurements in each band using the reference-band positions and shapes. Each patch is processed completely independently. Figure 4 provides a summary of the processing flow.

The first step is to simply run our maximum-likelihood detection algorithm (Section 4.7) independently on the coadd image from each band with a 5$\sigma$ threshold. As in CCD Processing, we use a circular Gaussian smoothing filter matched to the RMS size of the effective coadd PSF. This generates for each band a set of above-threshold regions we call Footprints. Each Footprint contains one or more positions corresponding to peaks in the smoothed image. We also perform one more round of background subtraction (Section 4.6) to account for the new, deeper detections moving out of the definition of the background. Unlike previous background subtractions, we use a single bin for each 4k$\times$4k patch and simply subtract that constant.

In the next step, we merge the Footprints and the peak positions within them across bands. The bands are processed sequentially in what we call their priority order, which we define as irzyg for the broad-band filters, followed by the 921, 816, 1010, 387, and 515 (nm) narrow-band filters. This ordering is simply a heuristic, but it approximately goes from higher S/N bands (for typical objects) to lower S/N. The algorithm starts with the Footprints and lists of peaks from the highest-priority band (i), and proceeds as follows for the other bands.

We merge the Footprints from the previous step with any overlapping Footprints from the current band, defining the merged Footprint to contain the union of all pixels contained by all Footprints that went into the merge. Any Footprints in the current set that now overlap are similarly merged. As a result, this step may decrease the number of distinct merged Footprints, but the total area covered by all Footprints can only increase. 2. 2.

If a peak in the new band is at least $1^{\prime\prime}$ from all peaks from the previous step, we add a new peak. 3. 3.

If a peak in the new band is less than $0.3^{\prime\prime}$ from the nearest peak from the previous step, we mark the peak as having been detected in the new band while maintaining the position from the previous step. 4. 4.

If a peak in the new band is between $0.3^{\prime\prime}$ and $1^{\prime\prime}$ from the nearest peak from the previous step, it is ignored: we cannot conclusively identify the peak as the same source, and we let the higher-priority band’s set of peaks take precedence.

After the merging is complete, we attempt to cull peaks that are likely spurious detections on the outskirts of bright objects, brought above the detection threshold by flux from the bright object. We cull peaks that were detected in only one band that belong to Footprints in which there are 30 brighter peaks (in that band) or in which more than 70% of the peaks are brighter. This culling is intended to be conservative; it preserves many spurious detections to avoid any significant chance of removing real objects. The final output of this step is a catalog of merged Footprints and peaks that is consistent across bands.

We return to processing the images for each band independently in the third step. We run our deblending algorithm (Section 4.8) on each Footprint from the merged catalog. This creates an object catalog containing an undeblended parent object for every Footprint and a deblended child object for every peak (in all subsequent steps, the list of objects includes both). Because the set of Footprints and peaks is the same for all bands, we guarantee that the catalog contains the same objects in every band, but the deblending is otherwise completely independent in each band. After deblending, we run the full suite of source measurement algorithms (Section 4.9) on all objects, yielding independent measurements in every band for each object. We then resolve duplicates in the overlap regions between tracts and patches by setting the detect_is-tract-inner and detect_is-patch-inner flags on the version of each duplicate that is closest to a tract and/or patch center (respectively). This catalog of independent per-band object measurements (deepCoadd_meas) is one of our most important output data products. Note that many measurements in this catalog will be performed in a band in which the object was not even detected, with predictably poor results for algorithms with many degrees of freedom.

To address this limitation, the fourth step in Coadd Processing determines a reference band for each object by considering both the bands in which it was detected and its PSF-weighted S/N (see Section 4.9.5). The reference band for an object must be one with sufficient S/N to yield reliable measurements (or the most reliable measurements possible) even for algorithms with many degrees of freedom, but it is not simply the highest S/N band; we also wish to maximize the number of objects with the same reference band. To that end, we consider each band in the same priority order defined in the second step of Coadd Processing, and define that band to be the reference band for an object if all of the following are true:

• •

the object is marked as detected in this band;

• •

the object’s PSF S/N in this band is at least 10 or no band has PSF S/N greater than 10 and this band’s PSF S/N is at least 3 higher than the last band considered;

• •

the PSF (Section 4.9.5, Kron (Section 4.9.7), and CModel photometry (Section 4.9.9) algorithms succeeded in this band or did not succeed in any bands.

In the SSP Wide layer, 57% of objects have $i$ as their reference band, followed by 21% for $r$, 13% for $g$, 6% in $z$, and 3% in $y$. In the UltraDeep layer, approximately 1% of objects use one of the narrow-band filters as their reference band, with the other bands following approximately the same proportions as the Wide layer.

The final step is another run of the source measurement suite, but this time in forced mode: we hold all position and shape parameters fixed to the values from the previous measurement in the reference band. This ensures that the forced measurements are consistent across bands and use a well-constrained position and shape, which is particularly important for computing colors from differences between magnitudes in different bands. Consistent colors also require accounting for the different effective coadd PSFs in different bands. The PSF and CModel photometry algorithms convolve a model for the true object’s morphology with the effective PSF model at the position of the object, which maximizes S/N but depends on the correctness of the assumed underlying morphology. The seeing-matched variants of aperture and Kron photometry algorithms instead convolve the coadd images prior to measurement to match their PSFs, which reduces S/N but is insensitive to morphology. We refer the reader to Section 4.9 for the details of all of these algorithms. This final step in Coadd Processing produces the deepCoadd_forced_src dataset.

4 Algorithms

4.1 Instrument Signature Removal

The Instrument Signature Removal (ISR) component of the pipeline seeks to correct the raw images for features introduced by the CCD camera. For ISR we use the ip_isr package, which contains multiple elements that can be used in a standard or custom order. For HSC, we use the following elements, applied in order:

• •

Saturation detection: pixels exceeding the maximum correctable level are flagged as SUSPECT, while pixels exceeding the saturation level (either the full well depth or the range of the analog-to-digital converter; Miyazaki et al. in prep.) are flagged as SAT.

• •

Overscan removal: we fit (with a single rejection iteration) the overclock pixels along the rows with a 30th order Akima spline and subtract the fit from the science pixels. This leaves zero-length exposures (bias frames) generally featureless except for features manifesting as a function of column.

• •

CCD assembly: the amplifier images comprising the CCD are assembled (applying necessary flips) into a single image for each CCD.

• •

Bias correction: we subtract an average of multiple (typically around 10) overscan-corrected zero-length exposures. These overscan-corrected bias frames are generally featureless except for glowing pixels on ccd=33 and a ridge line running the length of the rows between the second and third amplifiers.

• •

Linearity correction: we apply a cubic polynomial (with coefficients determined per amplifier) to correct for non-linearity. The polynomial coefficients are determined from a fit to a series of dome flats with different exposure times.

• •

Crosstalk correction: we correct for intra-CCD crosstalk using a single correction matrix for all CCDs, and flag pixels for which the correction was significant. The correction matrix was determined during instrument commissioning from the ratios of the responses of different amplifiers to bright sources on the sky.

• •

Brighter-fatter correction: see Section 4.2.

• •

Dark correction: we subtract an average of multiple (typically 3–5) long exposures with the shutter closed, scaled by the “darktime” (time since CCD wipe). The dark is generally featureless except for glowing pixels on ccd=33.

• •

Flat correction: we divide by the dome flat frames. The dome flat is the average of multiple (typically around 10) exposures, normalized to a mean of unity over the camera. We use dome flats because they are easy to obtain, easy to repeat with the same conditions, and the incandescent lamp spectrum roughly corresponds to that of our science sources (at least, better than does the night sky). Dome flat sequences are typically acquired at least once for each available filter during a run.

• •

Defect masking: we flag as BAD pixels which we have identified as defects. The HSC CCDs are of high cosmetic quality, and there are few defects, although there are a few amplifiers (out of 448) that have gone bad since commissioning.

• •

Fringe correction: for exposures acquired using the red $y$-band or NB0921 filters, we subtract a scaled fringe frame. The fringe frame is the average (with objects masked, and rejection of outliers) of science exposures from a variety of fields observed throughout one night. The scaling is determined by measuring the average value within many ($10^{4}$) small regions over the image on the science exposure and the fringe frame.

The ideal result is an image for which the counts are proportional to flux from the sky when conditions are photometric. This ideal is not actually achieved by the ISR component because the flat-fields inevitably contain scattered light and so do not perfectly reflect the response of the camera. Furthermore, the Subaru telescope currently uses multiple dome lamps mounted on the top-end near the flat-field screen which, combined with the vignetting in the camera, means that different CCDs respond to different combinations of the lamps and therefore the flat-field pattern can change with time as the lamps age differently. These shortcomings are addressed by the Joint Calibration component (Section 3.2) of the pipeline.

4.2 Brighter-Fatter Correction

It has been recently observed for thick-depleted CCDs that the size and shape of the PSF depend upon the intensity of the source (Antilogus et al., 2014; Gruen et al., 2015). The effect has been attributed to the fact that a lateral electric field is generated as charge accumulates in the potential well. This electric field causes charges to be laterally displaced, altering the field lines and causing brighter stars to be larger and more elliptical than faint stars. This causes problems for science cases like weak lensing that rely on using bright stars to estimate the PSF model for fainter galaxies.

We have devised a correction strategy in Coulton et al. (in prep.) to reapportion the flux. In this strategy we assume that the electric field has zero curl and can thus be written as the gradient of a potential. This potential is determined from a fit to the noise correlations in flat field data. In Figure 5, we show the performance of the correction algorithm on all of the $i$-band single epoch data in the S16a internal data release by comparing the difference between the measured and model PSF sizes and ellipticities with and without the correction. We also include a plot of the difference between the corrected stellar images and the PSF model where we have stacked the stars together in magnitude bins. The correction reduces the large magnitude-dependent bias for the size and ellipticity444Ellipticity is defined throughout this paper using the distortion parameterization; see Appendix C. There is also an overall magnitude-independent bias that is not corrected. Further work in Coulton et al. (in prep.) was able to correct some of these biases by narrowing the magnitude range of selected stars and improving the recovery of the potential, but these changes were not available in the S16a internal release or PDR1. Since we are comparing to the PSF model, there may be additional errors coming from PSF modeling itself (see Section 4.3). It should be noted that these corrections are nevertheless sufficient to meet the requirements for weak lensing science as demonstrated in Mandelbaum et al. (in prep.).

4.3 PSF Modeling

To perform PSF modeling we use a version of the PSFEx software (Bertin, 2013) modified for use as a library in the LSST DM codebase. The modifications remove the necessity of running SExtractor beforehand, but we have not implemented the full set of features available with the nominal version; we limit its capability to performing a polynomial fit of the PSF as a function of position on a single CCD. The modifications allow us to separate out the task of selecting stars from modeling the PSF. The standard version performs both of these tasks, while we perform the star selection with our own algorithm and leave only the modeling to PSFEx.

We choose the set of stars to feed PSFEx based on a clustering analysis of the measured sizes of objects in each CCD (measured via second image moments; see Section 4.9.4). We use a k-means clustering algorithm on the sizes with a fixed set of four clusters to isolate the stellar locus, after first limiting the source catalog to objects brighter than 12500 counts (approximately magnitude 22.3 in $i$). The first cluster is initialized with the position of the 10th percentile from the sizes of all the measured objects. The subsequent three clusters are then initially spaced by intervals of 0.15 pixels. Each object in the catalog is assigned to the cluster whose center is closest to its measured size, and as we iterate the position of each cluster is computed as the average size of its members. We iterate until objects do not move from one cluster to another. The objects from the cluster with the smallest mean size are used as the PSF candidate sample. A final rejection is done on the sample by clipping objects outside the interquartile range centered on the median. Figure 6 shows an example of the four different clusters as a function of size and magnitude for a single CCD. The red objects are those selected as stars. As a cross-validation sample we reserve 20$\%$ of the objects from the modeling. Figure 7 shows the number of stars selected per CCD for the $i$-band visits. On average, we select $\sim$72 stars to be used in the PSF modeling on each CCD.

Postage stamps (41$\times$41 pixels) are extracted from the images around each star selected for use in PSF estimation. Since we have not run deblending yet, we mask out secondary peaks (see Section 4.7) and regions suspected to be influenced by neighboring objects. These stamps are then fed to PSFex as the initial set of stars. PSFEx is configured to use BASIS_TYPE=PIXEL_AUTO which chooses the sampling of the PSF model as the native pixel size, unless the seeing goes below 3 pixels ($\sim$0.5\arcsec), in which case the model is oversampled by a factor of two. We use BASIS_NUMBER=20, indicating that the full modeling is done on the central 20$\times$20 pixels (or 10$\times$10 when oversampling). For the flux normalization, we use an aperture flux with 12 pixel radius, which is large enough to encompass the majority of the flux and still be high S/N. We fit the PSF model to each CCD completely independently, using a second order polynomial to interpolate between stars the values of each pixel in the kernel.

To assess the quality of the PSF modeling we fit Gaussian moments as described in Section 4.9.4 to the stars and the PSF model and compare their sizes and ellipticities. Another useful quality metric is the correlation function of PSF model ellipticity residuals (Rowe, 2010),

[TABLE]

where the sum is over all pairs separated by a given distance $r$. Using the coadded PSF models described in Section 3.3.1, we can run the same test on the stars and PSF models on the coadds. To select a set of stars on the coadds we require that a given star be selected as a PSF star in $\geq$20$\%$ of the individual visits to allow for chip gaps, etc.

Figure 8 shows the size and ellipticity residuals for the visit PSF stars, visit reserved stars, and coadd stars. Figure 9 shows the ellipticity and the size and ellipticity residual patterns averaged over the focal plane. These plots show that there is un-modeled structure in the outer regions of the focal plane, especially for $e_{1}$. Increasing the polynomial order to four removes some of this structure, but we were concerned that some visits would have too few stars to fit the additional degrees of freedom, and we have hence kept the polynomial order at 2. In addition, the reserved stars are systematically worse than the stars used in the fit, suggesting potential problems in the interpolation or with over-fitting the data. The coadd PSFs also perform slightly worse than the per-visit PSFs; naively, one would expect the size of the stars on the coadd to be larger relative to the stack of the PSF models due to astrometric scatter. However, Figure 8 shows that the stars on the coadd are systematically smaller, whereas the stars on the individual visits are larger. We were unable to identify the cause of this discrepancy. Figure 10 shows the $\rho_{1}$ values averaged over all the objects in each star set and the $\rho_{1}$ values for each visit evaluated at 1\arcmin.

One serious problem we found was that the PSF modeling becomes much worse as the seeing gets better. The lower right panel of Figure 8 and the right panel of Figure 10 show this explicitly: as the PSF FWHM decreases, the size residual and $\rho_{1}$ increase significantly. One problem with these visits is that the initial processing (see Section 3.1) uses an initial guess of the PSF size of 1\arcsec, which is too large for the good-seeing visits. This causes the star selector to behave strangely and allows a large number of galaxies into the sample. Fixing this does not resolve all problems with good-seeing data, however, and we have been unable to determine the cause of the remaining issues (though we believe them to be internal to PSFEx). In the S15a, S16a, and PDR1 data releases (see Aihara et al., 2017b), we chose to reject problematic visits entirely when building the coadd. After the processing for these releases was completed, we identified more regions and visits with good seeing data that had not originally been identified, but caused significant biases in the PSF modeling. There is one particularly large region in the $i$-band GAMA09 field that users should be aware of. Mandelbaum et al. (in prep.) shows PSF diagnostics for individual fields in the HSC-SSP Wide layer, including which parts of these fields pass the PSF model quality cuts for weak lensing.

In addition to the overall poor performance in good seeing, the lower right panel of Figure 8 shows an abrupt change at $\sim 0.5$\arcsec where we change from the native pixel sampling to the oversampled basis. We attempted to improve the performance in worse seeing by using the oversampled basis, but found that the additional degrees of freedom caused other problems. In Mandelbaum et al. (in prep.), a different set of diagnostics more directly related to weak lensing analysis is performed. In future releases, we plan to improve the PSF modeling code to handle undersampled images so that we do not need to reject our best seeing data. We also plan to automatically reject CCDs with problems so as to eliminate the need to look by hand for problematic visits.

4.4 Cosmic Ray Detection

The cosmic ray (CR) removal algorithm in HSC is based upon that used successfully in SDSS (Lupton et al., 2001), although the CRs detected by the 200$\mu m$-thick CCDs employed by HSC (Miyazaki et al., in prep.) have quite different morphologies. The major sources of CRs in HSC are electrons (which have long worm-like tracks, often showing short side-branches; see the upper panel in Figure 11) and muons, which are straight and wedge-like555The wider part of the track is generated at the top of the CCDs where diffusion of carriers within the silicon after photo-conversion is more important. (lower panel in Figure 11).

The first-pass algorithm is to search for all pixels which satisfy a series of conditions:

That the candidate bad pixel $p$ not be adjacent to a saturated pixel. 2. 2.

That $p$’s intensity $z$ exceed the locally-determined background (actually the mean of pairs of neighboring pixels) by $n\sigma$ where $\sigma^{2}$ is the sky variance. We take $n=6$. 3. 3.

That the gradients near the pixel exceed the band-limit imposed by the PSF; specifically we require that the pixel be part of a peak which is sharper than the center of a star centered in a pixel. Allowing for noise, this condition becomes

[TABLE]

where $c$ is a constant, $N(z)$ is the standard deviation of $z$, $\hbox{PSF}(d)$ is the value of the PSF at a distance $d$ from the center of a star, and $\bar{z}$ is the average of two pixels a distance $d$ away from our pixel. We have found that in practice we must multiply $\hbox{PSF}(d)$ by some fiddle factor, $c_{2}<1$, to avoid flagging the cores of stars as ‘cosmic rays’. We use $c=2.5$ and $c_{2}=0.6$.

These conditions are applied sequentially to the pixel being studied using the four pairs of neighboring pixels (NS, EW, NW-SE, and NE-SW, $d=1,1,\sqrt{2},\sqrt{2}$). The candidate cosmic ray must exceed condition 2 for all four pairs of neighbors, and condition 3 for at least one pair. The thinking behind this choice is that while most cosmic rays contaminate more than one pixel, they pass through the CCD in a straight line so almost all pixels have at least one pair of uncontaminated neighbors.

Once a cosmic ray contaminated pixel is identified, its location is noted and its value is replaced by an interpolation based on the pair of pixels that triggered condition 3; the interpolation algorithm used is the same as that used for fixing sensor defects (Section 4.5). This removal of contaminated pixels as they are found makes it easier to find other pixels affected by the same cosmic ray hit.

Once the entire frame has been processed, the pixels identified individually as being contaminated by cosmic rays are assembled into cosmic ray ‘events’ of contiguous pixels. Each such event must contain more than a minimum number of electrons (not DN); we have adopted a threshold of 150$e^{-}$.

We then go through the frame again, looking at pixels adjacent to these cosmic ray events. Processing is identical, except that we set $c_{2}=0$ for these extra contaminated pixels. The major difference from the algorithm used by the SDSS Photo Pipeline is that HSC iterates, repeating this search for extra pixels 3 times. CR-contaminated pixels found during this process are, of course, removed before the next iteration.

4.5 Bad Pixel Interpolation

Linear Predictive Codes (LPC, e.g. Press et al., 2007) are a special case of Gaussian Processes that estimate the value of a signal, in our case an image, based on other data points and a knowledge of the signal’s statistical properties. Given a stationary process

[TABLE]

where $s_{i}$ is the signal, and $n_{i}$ the noise, we may predict the value of the signal $s_{p}$ at some point $p$ as

[TABLE]

it is clear that the estimator is unbiased if $\sum_{j}D_{pj}=1$. It turns out that the coefficients $D_{pj}$ are determined by the autocorrelation function of the data and their signal-to-noise.

We can use LPC to interpolate over defects in images by interpreting $p$ as the index of a bad pixel. Assuming that the image consists purely of point sources, the autocorrelation is just the square of the PSF which we approximate as a Gaussian $N(0,\alpha^{2})$, so the autocorrelation is $N(0,2\alpha^{2})$.

In order to apply this to interpolation, we set the noise to be infinitely larger in the bad columns than the good. Considering the case of infinite signal-to-noise ratio and restricting ourselves to only $5$ terms (centered on the bad column) in estimating the $D|_{p=i}$, the resulting weights for a one-dimensional interpolation with $\alpha=1$ are $\{-0.274,0.774,0.000,0.774,-0.274\}$; the bad pixel has of course a weight of $0.000$.

Image mask bits are set to indicate both that a pixel was interpolated and the reason why (e.g. cosmic ray, sensor defect, saturation, etc.) pixels. Generally only objects whose centers were interpolated (flags_pixel_interpolated_center) should be considered to have unreliable measurements.

4.6 Background Subtraction

In background subtraction, we attempt to model and subtract a smooth, slowly-varying flux distribution from an image. This flux distribution has contributions from the sky, diffuse optical ghosts from bright stars, scattered light, and sources below the detection limit. In practice the subtracted background also includes light from astrophysical objects we would ideally like to leave in the image (or at least model and subtract explicitly), such as intra-cluster light in galaxy clusters, galactic cirrus, and the extended wings of the PSF around bright stars. The fundamental challenge in background subtraction is to subtract smaller-scale patterns of scattered light while avoiding over-subtraction around bright astrophysical objects, which biases the photometry of both the bright object and its neighbors.

Our algorithm is a simple combination of averaging in spatial bins and model-fitting:

We average (unweighed mean, with iterative 3$\sigma$ clipping) pixels in 128$\times$128 bins, ignoring any pixels belonging to detected objects. We also compute the variance in each bin. 2. 2.

We fit a 6th-order 2-d Chebyshev polynomial to the average pixel values over the image (typically a CCD or patch), using the average position of the non-masked pixels as the center of a bin. We weigh each bin by its inverse variance in the fit, ensuring that heavily-masked bins cannot strongly affect the fit.

The best-fit Chebyshev polynomial is then subtracted from the original image. Because the polynomial is not required to exactly reproduce the value of any bin at its center, the order of the polynomial effectively sets the spatial scale (relative to the image size) of the features that are included in the background estimate.

Our code also supports using linear or spline interpolation to estimate the value of the background between bin centers. When interpolation is used, the effective smoothing scale of the background estimate is determined entirely by the bin size; using 512$\times$512 bins yields behavior similar to the 6th-order Chebyshev fit on the 4k$\times$2k HSC CCDs. Overall we find that Chebyshev fitting is better at avoiding over-subtraction than large-scale bins with spline interpolation, largely because it is less sensitive to the placement of bin boundaries relative to the positions of bright objects.

Because we rely on masking detected objects to prevent them from biasing the background estimate, but the detection threshold is determined relative to the background, background subtraction and source detection are not independent operations. As described in Sections 3.1 and 3.4, we subtract an initial background in CCD Processing before detecting any sources and then re-estimate and subtract the background every time we push our detection limit fainter (twice in CCD Processing, and once more in Coadd Processing). The final round of background subtraction on coadds uses a single bin for the full 4k$\times$4k patch and simply subtracts a constant; we assume most spatially-varying features have been subtracted prior to coaddition, and simply scale the coadd to account for the fact that the average flux of new detections should no longer contribute to our definition of the background.

4.7 Detection

Sources are detected on both individual exposures and coadds by convolving the image with a smoothing filter and applying a simple threshold to the result. We consider even a single pixel above threshold as a detection (unlike e.g. SExtractor; Bertin & Arnouts 1996). When the smoothing filter is the transpose of the PSF, this approach is formally optimal (as we will show in Section 4.7.1) for isolated point sources in the limit where the noise is dominated by the sky and the background is known exactly. This is at least approximately the case for the sources at the edge of our detection limit (which are mostly barely-resolved galaxies).

We use a circular Gaussian approximation matched to the RMS width of the PSF instead of the full PSF model for performance reasons. Any decrease in sensitivity due to this approximation can easily be compensated by a small decrease in the threshold, and a cut on PSF flux S/N can then be used to select true detections at the original threshold. In SSP processing we simply use the nominal threshold (generally $5\sigma$) directly; the loss in S/N is usually less than 1% in the HSC-SSP Wide layer, and is always less than 3%.

This approach is also in general non-optimal when applied naively to non-optimal coadds, as the HSC pipeline currently does – the smoothing should instead be done separately for each image and those images combined to form a detection. The corresponding loss in detection S/N is just the loss due to suboptimal coaddition reported in Section 3.3 ($1.7\%\pm 1.5\%$ in the HSC-SSP Wide).666A general-purpose optimal coadd can be constructed by smoothing the per-epoch images with their PSFs, averaging the images, and then deconvolving the coadded image to whiten the noise (e.g. Zackay & Ofek, 2015). Because detection works on smoothed images, optimal coaddition for detection does not require this deconvolution but is otherwise mathematically equivalent.

The output of our detection algorithm is a set of above-threshold regions we call Footprints. In real data, object detections may of course overlap, so each Footprint also contains the list of peaks in the smoothed image that are found within it. In general we assume each peak corresponds to a distinct astrophysical object that should be separated from its neighbors by the deblender (Section 4.8). We also grow (dilate, in the language of mathematical morphology) the Footprints by the RMS width of the PSF after detection to make them better represent the region occupied by an object; if a single smoothed pixel above threshold corresponds to a detected point source, we know that source covers approximately the area of the PSF. If growing two or more Footprints causes them to overlap, they are merged (and their peak lists concatenated), so the Footprints in the final set are non-overlapping and individually simply-connected.

4.7.1 Theory

Given an image $z(\bm{r})$, locally constant per-pixel Gaussian noise with RMS $\sigma$, and a locally constant flux-normalized PSF $\phi$, the log likelihood of a point source with flux $\alpha$ at position $\bm{\mu}$ is given by

[TABLE]

where the index $i$ runs over pixels.

This represents a Gaussian (albeit an unnormalized one) probability in $\alpha$, which we can emphasize by expanding $L$ as

[TABLE]

with $\hat{\alpha}(\bm{\mu})$ defined to null the first-order term:

[TABLE]

where $A$ is the effective area of the PSF:

[TABLE]

Like the PSF, $A$ can be considered spatially constant here even though it may vary on scales much larger than a source. Ignoring the term $L(\hat{\alpha}(\bm{\mu}),\bm{\mu})$ because it does not depend on $\alpha$, we can now read off the maximum likelihood flux as $\hat{\alpha}(\bm{\mu})$ and its variance as

[TABLE]

which yields a signal-to-noise ratio of

[TABLE]

We can thus form an image of $\nu$ by convolving the original image $z$ with the transpose of its PSF, divide by $\sigma\sqrt{A}$, and threshold at our desired significance.

To optimally detect extended sources, we could repeat the above derivation after replacing $\phi$ with a PSF-convolved model of the target source (or a sequence of such models designed to approximately cover the range of possible morphologies). We do not currently include an additional detection filter for extended objects in the HSC pipeline; a very preliminary experiment with a larger filter had a small effect on the number of detections. This possibility will be investigated more carefully in the future.

4.7.2 Temporary Background Over-Subtraction

For most groups of objects, the fact that our detection algorithm is derived for isolated objects is not a serious problem. Flux from one object will artificially increase the detection significance of neighboring pixels, but this has no effect at all on the detectability of neighboring objects that would have already been above the detection threshold. It can also lead to spurious detections due to artificially elevated noise peaks, however, and these are extremely common and problematic around the brightest stars and galaxies.

The addition of spurious objects to the catalog is not in itself a serious issue; it is relatively easy to filter out the vast majority of these spurious detections via cuts on S/N and our blendedness parameter (Section 4.9.11). The problem is how they affect their real neighbors: these peaks always “steal” some flux from their neighbors in the deblending process, biasing their photometry low, and the large number of peaks near bright objects greatly increases the chance of more catastrophic deblender failures due to linear peak alignments (see Section 4.8). It is thus imperative that we remove as many spurious peaks as possible around bright objects prior to running the deblender.

One step to remove spurious peaks has already been discussed in Section 3.4: when merging Footprints across bands, we cull the faintest peaks from large blends when they are detected in only one band. Another way to mitigate the problem is actually to do worse (in a sense) at subtracting the background: if the spatial scale of background estimation is too small, we will subtract astrophysical flux from the wings of bright objects along with the sky background, which ensures that flux does not artificially increase the detection significance. Doing so biases the photometry (and other measurements) of the bright objects, of course. It also inevitably leads to over-subtraction of the sky around bright objects, which can in turn bias the measurements of fainter objects in the over-subtracted region.

To obtain most of the benefits of aggressive background subtraction while avoiding most of the drawbacks, we perform a temporary small-scale background subtraction just before detecting peaks (but after detecting Footprints); we add the temporary background back into the image after the detection step is complete. We typically use a bin size of 64$\times$64 pixels with spline interpolation in the temporary background, as compared to a bin size of 128$\times$128 with additional polynonmial smoothing in the permanent background subtraction (see Section 4.6). This is not a particularly good model for the mostly circular (or at least elliptical) flux profiles of bright objects, of course, but it is sufficiently flexible to subtract off most of the flux in the outskirts of the brightest objects, and it is easy and fast to fit.

Unfortunately, due to a mistake in the configuration files, temporary local background estimation was not enabled in PDR1 or the internal releases it was derived from. It will be enabled in all future data releases.

4.8 Deblending

4.8.1 Overview

In a deep optical imaging survey such as the HSC-SSP, a large fraction of objects are blended with their neighbors: they have overlapping isophotes at the surface brightness threshold used for detection. In the Wide layer of the SSP, 58% of all objects are members of blends; this increases to 66% and 74% in the Deep and UltraDeep layers, respectively (and the latter has not yet been observed to full depth). Our detection algorithm (Section 4.7) represents a blend as a Footprint containing multiple peaks. The job of our deblender algorithm is to apportion the flux in that Footprint to the peaks, creating a “child image” for each peak that allows it to be measured as a distinct source. While this would ideally be done in a fully consistent way across all filters, at present our deblending algorithm processes each filter independently, with cross-band consistency enforced only by a consistent set of input peaks and Footprints (as described in Section 3.4).

We store the child images in data structures we call HeavyFootprints. These represent an image with irregular boundaries by combining a Footprint that describes the boundary with a 1-d array of child image pixel values, concatenated row-by-row. These are used by the procedure described in Section 4.9.1 to construct a larger clean measurement image for each child object.

The hard boundaries we have defined here between detection, deblending, and measurement are themselves an algorithmic choice. In crowded stellar fields this separation of concerns would be a poor choice, as deblending and measurement are essentially the same operation (fitting multiple point source models simultaneously) and detection must be iterated with subtracting the stars that were detected and modeled in previous steps (e.g. Stetson, 1987). In uncrowded fields dominated by galaxies it proved very effective as part of the SDSS Photo Pipeline (Lupton et al., 2001), however, and we have adopted both this general approach and most of the details of the SDSS algorithm for the first generation of the LSST/HSC deblender. As we will discuss in Section 4.8.3, the blending problem in HSC-SSP data is much more difficult than it was in SDSS, and it is already clear the SDSS algorithm is not adequate at typical HSC-SSP depths. It is not yet clear whether the overall approach of separating detection, deblending, and measurement is itself a limiting factor.

4.8.2 Algorithm

The fundamental ansatz in the SDSS and HSC/LSST deblenders is that objects have approximate 180-degree rotational symmetry about their peak flux position: when some part of an object’s morphology is confused with that of a neighbor, we assume we can recover it by looking on the opposite side of the peak. More formally, we define a template image $T_{i}(\bm{r})$ for peak $i$ at position $\bm{p}_{i}$ from the image $z(\bm{r})$ as

[TABLE]

Note that $2\bm{p}_{i}-\bm{r}$ is just the reflection of $\bm{r}$ to the opposite side of $\bm{p}_{i}$. The use of the minimum reflects the fact that for objects with the assumed symmetry (and in the absence of noise) the lower of the two values is the one least affected by flux from neighbors, which is always nonnegative.

We then fit a linear combination of the templates to all pixels in the blend (assuming uniform variance), yielding a best-fit scaling parameter $\alpha_{i}$ for each:

[TABLE]

It is tempting to use the scaled template $\alpha_{i}T_{i}(\bm{r})$ directly as a deblended image for peak $i$, but we can construct a better one from the relative contributions of the scaled templates to each pixel. Our final deblended images $C_{i}(\bm{r})$ are defined as

[TABLE]

This final step yields two important improvements relative to the scaled templates:

• •

it exactly conserves the flux in each pixel, and by extension, the whole blend, because the sum over all $C_{i}$ is $z$ by construction;

• •

it allows objects that are not symmetric to retain their original morphology where it is not confused with a neighbor, because the ratio in Eqn. (13) approaches one when no neighbors contribute.

In both the SDSS and HSC/LSST deblenders, templates that closely resemble the PSF model are replaced with the PSF model before fitting for the scale factors $\alpha_{i}$. This eliminates unnecessary degrees of freedom, making the final result more robust to noise and crowding.

The SDSS deblender additionally included the ability to drop templates that appeared to be too similar (i.e. their normalized dot product was too close to one). This was important for avoiding “shredding”, in which a large, morphologically complex galaxy is split into many spurious children (see, e.g. Figure 12). This is impossible in the current HSC processing flow, however, because the deblender is run independently in every band (see Section 3.4) and we have no way to guarantee that the peaks removed would be consistent across bands. Small-scale tests of this feature on HSC data suggest that it may not be as helpful for HSC data processing as it was for SDSS, overall – while it does reduce reduce shredding in large galaxies, it also merges groups of objects that appear to be distinct. A comparison of HSC and SDSS measurements of the same galaxies in Aihara et al. (2017b) provides an estimate of the magnitude of the problem; about 15% of galaxies with $i<19$ appear to be shredded in HSC processing (though this analysis depends on the assumption that the SDSS deblending and measurements are correct).

Every blend processed by the deblender begins as a single parent record in our catalog with a Footprint with multiple peaks (deblending is a no-op for sources whose Footprint contain only one peak). For each deblended peak, it adds a new child record to the catalog and attaches a HeavyFootprint containing the deblended pixel values $C_{i}(\bm{r})$. We set two fields on all records to describe these relationships:

parent:

The ID of the record that this source was deblended from, or zero if it was not blended (either because it is the original blend or was never blended).

deblend_nChild:

The number of children created by deblending this object (zero for records that represent child objects).

These can be used to identify two particularly useful categories of object:

• •

Isolated objects that were never blended have both parent=0 and deblend_nChild=0. This criterion can be used to identify a very secure (albeit small and biased) sample of objects whose measurements are definitely not affected by blending.

• •

Unblended objects that were isolated or have been deblended have deblend_nChild=0. These objects can be used for most science analyses when other cuts (on e.g. blendedness; see Section 4.9.11) are used to remove sources with especially poor deblending.

The detect_is-primary flag is set for objects that have deblend_nChild=0 as well as both detect_is-tract-inner and detect_is-patch-inner (Section 3.4), and hence can be used to select a sample of objects with no duplicate processing.

4.8.3 Results

The single biggest failure mode of the deblender occurs when three or more peaks in a blend appear in a straight line, as in Figure 13. When building the template for the middle peak, the $\min$ in Eqn. (11) will have to choose between pixels that are both affected by neighbors, yielding a poor template. In SDSS, this alignment was sufficiently rare that the overall performance of the deblender was satisfactory. In HSC data, this problem is dramatically more prevalent, simply because blends are both more common and more complex at HSC depths.

Overall, the HSC deblender performs adequately in most areas of the survey, where blending is common but not severe and most galaxies in blends are small and not obviously morphologically complex (e.g. Figures 14 and 15). But it fails dramatically in the cores of most galaxy clusters and in the neighborhood of nearby galaxies or bright stars. These constitute a significant fraction of our catalog, and they (especially in the case of galaxy clusters) are of particular importance to major HSC science goals.

While a new deblender is currently being developed for future releases, at present we provide several ways to work around the current deblender’s problems. Most of these are discussed elsewhere in the paper:

• •

We eliminate some less significant peaks before deblending, both via temporary local background over-subtraction in detection (Section 4.7.2) and when merging peaks across bands (Section 3.4).

• •

We provide a blendedness metric (Section 4.9.11) that can be used to identify objects in the kind of severe blends the deblender does not handle reliably.

• •

We run variants of the aperture photometry (for small apertures only) and Kron photometry measurements at the positions of the children using the undeblended pixels (see Section 4.9.1).

We refer the reader to Murata et al. (in prep.) for a detailed investigation of the impact of blending (and deblending) on the measurements produced by the HSC Pipeline.

4.9 Source Measurement

Source measurement here refers to the suite of algorithms that characterize the properties of individual astronomical sources. Most of these algorithms are run in both CCD Processing (Section 3.1) and Coadd Processing (Section 3.4). Our source measurement system is highly extensible; it iterates over all sources in an image and runs a configurable suite of measurement plug-ins on each. These plug-ins can be implemented outside the main suite of pipeline packages, allowing pipeline users to easily add new measurement algorithms to the system. We do not attempt to classify objects prior to measurement; while some algorithms may only be appropriate for certain types of sources, we run the same measurement algorithms on all sources.

Source measurement algorithms that perform photometry can also be run in forced mode, which fixes the position and shape parameters at values measured previously (generally on another image). This is used in the last stage of Coadd Processing to measure consistent fluxes across all bands, with positions and shapes defined in one band (see Section 3.4).

Most of the HSC Pipeline’s algorithms measure at least one of three quantities: centroids, shapes, and fluxes. While there is no requirement that a plug-in measure any of these, those that do can be configured to fulfill a measurement slot – a label that indicates a recommended version of a particular measurement, allowing it to be used by other plug-ins as an input. We have a single centroid slot and a single shape slot, but several flux slots, reflecting the fact that different kinds of photometry are not expected to yield equivalent results.

For comparisons between different photometry algorithms and characterization of the HSC Pipeline’s overall photometric performance, we refer the reader to Huang et al. (in prep.) and Aihara et al. (2017b).

The details of the measurement system are described in the following subsections. We begin by describing two aspects of the system common to all measurement algorithms: how we use deblender outputs (Section 4.9.1) and how we measure and apply aperture corrections (Section 4.9.2). Centroid algorithms are discussed in Section 4.9.3, shape algorithms in Section 4.9.4, and fluxes in Sections 4.9.5-4.9.9. Two additional algorithms for classifying sources and estimating the impact of blending on an object are discussed in Sections 4.9.10 and 4.9.11, respectively.

4.9.1 Neighbor Replacement

While the deblender discussed in Section 4.8 is responsible for apportioning flux between neighboring sources in the pixels where they overlap, we do not run source measurement plug-ins directly on these deblended pixel values. The pixels processed by the deblender are only those included in the original detection Footprint, but many measurement plug-ins require pixels beyond those regions. Some plug-ins may in fact utilize pixels belonging to Footprints other than that of the source being measured. To generate an image that permits deblended measurement beyond a source’s Footprint, we use the following procedure:

We replace all Footprints in the image with random Gaussian noise with the same variance as the original noise in those pixels. 2. 2.

We insert the deblended pixels for a particular source to be measured back into the image (replacing the noise pixels). 3. 3.

We run all measurement plug-ins on the current source. 4. 4.

We re-replace the Footprint of the current source with noise. 5. 5.

We repeat steps 2-4 for all sources.

The obvious flaw in this procedure is that the flux in pixels between Footprints is counted multiple times, as it is “assigned” to each of the sources in its neighborhood; while this flux is by definition below the detection threshold, it cannot be considered to be purely background. However, nearly all practical source measurement algorithms (and all the ones we use) down-weight pixels far from the center of the source, so this flaw is not nearly as damaging as it may appear at first glance. Moreover, counting this flux multiple times is unquestionably better than not counting it at all in the limit where sources are far apart. Finally, if multiply-counted flux beyond Footprint boundaries does prove a problem, a simple solution exists: we can simply grow the size of the Footprints (either by a fixed radius or by extending to a lower surface brightness threshold) before running the deblender.

As a partial hedge against bad deblending (galaxy “shredding” in particular), we run the suite of measurement plug-ins on both child and parent sources and include both in the resulting catalog. Child measurements use the deblender outputs to provide the per-source deblended pixel values inserted into the image in the procedure described above. Parent measurements instead interpret all pixels within a contiguous above-threshold region as belonging to a single source. As no deblending is needed in this case, the per-source pixel values are simply the original pixel values of the image. We encourage most science analysis on HSC pipeline results to use only child sources, but the parent measurements may be better for extremely bright sources whose neighbors are sufficiently faint that they can simply be ignored.

We also run variants of the Kron photometry (Section 4.9.7) and fixed aperture photometry (Section 4.9.6) algorithms for each child object on PSF-matched images that have not been deblended, to further guard against deblending failures. This is described in Section 4.9.8.

4.9.2 Aperture Corrections

Most consumers of photometry measurements assume that they represent the total flux of each source. Many practical photometry algorithms directly measure only a fraction of the flux, with most attempting to measure the same fraction for all sources. When that fraction is indeed constant, the conversion to total magnitudes is implicit in the photometric calibration: the magnitude zero-point maps our (fractional) flux measurements to the (total) magnitudes in the photometric reference catalog.

When applied to stars and other point sources, the fraction of flux measured directly by most photometry algorithms depends on the PSF. While some algorithms utilize the PSF model to limit this dependency, our PSF model only extends to $\sim\!10\sigma$, where $\sigma$ is the second-moment radius of the core of the PSF (see Section 4.9.4). While this is sufficient for centroid, shape, and morphology measurements, it ignores a non-negligible fraction of flux at large radius.

To enforce consistent flux fractions, we measure and apply aperture corrections, using the following procedure:

We run all photometry algorithms on a sample of securely classified, moderately bright unsaturated stars (typically the same stars used for PSF modeling, as described in Section 4.3). 2. 2.

For each star, we compute the ratios of each flux algorithm measurement to the measurement used for photometric calibration (by default, a 4\arcsec diameter circular aperture flux). 3. 3.

We interpolate these ratios across each CCD using 2nd-order Chebyshev polynomials. This suite of interpolated ratios is what we actually call the aperture corrections. 4. 4.

We apply the aperture corrections to per-exposure source measurements by multiplying the fluxes measured by each algorithm with the aperture correction for that algorithm. 5. 5.

We generate aperture corrections on coadds at the position of each coadd detection by averaging the aperture corrections of the coadd’s constituent CCD images at that position, using the same weights used to build the coadd. The coadd aperture corrections are then applied to coadd flux measurements.

The value of the aperture correction is typically about 0.97 (i.e. it represents a 3% change in the flux) for PSF photometry (Section 4.9.5) in the HSC-SSP Wide survey. CModel (Section 4.9.9) aperture corrections are approximately the same, while those for Kron photometry (Section 4.9.7) are typically about twice as large (1.06), in the opposite direction.

Unlike PSF coaddition (Section 3.3.1), aperture correction coaddition via a simple weighted sum is not exact. Consider two photometry algorithms $A$ and $B$, with measurements on the $i$th epoch $f_{i}^{A}$ and $f_{i}^{B}$. The measurements we would make on the coadd generated by summing images with weights $w_{i}$ are then

[TABLE]

If we use $B$ for photometric calibration, the per-exposure aperture correction for $A$ is

[TABLE]

Similarly, the true coadd aperture correction $r_{\mathrm{coadd}}$ is defined as the ratio of the coadd measurements:

[TABLE]

We can now write the denominator in terms of the per-epoch measurements of $B$ and the per-epoch aperture corrections $r_{i}$:

[TABLE]

But calibrating to $B$ also implies that (in the absence of noise), we measure the same value for $B$ on all images: $f_{i}^{B}=f_{j}^{B}=f_{\mathrm{coadd}}^{B}$ for all $i,j$. That implies

[TABLE]

The simple weighted sum of per-epoch actually used by the pipeline is the first-order Taylor expansion of this about $r_{i}=1$:

[TABLE]

In practice, the difference is typically negligible ($2.5\times 10^{-4}$ on average in the HSC-SSP Wide layer). In approximately 3% of our coadd patches (those with the largest variation in seeing among the input images) it exceeds $10^{-3}$, however, and for future data releases we will use the exact correction.

A more difficult problem in our (or any) approach to aperture corrections is galaxy photometry. Galaxy photometry measurements can miss flux due to both the same large-radius PSF flux that affects point sources as well as per-galaxy morphology differences. Because we have no source of ground truth for the latter, we have no way to independently measure the former. Our approach thus far has simply been to apply the aperture corrections for stars to galaxies as well, which should adequately correct photometry for galaxies near the size of the PSF and be better than no correction for large galaxies. In any case, we suspect (but cannot easily demonstrate) that systematic errors in galaxy fluxes due to per-source morphology and/or surface-brightness differences dominate over the errors due to incomplete correction for the wings of the PSF.

One aspect of this approach that may appear to be a serious problem (but in fact is not) is our reliance on a finite 4\arcsecaperture for calibration. This aperture is not large enough to capture all of the flux in the wings of the PSF, and hence our aperture corrections do not force all fluxes to a constant missing flux fraction. However, our final photometric calibration (see Section 3.2) includes a spatially-varying per-visit multiplicative term that can “soak up” any differences in extra flux in the wings of the PSF beyond the aperture. This is not merely a first-order correction; if the functional form of the spatial variation is sufficiently general, a multiplicative calibration correction and aperture correction are completely degenerate. This does mean, however, that our photometric calibration coefficients cannot be interpreted naively as a measure of the combination of transparency and exposure time.

4.9.3 Centroids

The first algorithm run by the HSC Pipeline’s measurement system is a centroider, as the resulting position is used by later shape and photometry algorithms. While the HSC Pipeline contains multiple centroid algorithms and implementations, we have thus far only used one in production: an approximate maximum-likelihood algorithm first developed for the SDSS Photo Pipeline (Pier et al., 2003).

Basically, we note that in the limit that the noise is dominated by the background the maximum-likelihood estimate of the position of a star is the peak of the image correlated with the PSF (approximated as a Gaussian, as we do in Section 4.7); we then use parabolic interpolation in this smoothed image to find the position of the likelihood peak. When the object is significantly larger than the smoothing filter we bin the original image and re-smooth (effectively doubling the smoothing length); this is repeated until the smoothing filter and the object size roughly agree.

Unlike the SDSS algorithm, we do not de-bias centroid measurements to account for asymmetry in the PSF, as this is largely degenerate with the astrometric solution as long as the PSF asymmetry is approximately constant at the smallest spatial scales on which the astrometric solution varies.

4.9.4 Shapes

HSC Pipeline’s “shape” algorithms produce an ellipse that characterizes the extent of the object without correcting for its convolution with the PSF. This is obviously not the only way to characterize an object’s morphology or size (or even the most useful for science), but most other characterizations are highly algorithm-specific and hence hard to interpret generically. More importantly, this PSF-uncorrected ellipse is what is typically most useful to downstream flux measurement algorithms. Algorithms (both those that produce this ellipse and those that do not) are free to define additional outputs, and these frequently do include PSF-corrected shapes such as those used for weak gravitational lensing.

All of the HSC Pipeline’s current shape plug-ins compute the second moments of the image of the source about its centroid, using a Gaussian weight function. Given image $z(\bm{r})$ (with position $\bm{r}$ defined about its centroid), the raw moments (a $2\times 2$ matrix) are:

[TABLE]

where $\bm{C}$ is the covariance matrix of the Gaussian weight function. Using a weight function is necessary to keep the noise in the measurement from diverging, but it also biases the result. We can see this by inserting a Gaussian image with known moments $\bm{Q}$ into Eqn. (21), and taking the limit of small pixels (so the sums become integrals):

[TABLE]

To correct for this bias (only approximately, because real sources are not Gaussian), we can solve Eqn. (24) for $\bm{Q}$:

[TABLE]

This provides a way to correct the measured raw moments $\bm{M}$ given the moments of the weight function $\bm{C}$, which may be a predefined constant or a quantity computed iteratively.

The three unique elements of the symmetric $2\times 2$ matrix $\bm{Q}$ parameterize an ellipse, and this is the form used to record the outputs of our shape algorithms. See Appendix C for information on how this parameterization relates to others.

Two of the three shape algorithms currently implemented in the HSC Pipeline use adaptive Gaussian moments, in which the weight function’s moments $\bm{C}$ are iteratively matched to the (corrected) measured moments $\bm{Q}$. This makes the correction formula Eqn. (25) particularly simple: $\bm{Q}=2\bm{M}$. More importantly, it produces the highest (for a Gaussian weight function) S/N measurement, assuming that the iterative procedure converges. This is not guaranteed, especially for low surface-brightness objects, where noise can make the raw moments singular. Nearly singular moments are in some respects even more problematic, as these can yield to spuriously large ellipses for small objects, with no straightforward indicator of failure, which can then compound by increasing the size of the weight function. Our two adaptive moments algorithms, SdssShape and HsmMoments, differ primarily in the heuristics they apply to stabilize the iteration and test for convergence. SdssShape is a reimplementation of the algorithm used in the SDSS Photo Pipeline (Lupton et al., 2001), while the HsmMoments algorithm is simply a thin wrapper around the HSM (Hirata & Seljak, 2003) implementation in GalSim (Rowe et al., 2015). The two algorithms fail on approximately the same number of sources (16% for SdssShape and 17% for HsmMoments), though this can depend on the details of how the sample is otherwise selected. The two algorithms fail on a different population of sources (only approximately 8% of sources are not measured successfully by either algorithm) and hence we run both in HSC-SSP data release productions. Most of these failures occur when the moments matrix $\bm{Q}$ becomes singular. In addition to adaptive moments, the HSM algorithms also include several approaches to estimating PSF-corrected ellipticities for the purpose of estimating shear due to weak gravitational lensing. These are described more fully in Hirata & Seljak (2003) and Mandelbaum et al. (in prep.).

The HSC Pipeline’s third shape algorithm, SimpleShape, uses a circular Gaussian weight function with fixed (but configurable) radius. This makes iteration unnecessary, removing one source of catastrophic failure. Singular moments are still possible, of course, and the resulting measurement is more noisy than the adaptive methods. We use SimpleShape primarily when processing data from HSC’s 8 out-of-focus wavefront sensors as described in Furusawa et al. (in prep.).

The CModel photometry algorithm (Section 4.9.9) also measures multiple ellipses for each source (corresponding to different profiles), but these ellipses are approximately PSF-corrected (assuming analytic forms for the source morphology) and hence cannot be used as a “slot” shape that feeds other algorithms.

4.9.5 PSF Photometry

The algorithm for PSF photometry given a measured centroid can be derived from two different starting premises: maximum-likelihood fitting and matched-filter signal processing. In the former, we consider an image $z_{i}$ with Gaussian noise $\sigma_{i}$ and assume a one-parameter model that simply multiplies the PSF model $\phi_{i}$ (shifted to the same sub-pixel centroid as the source, which we consider fixed) with an amplitude parameter $\alpha$. The likelihood is then

[TABLE]

Differentiating the logarithm of $P$, setting the result to zero, and some algebra yields the familiar solution for linear least squares:

[TABLE]

In the matched-filter version, we simply compute the inner product of the data (the image) and the filter (the PSF model, shifted to the centroid of the source, as before), and divide by the effective area of the filter:

[TABLE]

When the per-pixel noise is constant, as is the common case for faint objects where noise from the sky dominates, these two derivations agree. When the per-pixel noise is not constant, as is the case for bright objects whose own photon noise dominates, both versions yield the correct result, as can be seen if we insert $z=\alpha\phi_{i}$ into Eqns. (27) and (28). The noise properties of the measurements are not the same, however; the variance in the maximum-likelihood measurement is

[TABLE]

while for the matched filter estimate it is

[TABLE]

These once again reduce to the same estimate when sky noise dominates, but the maximum-likelihood measurement has lower noise when it does not. In the bright-object regime, however, systematic errors from incorrect PSF modeling are a much bigger concern than noise. To illustrate the effect such an error has on these estimates, we consider a true image

[TABLE]

where $\epsilon$ represents the deviation of our model from the true PSF, along with noise of the form

[TABLE]

where $b$ is the level of the background (before it is subtracted). If we work in units of electrons, this is the natural combination of Poisson noise from both the source and the background in the limit where both are large enough that we can approximate the Poisson distributions as Gaussian. Inserting these into Eqn. (27) and Eqn. (28) yields, respectively,

[TABLE]

and

[TABLE]

Both formulae yield the incorrect flux, of course, but the damage done to the maximum-likelihood estimate is much more severe: it no longer responds linearly to the true flux. The matched-filter estimate, in contrast, remains linear, and can be corrected completely by aperture corrections if the error in the PSF model is constant or smoothly-varying over the image.

As the above discussion suggests, PSF photometry in the HSC Pipeline always uses the matched-filter formula Eqn. (28), and considers the per-pixel noise only when estimating the uncertainty in the flux. We also treat the centroid as fixed at the position measured by the centroider (see Section 4.9.3) and do not attempt to fold centroid uncertainties into the PSF flux uncertainty; even when the S/N is low, additional uncertainty introduced by the centroid is subdominant. We use Lanczos interpolation to shift the PSF model image to the centroid of the source, and when the PSF model is defined on an oversampled pixel grid relative to the image, we perform the shift before down-sampling to the image pixel grid.

4.9.6 Fixed Aperture Photometry

The HSC Pipeline measures a set of aperture fluxes with fixed radii. The radii are configurable, but by default are set to am approximately logarithmically-spaced sequence from 1-23\arcsec in diameter. These fluxes are by default not aperture-corrected, because their intent is to measure a crude, PSF-uncorrected radial profile, not a total magnitude.

However, not applying aperture corrections to these fluxes also puts them in a different photometric system from our other photometry. As discussed in Section 4.9.2, we use the 4\arcsec diameter aperture for photometric calibration, and hence define it to be the total flux. Apertures larger than 4\arcsec thus typically appear to have more flux than the “total” flux measured in the aperture-corrected photometric system. In the future, we plan to address this by estimating the curve-of-growth from the brightest stars and interpolating this in the same way we interpolate our current aperture corrections; we could then tie the radial profile to the aperture corrections at infinity. This would require using saturated stars to get a sufficiently large signal-to-noise ratio in the largest aperture fluxes, which would in turn require more sophisticated algorithms for aperture photometry. An approach that worked in SDSS but has not yet been implemented for HSC/LSST is to divide each aperture both radially and azimuthally into segments of annuli; this allows discrepant segments (due to e.g. bleed trails) to be removed with outlier rejection and more robust uncertainties to be estimated from the statistics of the segments.

Instead, aperture photometry measurements in the HSC Pipeline are simply direct integrals over the circular area. For small radii, they are indeed integrals, not simple sums – we use the approach of Bickerton & Lupton (2013) to exactly integrate over sub-pixel regions, avoiding the discretization artifacts that can otherwise affect small-aperture fluxes. We just use naive sums for larger apertures, where the discretization noise is negligible compared to photon noise and the performance difference between the two approaches is significant.

4.9.7 Kron Photometry

The Kron radius (Kron, 1980) of a galaxy is defined as

[TABLE]

where $z$ is the surface brightness and $r$ is the distance to the center of the galaxy.

Ideally, the integration area $\cal A$ would extend to infinity, but in practice the integrals do not converge very well (the numerator diverges logarithmically if the profile falls as $r^{-3}$).

Once $R_{\mbox{\tiny Kron}}$ has been determined, the Kron flux is the flux enclosed in a region of radius $f_{\mbox{\tiny Kron}}R_{\mbox{\tiny Kron}}$; a common choice is $f_{\mbox{\tiny Kron}}=2.5$.

If the object is elliptical things are a bit more complicated. We can still measure $R_{\mbox{\tiny Kron}}$, but its meaning is unclear; HSC uses 2nd-moment shapes (Section 4.9.4) to set the axis ratio and position angle, and then defines the Kron ellipse to have area $\pi R_{\mbox{\tiny Kron}}^{2}$.

The HSC algorithm begins by defining the area, ${\cal A}$, over which we evaluate $R_{\mbox{\tiny Kron}}$. We start with the adaptive moments defined in Section 4.9.4 and take the initial $\cal A$ to be the ellipse with the same shape and orientation but with area $\texttt{nSigmaForRadius}\times\sqrt{\det\bm{Q}}$ (by default nSigmaForRadius = 6). The adaptive moments code fails for some objects; in these cases we take the initial ellipse to be that of the PSF model.

With this definition of $\cal A$ we evaluate Eqn. (35) and adjust $\cal A$ to have area $\pi\times\texttt{nSigmaForRadius}^{2}\times R_{\mbox{\tiny Kron}}^{2}$. In theory we could iterate, but in practice we have found that a single update is sufficient.

There are some failure modes to keep track of; in particular if $\cal A$ reaches the edge of the image we give up and set kron_flags_edge. If either the numerator or denominator in Eqn. (35) is non-positive (as can happen as $z$ is not guaranteed to be positive in the presence of noise) or the estimated $R_{\mbox{\tiny Kron}}$ is less than the PSF’s $R_{\mbox{\tiny Kron}}$ we set kron_flags_usedPsfRadius and set $R_{\mbox{\tiny Kron}}$ to the Kron radius of the PSF.

With a value of $R_{\mbox{\tiny Kron}}$ in hand we can measure the Kron flux as the flux included in an elliptical aperture with determinant radius (see Appendix C) $\texttt{nRadiusForFlux}\times R_{\mbox{\tiny Kron}}$; by default $\texttt{nRadiusForFlux}=2.5$ as suggested above. The only error that can be encountered at this stage is when this aperture extends outside the image, in which case kron_flags_edge is set and no flux is reported. The error on the Kron flux is purely that due to the pixel noise (i.e. we do not attempt to propagate the error in estimating the radius, which would be correlated with the errors in the pixel values).

4.9.8 Afterburner Photometry: PSF-Matched Apertures Without Deblending

As we discussed in Section 4.8, our deblender performs adequately when the number of objects in a blend and the amount of overlap between them are both small. In dense regions such as galaxy clusters, however, it can produce catastrophically bad results that make it essentially impossible to accurately photometer the deblended children. As an alternative, our afterburner777The “afterburner” name is a historical artifact referring to when these algorithms were run relative to the rest of the pipeline, but it has caught on in the HSC collaboration and is used in other HSC papers (e.g. Aihara et al., 2017b), and hence it is helpful to use it here as well for clarity. photometry algorithms perform measurements at the position of each child object on the original (not deblended) image. We run both Kron photometry and fixed circular aperture photometry in this mode.

The measurement image used by the afterburner algorithms is not the original image, however; using aperture photometry to measure consistent colors requires images with the same PSF in all bands. We thus approximately match the coadd images in each band to predefined circular Gaussian PSFs with full-width at half maximums (FWHMs) of 3.5, 5.0, and 6.5 pixels (0.58, 0.84, 1.09\arcsec), by convolving with a Gaussian matching kernel derived by assuming the original PSF is Gaussian. More precisely, if the target PSF is defined as a Gaussian with covariance matrix

[TABLE]

and the 2nd moments matrix (see Section 4.9.4) of the PSF model at the position of the source to be measured is $\bm{Q}_{\mathrm{psf}}$, then the matching kernel we use is simply a Gaussian with covariance $\bm{Q}_{\mathrm{match}}=\bm{Q}_{\mathrm{target}}-\bm{Q}_{\mathrm{psf}}$. When $\det{\bm{Q}_{\mathrm{match}}}\leq 0$ (i.e. the original PSF is larger than the target PSF) we do not match the image at all; the fact that we match to a sequence of increasingly larger target PSFs ensures that at least one set of afterburner measurements should exist for nearly all objects. The matching is performed locally for each object, so we assume the PSF (and hence the matching kernel) is constant over the pixels to be used in the measurement. We never produce a PSF-matched image that covers the full area of a patch.

Because the matching process makes the PSF broader, it actually increases the degree to which objects are blended, and hence the only measurements even approximately unaffected by blending are those measured with a very small aperture. These will of course only measure the colors of the cores of objects, and even they rely on blending not being too severe, but they nevertheless appear to be our most reliable measurements of galaxy colors (Tanaka et al., 2017; Aihara et al., 2017b). Larger fixed aperture and Kron photometry measurements on these PSF-matched images are likely to be strongly affected by blending, but can still be quite useful for objects that were not blended or whose flux dominates that of their neighbors. They are also useful when the blending involves galaxies with similar colors at the same redshift, as is the case in the galaxy cluster cores where our deblender’s performance is at its worst.

4.9.9 CModel Photometry

Our primary algorithm for general galaxy photometry is a modified version of the approach used in the SDSS Photo Pipeline (Abazajian et al., 2004). It fits multiple PSF-convolved galaxy models in a sequence designed to approximate a bulge-disk decomposition or Sérsic model (Sérsic, 1963) fit while minimizing the number of degrees of freedom in any individual fit. Roughly speaking, the sequence is:

Fit an elliptical model with an exponential profile (a Sérsic profile with index $n=1$) to the source, with both the ellipse parameters and the amplitude free. 2. 2.

Fit a de Vaucouleurs model (de Vaucouleurs 1948; a Sérsic profile with index $n=4$) with the same parameters free (allowing the fit to yield different values for these parameters). 3. 3.

Fit both models simultaneously, keeping their ellipses fixed at the results from the previous two fits, and allowing only the two amplitudes to vary.

The amplitudes from all three of these fits are useful flux measurements. The last is not well-motivated statistically, but worked well empirically in SDSS as a crude linear interpolation between exponential and de Vaucouleurs models. This model is typically a worse fit to real galaxy morphologies than a Sérsic profile, but it can be fit much more robustly and efficiently. This is especially true for low S/N and/or poorly-resolved galaxies.

Because all of these models are convolved with the PSF of the image being measured, CModel fluxes are PSF-corrected to the extent the models are sufficient to describe the underlying galaxy morphology. Given the nearly limitless range of morphologies real galaxies have, the models are obviously inadequate at some level, but in practice we have seen no indication that biases due to this sort of mismatch are significantly affecting the CModel colors. This may be masked by other failures in CModel, however; as our photometry algorithm with the most degrees of freedom, it has a tendency to respond very poorly to even minor deblending failures.

We also expect any model biases in CModel photometry to be much less significant for colors (ratios of fluxes) than the fluxes themselves. It is difficult to decide how far out to extrapolate an analytic galaxy profile even for nearby objects that are being modeled with significant human intervention, and these poorly-constrained tails can have a large influence on the integrated flux. The best we can hope for in fully-automated fitting of faint, poorly-resolved sources is some degree of consistency, and even this is difficult to verify because we lack any kind of ground truth. Colors need not correspond to total fluxes, however; a color measurement that systematically gives too much or too little weight to some stellar populations in a galaxy will still correspond to some physically meaningful star formation history as long as the measurements are consistent across bands (note that this is true even in the presence of large population gradients). We can almost guarantee this in forced photometry: we use the same model in all bands, allowing only the amplitude parameters to vary. If we had the same PSF in every band or complete models, consistent colors could be guaranteed; because the PSF is not the same in all bands and we correct for these differing PSFs only by convolving our (imperfect) models, colors may still be slightly inconsistent.

While the CModel algorithm produces PSF-corrected ellipse shapes, the fitting algorithm has been tuned for photometry, and we do not recommend using these ellipses for weak lensing shear measurements. Similarly, half-light radii derived from these ellipses are probably not reliable for objects close to the PSF size, and even the radii of well-resolved objects are more sensitive than fluxes to our simplistic choice of models.

The technical details of the CModel algorithm are described more fully in Appendix B. This includes discussion of a serious bug in PDR1 that strongly distorts the structural parameters measured by CModel (without significantly damaging – and in some cases probably improving – measured colors for most objects).

An analysis of CModel performance on simulated galaxies in Huang et al. (in prep.) indicates that CModel performs well for photometry down to $i\sim 25$ despite the algorithm’s aforementioned limitations. As noted by Aihara et al. (2017b), however, CModel colors are outperformed by the PSF-matched small-aperture “afterburner” photometry (Section 4.9.8) in yielding a lower scatter in the red sequence of galaxies in clusters (Aihara et al., 2017b, Figure 19), which is one of the only tests of galaxy photometry we have that does not rely on simulations. The afterburner photometry also seems to yield better photometric redshifts (Tanaka et al., 2017). It is unclear whether this is due to CModel’s sensitivity to deblender failures (which the afterburner fluxes avoid), the number of degrees of freedom in the algorithms (CModel has many more free parameters that must be constrained), or astrophysical differences between the stellar populations in the cores of galaxies (which the afterburner measure) and those in the outskirts (which are also included in the CModel photometry).

4.9.10 Star/Galaxy Classification

Separating stars from galaxies in images relies chiefly on a measure of how different an object’s morphology is from that of a point source. Our primary measure of this is called extendedness, and it is defined as the difference between the PSF magnitude and the CModel magnitude. For a star, the expected value of this difference is zero because the scale size of the best-fit galaxy model will approach zero before PSF convolution and hence the PSF-convolved model will be very close to the PSF itself. For resolved galaxies the expected value of the difference in magnitudes is significantly different from zero, and in most cases it is positive, because PSF photometry typically underestimates the flux of extended objects.

To illustrate extendedness-based star/galaxy classification, Figure 16 shows a diagram of the PSF-CModel magnitude difference vs. CModel magnitude in $i$ for objects detected in the COSMOS UltraDeep field. These objects were matched to the HST/ACS catalog of Leauthaud et al. (2007). With significantly higher resolution (0.1\arcsec FWHM PSF) and comparable depth, we can effectively use this HST/ACS catalog to provide “ground truth” star/galaxy labels: the stars in the figure are colored blue, and the galaxies red. Figure 17 shows the purity and completeness of the samples of stars (right panel) and galaxies (left panel) obtained from the cuts in extendedness denoted by the horizontal lines in Figure 16. Extendedness seems to work well down to $i\sim 24$, where the ratio of stars to galaxies shrinks dramatically and most galaxies become smaller than the PSF. The classification_extendedness column in our catalogs is obtained by setting a hard cut (at $m_{\mathrm{psf}}-m_{\mathrm{cmodel}}=0.0164$) on this magnitude difference, but users of the catalog may be able to obtain a more suitable classifier for their science by making a different cut.

To produce more accurate star/galaxy classifications, we combine colors measured from CModel fluxes with extendedness via supervised machine learning, using the HST/ACS catalog of Leauthaud et al. (2007) to provide labels for the training set. We estimate probability densities for both stars and galaxies using Extreme Deconvolution (XD; Bovy et al., 2011), an algorithm based on mixtures of Gaussians that estimates the probability density from multivariate data, taking into account measurement errors. The probability densities are then combined with a flat prior to compute a posterior probability of being a star.

Because the true probability of begin a star is a strong function of magnitude (largely because galaxy counts rise faster with magnitude than star counts), we divide the catalog 4 magnitude bins in $i$: [18, 22], [22, 24], [24, 25], and [25, 26]. An XD fit in colors and extendedness is then produced in each bin independently. T he left panel in Figure 18 shows that the posterior probabilities inferred from the XD density distributions do behave like a probability; the right panel plots the posterior probability vs CModel magnitude in $i$ for objects detected in the COSMOS UltraDeep layer. Figure 19 shows the purity and completeness of samples of stars (right panel) and galaxies (left panel) obtained from cuts in the posterior. While most internal data releases have included a pstar column with the posterior probability of belonging to the star class estimated with this procedure, the machine learning classifier has not yet been run on PDR1 and hence the pstar column is not yet available; it will be added in a future incremental release.

4.9.11 Blendedness Metrics

Given the central importance (and fallibility) of the deblender, it is useful to provide some measure of the extent to which each object is blended. The metric we have developed for this purpose essentially measures the fraction of the total flux in the neighborhood of a source that belongs to its neighbors. This neighborhood does not have a hard boundary; it is defined by a Gaussian weight function with its ellipse taken from the source’s shape (see Section 4.9.4).

More formally, given an image $z(\bm{r})$ and the source moments matrix $\bm{Q}$, the Gaussian-weighted flux is

[TABLE]

where $k(\bm{Q})$ is an ellipse-dependent normalization we can safely ignore because it will cancel out in the next step. With $z_{c}$ the child image (with neighbors replaced with noise as per Section 4.9.1) and $z_{p}$ the parent (original) image, the raw blendedness $b_{\mathrm{raw}}$ is defined as:

[TABLE]

Note that we use the same shape $\bm{Q}$ in both measurements, and that this shape is itself derived from the child image.

To make sure the quantity is well-behaved in the presence of noise, we also compute what we call the absolute blendedness $b_{\mathrm{abs}}$, in which we use the absolute value of the image along with an attempt to correct the bias this introduces:

[TABLE]

where $d(z)$ is a de-biasing term derived in Appendix D.

Blendedness behaves as one would expect, given its name: for isolated objects, $z_{c}=z_{p}$, and $b=0$; for a faint object close to a much brighter object, $b$ approaches unity. It is also obviously dependent on the actual performance of the deblender, which is responsible for determining $z_{c}$ from $z_{p}$ for each child. This is most problematic when the deblender assigns too much flux to a child, which both artificially increases that child’s signal-to-noise ratio and decreases its blendedness, making it less likely that cuts on either will successfully remove poor measurements. While this is a common deblender failure mode, we have nevertheless found that a conservative blendedness cut can effectively mitigate the problem, because it occurs most frequently for faint objects on the outskirts of extremely bright objects, where $b$ is high even if it is underestimated. Mandelbaum et al. (in prep.) include only objects with $\log_{10}b<-0.375$, for instance, in the HSC-SSP weak lensing shear catalog. We refer the reader to Murata et al. (in prep.) for extensive tests of the blendedness parameter.

5 Future Work and Conclusions

Like the SSP survey, the HSC Pipeline is a work in progress. It is in nearly all respects a state-of-the-art optical imaging data reduction pipeline, and it has already enabled an impressively broad array of scientific discoveries (as demonstrated by the rest of this PASJ special volume). But it has also in some contexts been the limiting factor in taking full advantage of the exquisite data produced by the HSC camera and the Subaru telescope, highlighting the algorithmic challenges posed by extremely deep, ground-based optical imaging data (and in some cases the challenges posed by excellent seeing).

The biggest algorithmic challenge for HSC going forward is clearly deblending. Our algorithm, which was quite adequate at the brighter magnitude limit of the SDSS and is still probably the most sophisticated and ambitious approach currently used in production by any many survey, is simply not good enough to handle the more severe blending prevalent in deeper ground-based data. A completely new algorithm is currently in development that should give us the ability to experiment with much more varied (and frequently more restrictive) constraints than rotational symmetry, such as monotonicity (in the radial profile) and consistency across bands.

PSF modeling is also a serious concern, as we do not currently obtain PSF models accurate enough for weak lensing in the best-seeing images from PSFEx, and we do not believe this can be fixed with configuration changes. Instead of a purely empirical model that uses a pixel basis for PSF images and polynomial interpolation, our next-generation PSF modeling involves attempting to separate the optical and atmospheric components of the PSF in the spirit of Davis et al. (2016). This should provide better interpolation and more accurate PSF images (especially in good seeing, where the contribution from the optics is more important).

These issues also highlight the role of HSC as a pathfinder for LSST; LSST’s 18,000 $\mathrm{deg}^{2}$ wide survey will exceed the SSP Wide depth before its fourth data release, and its deepest fields may exceed SSP UltraDeep depths even sooner. As a result, its blending problems will ultimately be even more severe. Physical PSF models will also be important for LSST; its requirements on PSF accuracy are significantly tighter, and LSST’s optical design will produce larger chromatic PSF effects (Meyers & Burchat, 2015), which can probably only be modeled well with a more physically-motivated approach. Overall, by essentially commissioning an early prototype of LSST pipelines on HSC data and using the results for cutting-edge research, the HSC-SSP has revealed and clarified (and will continue to reveal and clarify) the most challenging algorithmic problems for LSST in a way that “data challenges” cannot; there is no substitute for having a collaboration of engaged scientists actually use the outputs of a pipeline to do research.

Efforts to improve many other aspects of the HSC pipeline are also underway, including new algorithms for background subtraction, artifact masking, and galaxy photometry. The many straightforward bugs noted in this paper that affected the processing of PDR1 are our highest priority; most of these are already fixed in the current version of the pipeline and will be included in the next internal release. With significant HSC data now available to the full LSST DM team from PDR1 and the algorithmic improvements made on the HSC Pipeline now integrated back into the LSST codebase, we expect this collaboration to be even more productive. The HSC Pipeline will continue to improve as the SSP survey progresses, advancing the state-of-the-art in image processing algorithms to enable cutting-edge science from one of the highest-quality datasets in astronomy.

{ack}

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.

The HSC Pipeline and SSP processing would not have been possible without both the software developed by LSST Data Management and the continuing support of DM scientists and engineers who are not members of the HSC collaboration. We thank the LSST Project for making their code available as free software at http://dm.lsstcorp.org.

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory.

This paper is based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan.

HM is supported by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. RM is supported by the Department of Energy Early Career Award Program. SM is supported by the Japan Society for Promotion of Science grants JP15K17600 and JP16H01089.

Appendix A Joint Calibration Mathematical Details

In this appendix we describe the detailed mathematical form of the model fit by our joint astromeric and photometric calibration pipeline (Section 3.2). In the following, super- or subscript $s$, $e$, and $c$ denote variables that vary for each star, exposure, or CCD chip, respectively.

The astrometric model consists of several composed transforms. The first of these relates the focal plane coordinates $(u,v)$ (in units of nominal pixels, i.e. 15$\mu m$) for a star to the CCD coordinates $(x,y)$ via the CCD offsets $X_{c}$ and $Y_{c}$ and rotations $\theta_{c}$ by

[TABLE]

We map also celestial coordinates $(\alpha,\delta)$ for each star to gnomonic projected coordinates $(\xi,\eta)$, with the projection centered at sky position $(p_{\alpha},p_{\delta})$:

[TABLE]

For sources not matched to a reference catalog object we solve for $\alpha^{s}$ and $\delta^{s}$ as part of the fit. For matched sources $(\alpha,\delta)$ are held fixed to the reference catalog positions.

We relate $(u,v)$ to $(\xi,\eta)$ via 2-d polynomials with coefficients $a$ for $\xi$ and $b$ for $\eta$; the differences between the model and data in $(\xi_{,}\eta)$ are thus

[TABLE]

and the sum of differences yields the full $\chi^{2}$ for the astrometric fit:

[TABLE]

We assume all centroid uncertainties are the same.

To minimize, we linearize each residual in all of the free parameters ($a^{e}$, $b^{e}$, $X_{c}$, $Y_{c}$, and $\theta_{c}$, as well as $\alpha^{s}$ and $\delta^{s}$ for sources not matched to the reference catalog), yielding the matrix equation

[TABLE]

with $\bm{r}$ the vector of parameters defined block-wise as

[TABLE]

and $\bm{J}$ the matrix of first derivatives with blocks

[TABLE]

with (dropping most $s,e$ superscripts for brevity)

[TABLE]

Note that $P$, $a$, and $b$ are flattened into vectors in the block matrix form; the ordering of the two polynomial dimensions and the index over exposures is unimportant as long it is consistent for all three quantities.

We initialize all free parameters with the values from the single-visit astrometric fit performed in CCD Processing (Section 3.1), and fit by iteratively solving Eqn. (60).

In the photometric fit, we solve for polynomial coefficients $f_{j,k}$, per-exposure magnitude offsets $dm^{e}$, and per-source true magnitudes $m_{0}^{s}$, with the latter held fixed at the reference catalog magnitude for sources matched to the reference catalog. We minimize

[TABLE]

where $T_{n}$ is the nth-order Chebyshev polynomial of the first kind. As with the astrometric fit, we initialize with the values from CCD Processing and proceed by iteratively solving a linearized version of the difference between the true magnitude $m^{s}_{0}$ and the per-exposure measured magnitudes $m^{s,e}$.

Both the astrometric and photometric fits are performed independently for each tract, using inputs from only those CCDs that overlap the tract; CCDs that overlap multiple tracts are fit multiple times. In the HSC-SSP Wide layer, each iteration of the astrometric fit involves an LU factorization of an approximately 6000$\times$6000 square matrix. The matrix for each photometric fit iteration is approximately 2500$\times$2500 (again in the Wide layer).

Appendix B CModel Algorithm Details

The elliptical exponential and de Vaucouleurs models fit in the CModel algorithm have the precise form

[TABLE]

where $\bm{r}$ is the pixel position (relative to the source centroid), $\alpha$ is the amplitude parameter, $k(r)$ is the radial profile, and $\bm{S}$ is the “generating transform” (see Appendix C) of the isophotal ellipse that contains half the flux of the model. The matrix $\bm{S}$ depends on the ellipse parameters. We use the $\{\eta_{1},\eta_{2},\ln r_{\mathrm{tr}}\}$ parameterization of the ellipse (again, see Appendix C) because any triple in $\mathbb{R}^{3}$ is a valid ellipse in this parameterization and hence we do not need a numerical optimizer that can handle parameter constraints.

The profiles $k(r)$ themselves are approximations to the exponential profile

[TABLE]

and de Vaucouleurs profile

[TABLE]

formed from a linear combination of Gaussians with different widths, as described in Hogg & Lang (2013). More specifically, we use their six-Gaussian lux and eight-Gaussian luv profiles. These are matched to the modified exponential and de Vaucouleurs profiles used in SDSS model fitting, which have softened cores and truncated wings (at four or eight half-light radii, respectively). Both of these modifications allow the multi-Gaussian approximations to better approximate the ideal profiles in the range of radii over which they can be best constrained by our data, and the truncation in particular helps to prevent spuriously large fluxes due to extrapolation.

The primary advantage of these multi-Gaussian approximations is that they allow fast convolution, if the PSF is approximated by a compatible analytic function. Rather than follow Hogg & Lang (2013) and Sheldon (2014) in using sums of Gaussians for the PSF, however, we use a sum of two Gauss-Hermite (“shapelet”) expansions with different scales. As demonstrated by Bosch (2010), using multiple Gauss-Hermite expansions provides significantly more flexibility than using a single one in representing functions with broader-than-Gaussian tails while maintaining fast convolution properties. We fit the Gauss-Hermite PSF approximation to our image-based PSF model (either the single-epoch PSFEx model or the effective coadd PSF described in Section 3.3.1) only after evaluating the image-based model at the position of the source.

When fitting for ellipse parameters in the separate exponential and de Vaucouleurs fits (steps 1 and 2 in Section 4.9.9), we combine the likelihood with a simple Bayesian prior on the ellipse parameters to form the objective function that we optimize. The prior is flat in amplitude and position angle, and is separable in $\ln r_{\mathrm{tr}}$ and $\eta=\sqrt{\eta_{1}^{2}+\eta_{2}^{2}}$. The prior on $\ln r_{\mathrm{tr}}$ is piecewise, transitioning from flat to a Student’s $t$ distribution at $\sim 0.3$\arcsec ($\ln\frac{r_{\mathrm{tr}}}{\mathrm{arcsec}}=-1$), which is also the peak of the Student’s $t$. The prior on the flux is intentionally flat, and we do not condition the prior for the ellipse parameters on the flux. Doing so may be worthwhile, but is potentially dangerous; flux measurements are expected to respond strictly linearly to changes in the true underlying brightness, and a prior on flux could disrupt that scaling.

The distribution in ellipticity and the outer piece of the distribution in radius are designed to be broader than the empirical distributions of galaxy sizes and ellipticities in galaxies in the HST COSMOS catalog of Leauthaud et al. (2007), as shown in Figure 20. Unfortunately, a bug in the relative weighting of the likelihood and prior in the version of the code used in the first HSC-SSP public data release resulted in the effective prior being much stronger than intended. The most significant result of this bug is a very sharp cutoff in the radius prior at the $\sim 0.37$\arcsec transition point. Ellipticities are also biased low (more circular), but this has a small effect on the flux. The radius bias yields an artificial “pile-up” of galaxies at $m_{\mathrm{psf}}-m_{\mathrm{cmodel}}\sim 0.6$ as shown in Figure 21. In at least one context, however, this biased flux is actually be an improvement; even minor deblend failures can cause CModel radii for faint galaxies to be spuriously large, and the sharp effective prior suppresses this effect. Because all bands are approximately equally affected by this bug, galaxy colors from CModel do not appear to be biased nearly as strongly as fluxes. As a result, we have been reluctant to fix this bug until we have another approach in place to address CModel sensitivity to deblender failures. A radius prior that is more informative but matched to the true distribution of faint galaxy radii is the most likely candidate. This would probably have to be at least conditioned on flux to allow bright (and typically larger) galaxies to see a different radius prior, and as we have noted above, it is not clear whether such a prior can yield a linear response of the measurement to changes in intrinsic brightness. Given the dependency of the measured flux on the best-fit radius and the difficulty of determining unbiased galaxy sizes, it is clear our current CModel fluxes are already at some level nonlinear. We believe this is true of all existing approaches to measuring galaxy photometry that infer a model or aperture size from the data.

Appendix C Ellipse Parameterizations

Many of the algorithms for photometry and shapes in the HSC Pipeline rely on different parameterizations of ellipse, including the second moments discussed in Section 4.9.4, the Kron photometry algorithm discussed in Section 4.9.7, and the CModel algorithm (Section 4.9.9 and Appendix B). We describe here the relationships between these parameterizations and our conventions for them for clarity and reader convenience.

For our purposes, all ellipse parameterizations are three numbers, as we generally consider the position of the ellipse (another two numbers) as a distinct quantity. The most intuitive parameterization is a combination of semi-major radius ($A$), semi-minor axis ($B$), and position angle ($\theta$), which we define here as the angle from the x-axis to the semi-major axis, measured counterclockwise.

In addition to its intuitive simplicity, the $\{A,B,\theta\}$ parameterization makes it easy to compute what we will call the ellipse’s “generating transform” $\bm{S}$: a linear transform that maps the unit circle to the ellipse:

[TABLE]

The generating transform is not unique; it may be multiplied on the right by any orthogonal matrix. In addition to the form we have defined above, a symmetric form where this matrix is the inverse of the rotation on the left is also particularly common. As described in Appendix B, the ellipse-generating transform can be used to map a 1-d radial profile to an elliptical model.

Most algorithms based on the moments of images (see Section 4.9.4) instead use the three unique elements of a symmetric positive definite $2\times 2$ matrix we will call $\bm{Q}$. For the 1$\sigma$ contour of an elliptical 2-d Gaussian, $\bm{Q}$ is just the covariance matrix. The three unique elements of $\bm{Q}$ are most easily related to the $\{A,B,\theta\}$ parameterization by the fact that $\bm{S}$ is the matrix square root of $\bm{Q}$:

[TABLE]

This also means that $A^{2}$ and $B^{2}$ are the eigenvalues of $\bm{Q}$, and the columns of the rotation matrix for $\theta$ are its eigenvectors.

Weak gravitational lensing frequently uses definitions of ellipticity comprised of two real numbers that can be interpreted as a single complex number. The magnitude of the complex ellipticity is some function of $A$ and $B$, and its phase is $2\theta$. Common functions to define the magnitude include the conformal shear

[TABLE]

the distortion

[TABLE]

and what we will call the shear888Both $\delta$ and $g$ are confusingly referred to as both ellipticity and eccentricity in the literature; we use “distortion” and “shear” because, to our knowledge, $\delta$ is never called shear and $g$ is never called distortion.

[TABLE]

All ellipticities in this paper use the distortion definition. The conformal shear $\eta$ is defined on $(0,\infty)$, and hence its real and imaginary parts

[TABLE]

are defined on $(-\infty,\infty)$, making them very useful parameters for working with unconstrained optimizers or probability distributions with infinite extents. $\delta$ and $g$ are defined on $(0,1)$ and their real and imaginary parts on the unit circle, making them less useful for these purposes, but they are easier to connect to gravitational lensing observables.

To represent a full ellipse using any of these definitions of complex ellipticity, we need to add a radius parameter. Either $A$ or $B$ may be used in this role, but we tend to instead use two definitions related to the determinant and trace of $\bm{Q}$

[TABLE]

The determinant radius $r_{\mathrm{det}}$ is proportional to the square root area of the ellipse, but we prefer the trace radius $r_{\mathrm{tr}}$ (or its logarithm) for fitting, as it remains nonzero (finite) when $B$ approaches zero but $A$ does not, which is common when fitting noisy, unresolved objects as in the CModel algorithm.

Appendix D Correcting for Bias in Absolute Blendedness

In order to compute a blendedness parameter that is well behaved in the presence of noise, we use the absolute value of the image, rather than the image itself while computing the blendedness parameter in Eqn. (41). For objects which are well above the noise level, this procedure should have no effect on the blendedness parameter, but using the absolute value will bias the blendedness for objects which are close to the noise level. We correct for such a bias using the following formalism.

Consider an object which results in a photon count $z(\bm{r})$ in a pixel located at ${\bm{r}}$ in the image plane. Under the assumption that the photon counts are Gaussian-distributed with mean $\mu$ and dispersion $\sigma$, the total integrated flux is given by

[TABLE]

where $A$ is the area over which the integral is performed, and we have assumed a uniform profile for the object. Instead, if we integrate $|z|$ over the same area, we obtain a biased answer

[TABLE]

This bias can be corrected if instead of $|z|$, we integrate $|z|-d$, where $d$ is given by

[TABLE]

Because that the objects in our images will have a non-uniform profile, we will use $d(z)$ instead of a constant $d$ in Eqn. (41).

Given the photon count $z(\bm{r})$ and the dispersion $\sigma$, we thus need to estimate $\mu$ in order to compute $d(z)$. We compute $\mu(z,\sigma)$ by imposing the condition that the expected photon count be greater than 0. Therefore,

[TABLE]

The original implementation (used in PDR1) of the bias correction for blendedness used an incorrect formula, resulting primarily in a remapping of the blendedness values for low S/N objects (though it was not strictly a remapping). Because all users of the blendedness parameters have relied on calibrating its values directly via e.g. via simulations (Murata et al., in prep.) instead of relying on its theoretical interpretation, we do not expect this bug to significantly affect science analyses.