Pan-STARRS Pixel Analysis: Source Detection and Characterization

The basic psphot analysis is divided into several major stages, as listed below.

• 1.  
  Image preparation. Load data, characterize the image background, load or construct variance and mask images.
• 2.  
  Initial source detection. Smooth, find peaks, measure basic properties with focus on the point sources to measure the PSF.
• 3.  
  PSF determination. Select PSF candidates, perform model fits, build PSF model from fits, select best PSF model class.
• 4.  
  Bright-source analysis. Fit sources with PSFs, determine PSF validity, subtract PSF-like sources, fit non-PSF model(s), select best model class, subtract model.
• 5.  
  Faint source analysis. Detect low-level sources, measure properties (aperture or PSF).
• 6.  
  Aperture corrections.Measure the curve of growth, spatial aperture variations, and background-error corrections.
• 7.  
  Output. Write out sources in selected format, write out difference image, variance image, etc., as selected.

In addition to this basic sequence, additional analysis steps may be performed. An "extended source" analysis mode is available to measure photometry and morphology of galaxies and other resolved sources. Forced photometry may be performed for both point-like and extended sources. A special mode is available for the photometry of sources detected in difference images. These different modes are discussed in their own sections below.

Table 1 lists the types of analyses performed by psphot, specifying which of the psphot use cases performs the given analysis. The table also provides a reference to the section of this paper in which the analysis is described. Not all analyses are relevant to all sources in all images. The table identifies those cases where the analyses are applied to only a subset of all sources.

psphot is highly configurable. Users may choose via the configuration system which of the above analyses are performed. This is useful for testing but also allows for specialized use cases. For example, the PSF model may already be available from external information, in which case the PSF modeling stage can be skipped.

Ultimately, all measurements of individual astronomical sources from psphot are reported in one of the tables in the PSPS database. As discussed in detail in Paper VI, measurements from individual exposures are available from the Detection table. Measurements of objects in the stacked images are stored in one of several Stack... tables, while the "forced" measurements from individual warp images are stored in tables beginning with ForcedWarp.

During the psphot analysis, there are a wide variety of conditions that are identified by the analysis software. As part of the output data for each detected source, two fields that encode these conditions as bit values in the two 32 bit integers are provided. The following two tables list the individual bit values in these two fields. These informational and warning bits are described in more detail later in this article.

Table 2 lists the flags recorded in the output field FLAGS. When data from psphot is loaded into a DVO database (Magnier et al. 2020b), these values are stored in the field Measure.photFlags and exposed in the public database (PSPS; Flewelling et al. 2020) in the fields Detection.infoFlag, StackObjectThin.XinfoFlag (where X is one of grizy), and ForcedWarpMeasurement.FinfoFlag. Table 3 lists the flags recorded in the output field FLAGS2. When data from psphot are loaded into a DVO database (Magnier et al. 2020b), these values are stored in the field Measure.photFlags2, and they are exposed in PSPS in the fields Detection.infoFlag2, StackObjectThin.XinfoFlag2 (where X is one of grizy), and ForcedWarpMeasurement.FinfoFlag2.

Table 3. Detection Flag Values #2 Reported by psphot

Flag Name	Flag Value	Description
PM_SOURCE_MODE2_DIFF_WITH_SINGLE	0x00000001	Diff source matched to a single positive detection
PM_SOURCE_MODE2_DIFF_WITH_DOUBLE	0x00000002	Diff source matched to positive detections in both images
PM_SOURCE_MODE2_MATCHED	0x00000004	Source generated based on another image
PM_SOURCE_MODE2_ON_SPIKE	0x00000008	$\gt 25 \%$ of (PSF-weighted) pixels land on diffraction spike
PM_SOURCE_MODE2_ON_STARCORE	0x00000010	$\gt 25 \%$ of (PSF-weighted) pixels land on star core
PM_SOURCE_MODE2_ON_BURNTOOL	0x00000020	$\gt 25 \%$ of (PSF-weighted) pixels land on burntool
PM_SOURCE_MODE2_ON_CONVPOOR	0x00000040	$\gt 25 \%$ of (PSF-weighted) pixels land on convpoor
PM_SOURCE_MODE2_PASS1_SRC	0x00000080	Source detected in first pass analysis
PM_SOURCE_MODE2_HAS_BRIGHTER_NEIGHBOR	0x00000100	Peak is not the brightest in its footprint
PM_SOURCE_MODE2_BRIGHT_NEIGHBOR_1	0x00000200	${{flux}}_{{\rm{n}}}/({r}^{2}{{flux}}_{{\rm{p}}})\gt 1$
PM_SOURCE_MODE2_BRIGHT_NEIGHBOR_10	0x00000400	${{flux}}_{{\rm{n}}}/({r}^{2}{{flux}}_{{\rm{p}}})\gt 10$
PM_SOURCE_MODE2_DIFF_SELF_MATCH	0x00000800	Positive detection match is probably this source
PM_SOURCE_MODE2_SATSTAR_PROFILE	0x00001000	Saturated source is modeled with a radial profile
PM_SOURCE_MODE2_ECONTOUR_FEW_PTS	0x00002000	Too few points to measure the elliptical contour
PM_SOURCE_MODE2_RADBIN_NAN_CENTER	0x00004000	Radial bins failed with too many NaN center bin
PM_SOURCE_MODE2_PETRO_NAN_CENTER	0x00008000	Petrosian radial bins failed with too many NaN center bin a
PM_SOURCE_MODE2_PETRO_NO_PROFILE	0x00010000	Petrosian not built because radial bins missing
PM_SOURCE_MODE2_PETRO_INSIG_RATIO	0x00020000	Insignificant measurement of Petrosian ratio
PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN	0x00040000	Petrosian ratio in the zeroth bin (likely bad)
PM_SOURCE_MODE2_EXT_FITS_RUN	0x00080000	We attempted to run extended fits on this source
PM_SOURCE_MODE2_EXT_FITS_FAIL	0x00100000	At least one of the model fits failed
PM_SOURCE_MODE2_EXT_FITS_RETRY	0x00200000	Trailed asteroid model fit was retried with new window
PM_SOURCE_MODE2_EXT_FITS_NONE	0x00400000	All of the model fits failed

Notes. These are saved in output catalogs as the field FLAGS2, in the DVO database as Measure.photFlags2, and in the public database as Detection.infoFlag2, StackObjectThin.XinfoFlag2 (where X is one of grizy), and ForcedWarpMeasurement.FinfoFlag2.

^a Not used for DR1 or DR2.

Download table as:  ASCIITypeset image

The first step is to prepare the image for detection of the astronomical sources. We need three separate images: the measured flux (signal image), the corresponding variance image, and a mask defining which pixels are valid and which should be ignored. The signal and variance images are represented internally as 32 bit floating point values. The variance and mask images may either be provided by the user, or they may be automatically generated from the input image, based on configuration-defined values for the image gain, read noise, saturation, and so forth. Within the IPP analysis, we normally use images that are equivalent to the digital numbers (scaled by the detrend images), but as long as the variance image is constructed in a consistent fashion, psphot can use images in electron, calibrated flux units, or other conventions (though this would require some tuning of configuration parameters). For the function-call form of the program, the flux image is provided in the API, and references to the mask and variance are provided in the configuration information. As in the standalone C program, the variance and mask may be constructed automatically by psphot.

The mask is represented as a 16 bit integer image in which a value of 0 represents a valid pixel. Each of the 16 bits define different reasons a pixel should be ignored, listed in Table 4. This allows us to optionally respect or ignore the mask depending on the circumstance. For example, in some cases, we ignore saturated pixels completely while in other circumstances, it may be useful to know the flux value of the saturated pixel. In addition, the mask pixels are used to define the pixels available during a model fit; those which should be ignored for that specific fit are "marked" by setting a special bit (MARK=0x8000). The initial mask, if not supplied by the user or library calls, is constructed by default from the image by applying three rules: (1) Pixels that are above a specified saturation level are marked as saturated. The level is specified by the camera format keyword CELL.SATURATION, which may specify a value or define a header keyword which in turn specifies the value in the image header. In the case of PS1 PV3, the header keyword MAXLIN specifies the saturation level for each chip (see Waters et al. 2020). (2) Pixels that are below a user-defined value (CELL.BAD = 0 for PV3) are considered unresponsive and masked as dead. (3) Pixels that lie outside of a user-defined coordinate window are considered nondata pixels (e.g., overscan) and are marked as invalid. (psphot recipe keywords XMIN, XMAX, YMIN, YMAX, all set to 0 for PS1 PV3—invalid pixels were specified for PS1 PV3 with a supplied mask image (see Waters et al. 2020).

The library functions used by psphot understand two types of masked pixels: "bad" and "suspect." Bad pixels are those that should not be used in any operations, while suspect pixels are those for which the reported signal may be contaminated or biased, but may be usable in some contexts. For example, a pixel with poor charge transfer efficiency is likely to be too untrustworthy to use in any circumstance, while a pixel in which persistence ghosts have been subtracted might be useful for detection or even analysis of brighter sources. Table 4 lists the 16 bit values used for PS1 mask images, along with their description (see Waters et al. 2020 for additional information).

An important point to note is that psphot does not attempt to interpolate or replace bad pixel values in the images before processing. The GPC1 images have quite extensive masking due to both defects and natural gaps between detectors and amplifier regions. On average, roughly 71% of the full usable field of view is covered with valid pixels (see Paper III for more discussion). Any attempt to interpolate bad pixels would be quickly overwhelmed by these extensive regions. Rather than attempt to fill in the bad pixels, we rely in the PS1 PV3 processing on the fact that regions on the sky were observed many times. Thus, it should be noted that model-fitting measurements (which can naturally ignore masked pixels) should generally be more reliable than aperture-like measurements for single exposures. Aperture-like measurements from the stacks do not suffer from this masking issue. See also the discussion of the PSF_QF and PSF_QF_PERFECT parameters for judging the impact of masking on a particular source (Section 4.5.3).

The variance image, if not supplied, is constructed by default from the flux image using the configuration supplied gain and read noise values to calculate the appropriate Poisson statistics for each pixel. The parameters are determined based on the camera format keywords CELL.GAIN and CELL.READNOISE, which in the case of PS1 PV3 refer to the header keywords GAIN and RDNOISE. In this case, the image is assumed to represent the readout from a single detector, with well-defined gain and read-noise characteristics. This assumption is not always valid. For example, if the input flux image is the result of an image stack with a variable number of input measurements per pixel (due to masking and dithering), the variance cannot be calculated from the signal image alone. It is necessary in such a case to supply a variance image that accurately represents the variance as a function of position in the image.

Some image processing steps introduce cross-correlation between pixel fluxes. An obvious case is smoothing, but geometric transformations that redistribute fractional flux between neighboring pixels also introduces cross-correlations. In the noise model, it is necessary to track the impact of the cross-correlations on the per-pixel variance. In the general case, this would require a complete covariance image, consisting of the set of cross-correlated pixels for each image pixel. Because a typical smoothing or warping operation may introduce correlation between 25 and 100 neighboring pixels, the size of such a covariance image is prohibitive.

Before sources are detected in the image, a model of the background is subtracted. The image is divided into a grid of background points with a spacing defined by the psphot recipe values BACKGROUND.XBIN, BACKGROUND.YBIN, set to 400 pixels (∼100'') for PV3. Superpixels of size BACKGROUND.XSAMPLE, BACKGROUND.YSAMPLE (2 × 2 for PV3) times larger than this spacing are used to measure the local background for each background grid point, thus oversampling the background spatial variations. In the interest of speed, a subset of IMSTATS_NPIX (10,000 for PV3) randomly selected unmasked pixels in these regions are used to determine the background. The background value for each superpixel is determined by fitting a Gaussian distribution to the histogram of pixel values.

If the image were empty of stars and only contained flux from a uniform background sky, we would expect the distribution to be Poisson distributed and in general in a high-enough signal range to be essentially Gaussian. We fit a symmetric Gaussian to all histogram bins within 15% of the peak bin value to determine the mean and standard deviation values for the background.

If, however, the sky is not empty of stars or other sources, and we have correctly masked the large majority of nonresponsive pixels, then we expect the flux distribution of the pixels to be asymmetric with a Gaussian core representing the sky and a tail to the high end representing the pixels with astronomical source flux contributions. We would like to determine the mean of the underlying Gaussian without suffering bias from the stellar flux. We thus perform a second Gaussian fit using an asymmetric subset of the histogram pixels, fitting those histogram bins that are left of the peak but for which the bin value is greater than 25% of the peak bin, or right of the peak but only using those bins for which the bin value is greater than 50% of the peak bin value.

If the fit to the asymmetric lower fraction of the curve is less than the symmetric fit but greater than the above lower bound of the full symmetric fit, then the lower fraction value is kept as the true mean sky value for this superpixel. Table 5 shows a comparison of this technique to several other methods to measure the sky background using simulated data with a range of stellar densities. The stellar density listed in the table is the number of stars per square degree at the 5σ detection limit in the lowest-density image. In our simulations, we find that as the stellar density rises to values typical in the Galactic plane regions, this technique results in a more accurate estimate of the background, though it still overestimates the background compared to the truth.

Table 5. Comparison of Background Measurement Methods

Density	True	Image	Image	Gauss	psphot
${\mathrm{log}}_{10}({{\rm{\deg }}}^{-2}$)	Sky	Mean	Median	Fit	Value
4.2	202.8	203.3	202.8	202.8	202.9
4.7	202.8	204.9	203.1	203.0	203.0
5.2	202.8	210.6	204.0	203.5	203.5
5.7	202.8	233.9	207.4	205.4	205.3
6.2	202.8	300.9	219.7	211.2	210.6
6.7	202.8	534.6	286.2	242.8	233.9

Note. Backgrounds were measured for simulated images with the given stellar density (at the low-density detection threshold) and known background level. The psphot technique is less biased at high stellar densities.

Download table as:  ASCIITypeset image

Bilinear interpolation is used to generate a full-resolution image from the grid of background points, and this image is then subtracted from the science image. The background image and the background standard deviation image are kept in memory from which the values of SKY and SKY_SIGMA are calculated for each source in the output catalog. For more details of the background subtraction, see the discussion in Section 3.11 of Waters et al. (2020).

Because the subtraction of the sky model suppresses larger-scale structures, features such as large galaxies that are comparable to the superpixel size are adversely affected by the subtraction. Photometry for galaxies larger than ∼30'' is unreliable as a result. The superpixel size used for the sky model in the PV3 analysis was chosen as a compromise between the need to follow bright features with small spatial scales and the desire to measure photometry of galaxies of sizes up to at least 30''. Features that we wished to suppress include both astronomical sources, such as bright nebulosity and the wings of bright stars, and non-astronomical sources, such as moonlight and other scattered-light sources. In some contexts, we have used a finer spacing for the background model, such as in the dedicated analysis of the photometry of the Andromeda galaxy, where we are only interested in stellar sources, and the analysis is otherwise badly affected by the background from this galaxy.

4.4.1. Peak Detection

The initial source detection step is focused on finding and identifying the brighter point sources. The goals are twofold: (1) to select sources that can be used to model the PSF and (2) to subtract the brighter sources so that fainter sources may be found throughout the image .

The sources are initially detected by finding the location of local peaks in the image. The flux and variance images are smoothed with a small circularly symmetric kernel using a two-pass 1D Gaussian. The smoothed flux and variance images are combined to generate a significance image in signal-to-noise units, including correction for the covariance, if known. At this stage, the goal is only to detect the brighter sources, above a user-defined signal-to-noise ratio (S/N) limit (PEAKS_NSIGMA_LIMIT = 20.0 for PV3). A maximum of PEAKS_NMAX (5000 for PV3) are found at this stage.

For an image with a Gaussian PSF of the same size, this method would represent the optimal detection algorithm, equivalent to a matched filter. At this stage, our goal is simply to detect the brighter sources, so the exact size and shape of the PSF are not critical. The detection efficiency for the brighter sources is not strongly dependent on the form of this smoothing function. Instead, our goal with the smoothing kernel is to reduce our sensitivity to pixel-to-pixel fluctuations in the location of the peak of the sources in the image.

The local peaks in the smoothed image are found by first detecting local peaks in each row. For each peak, the neighboring pixels are then examined, and the peak is accepted or rejected depending on a set of simple rules. The rules are defined so that we choose a unique set of peaks that are not immediately adjacent to other peaks. First, any peak that is greater than all eight neighboring pixels is kept. Any peak that is lower than any of the eight neighboring pixels is rejected. Any peak that has the same value as any of the other eight pixels is kept if the pixel X and Y coordinates are greater than or equal to the other equal-value pixels. This last rule means that a flat-topped region will result in peaks at the maximum X and Y corners of the region.

We use the 9 pixels that include the source peak to fit for the position and position errors. We model the peak of the sources as a 2D quadratic polynomial, and use a very simple biquadratic fit to these pixels. We use the following function to describe the peak:

$$\begin{eqnarray*}&&f(x,y)={C}_{00}+{C}_{10}x+{C}_{01}y+{C}_{11}{xy}+{C}_{20}{x}^{2}+{C}_{02}{y}^{2},\end{eqnarray*}$$

and write the chi-square equation:

$$\begin{eqnarray*}&&{\chi }^{2}=\displaystyle \sum _{i,j}{\left({F}_{i,j}-f(x,y)\right)}^{2}/{\sigma }_{i,j}^{2}.\end{eqnarray*}$$

By approximating the error per pixel as the Poisson error on just the peak, pulling that term out of the above equation, and recognizing that the values $X,Y$ in the 3 × 3 grid centered on the peak pixel have values of only 0 or 1, we can greatly simplify the chi-square equation to a square matrix equation with the following values:

$$\begin{eqnarray*}\left(\begin{array}{cccccc}9 & 0 & 0 & 0 & 6 & 6\\ 0 & 6 & 0 & 0 & 0 & 0\\ 0 & 0 & 6 & 0 & 0 & 0\\ 0 & 0 & 0 & 6 & 0 & 0\\ 6 & 0 & 0 & 0 & 6 & 4\\ 6 & 0 & 0 & 0 & 4 & 6\end{array}\right)\left(\begin{array}{c}{C}_{00}\\ {C}_{10}\\ {C}_{01}\\ {C}_{11}\\ {C}_{20}\\ {C}_{02}\end{array}\right)=\left(\begin{array}{c}\displaystyle \sum {F}_{i,j}\\ \displaystyle \sum {F}_{i,j}x\\ \displaystyle \sum {F}_{i,j}y\\ \displaystyle \sum {F}_{i,j}{xy}\\ \displaystyle \sum {F}_{i,j}{x}^{2}\\ \displaystyle \sum {F}_{i,j}{y}^{2}\end{array}\right).\end{eqnarray*}$$

Inverting the 3 × 3 matrix terms for C₀₀, C₂₀, and C₀₂, the location of the peak is determined from the minimum of the biquadratic function above and is given by

$$\begin{eqnarray}&&{x}_{\min }=({C}_{11}{C}_{01}-2{C}_{02}{C}_{10}){D}^{-1},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{y}_{\min }\,=\,({C}_{11}{C}_{10}-2{C}_{20}{C}_{01}){D}^{-1},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&D=4{C}_{20}{C}_{02}-{C}_{11}^{2}.\end{eqnarray} \tag{ 3 }$$

The resulting peak position, (${x}_{\min },{y}_{\min }$), is used as the default starting coordinate for the source. Later in the psphot analysis, improved measurements of the source positions are calculated as discussed below.

4.4.2. Footprints

The peaks detected in the image may correspond to real sources, but they may also correspond to noise fluctuations, especially in the wings of bright stars. psphot attempts to identify peaks that may be formally significant but are not locally significant. It first generates a set of "footprints," contiguous collections of pixels in the smoothed significance image above the detection threshold (PEAKS_NSIGMA_LIMIT). These regions are grown by a small amount to avoid errors on rough edges—an image of the footprints is convolved with a disk of radius FOOTPRINT_GROW_RADIUS (3 pixels for PV3). Peaks are assigned to the footprints in which they are contained (note by construction all peaks must be located in a footprint because the peaks must be above the detection threshold).

For any peak that is not the brightest peak in that footprint, it is possible to reach the brightest peak by following a sequence of the highest valued pixels between the two peaks. The lowest pixel along this (potentially meandering) path is the key col for this peak (as used in topographic descriptions of a mountain). If the key col for a given peak is less than FOOTPRINT_CULL_NSIGMA_DELTA (4.0 for PV3) sigmas below the peak of interest, the peak is considered to be locally insignificant and removed from the list of possible detections (see Figure 1). If more than one such path is possible, the path with the highest key col is used for this test. In the vicinity of a saturated star, the rule is somewhat more aggressive as the flat-topped or structured saturated top of a bright star may appear as multiple peaks with highly significant cols between them. However, this is an artifact of the proximity to saturation. Sources for which the peak is greater than 50% of the saturation value require the col to also be a fixed fraction (5%) of the saturation below the peak to avoid being marked as locally insignificant.

Zoom In Zoom Out Reset image size

Sometimes, it is useful to know if a source has a near neighbor that may be affecting the photometry. Three flag bits are used to identify such possible situations. Peaks that are not the brightest peak within a single footprint have the flag bit PM_SOURCE_MODE2_HAS_BRIGHTER_NEIGHBOR set. This is a fairly common situation. We also define the following ratio to compare the flux of the bright source to the flux of a neighbor scaled by intervening area: $R={f}_{{\rm{n}}}/{r}^{2}{f}_{{\rm{p}}},$ where ${f}_{{\rm{n}}}$ is the flux of the brightest neighbor in the footprint, ${f}_{{\rm{p}}}$ is the flux of the source of interest, and r is the separation between the two sources. If $R\gt 1$, the flag bit PM_SOURCE_MODE2_HAS_BRIGHT_NEIGHBOR_1 is set. If $R\gt 10$, the flag bit PM_SOURCE_MODE2_HAS_BRIGHT_NEIGHBOR_10 is set.

4.4.3. Centroid and Higher-order Moments

Once a collection of peaks has been identified, a number of basic properties of the sources related to the first, second, and higher moments are measured. These moments can be used for a crude classification of the sources. As discussed below, the second moments are used to select candidate stellar sources to be used in modeling the PSF and to identify "cosmic rays" and extended sources. The radial moment is used in the measurement of the Kron magnitudes (Kron 1980). The higher-order moments are provided primarily for image quality diagnostics.

In order to measure the moments, it is necessary to define an appropriate aperture in which the moments are measured. We also apply a "window function," down-weighting the pixels by a Gaussian, centered on the object, with size ${\sigma }_{w}$ chosen to be large compared to the PSF size, ${\sigma }_{\mathrm{PSF}}$. This window function reduces the noise of the measurement of the moments by suppressing the noisy pixels at high radial distance as well as by reducing the contaminating effects of neighboring stars. The choice of ${\sigma }_{w}$ and the aperture is an iterative process: for a given value of ${\sigma }_{w}$, the PSF stars will have a measured value of the PSF size, ${\sigma }_{\mathrm{PSF}}^{{\prime} },$ which is different from the true value due to the effect of the window function. The measured value of the PSF size will be biased high or low depending on both the signal-to-noise of the source and the size of the window function compared to the true PSF size.

These effects are illustrated in Figure 2 using simulated data. An image was generated with a PSF model matching the radial profile of the PS1 PSF model with ${\sigma }_{\mathrm{PSF}}$ corresponding to an FWHM of 14. For bright stars, as the window function ${\sigma }_{w}$ is increased, the measured FWHM rises from an initially underestimated value to meet the truth value. For faint stars, the measured value of the FWHM is initially underestimated as well. However, as the value of ${\sigma }_{w}$ increases, the measured FWHM for faint stars rises and then overshoots the truth value, while the scatter increases. Thus, for large values of ${\sigma }_{w}$, the result is both a poorly estimated FWHM for the image and a trend with the S/N of the star. We attempt to minimize the scatter and trends with instrumental magnitude at the cost of overall bias.

Zoom In Zoom Out Reset image size

In a real image, we do not know the true value of the PSF size. If we simply choose a very large window function and rely on bright stars, our estimate of the PSF size will be quite noisy. Compounding this problem are the two additional facts that (1) we do not know which are the real stars (as opposed to bright galaxies or possible image artifacts) and (2) the brighter stars are themselves subject to additional biases due to saturation and other nonlinear effects (c.f., "the Brighter-Fatter" effect; Antilogus et al. 2014; Gruen et al. 2015). To make a robust choice for ${\sigma }_{w}$, we choose a value such that the measured value of ${\sigma }_{\mathrm{PSF}}^{{\prime} }$ is 65% of ${\sigma }_{w}$. The resulting second-moment values are biased somewhat low (∼75% of the truth value for the PS1 PSF profile) but are relatively unbiased as a function of brightness.

To choose the value of ${\sigma }_{w}$, we try a sequence of values spanning a range guaranteed to contain any reasonable seeing values. The values are specified in the psphot recipe as PSF.SIGMA.VALUES and have the following values for PS1 PV3: (1, 2, 3, 4.5, 6, 9, 12, 18) pixels ∼(026, 051, 077, 115, 154, 23, 31, 46). For each of these ${\sigma }_{w}$ values, we then select candidate PSF stars based on the distribution of the measured ${\sigma }_{\mathrm{PSF}}^{{\prime} }$ in the two principal directions: ${\sigma }_{x,x}$ and ${\sigma }_{y,y}$ (see Section 4.5.2, below). For each test value of ${\sigma }_{w}$, we determine the ratio ${\rho }_{\sigma }=\tfrac{{\sigma }_{x}+\sigma y}{2{\sigma }_{w}}$, i.e., the ratio of the window size to the observed PSF size. We interpolate to find a value of ${\sigma }_{w}$ for which ${\rho }_{\sigma }$ is expected to be 0.65. We use an aperture with a radius of $4{\sigma }_{w}$ to select the pixels for the measurement of the moments.

Once ${\sigma }_{w}$ has been determined, moments are measured as defined below:

$$\begin{eqnarray}&&{x}_{0}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}{w}_{i}({f}_{i}-{s}_{i}){x}_{i},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{y}_{0}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}{w}_{i}({f}_{i}-{s}_{i}){y}_{i},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{M}_{{xx}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}{w}_{i}({f}_{i}-{s}_{i}){\left({x}_{i}-{x}_{0}\right)}^{2},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{M}_{{xy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}{w}_{i}({f}_{i}-{s}_{i})({x}_{i}-{x}_{0})({y}_{i}-{y}_{0}),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{M}_{{yy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}{w}_{i}({f}_{i}-{s}_{i}){\left({y}_{i}-{y}_{0}\right)}^{2},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{M}_{{xxx}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}}({f}_{i}-{s}_{i}){\left({x}_{i}-{x}_{0}\right)}^{3},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{M}_{{xxy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}}({f}_{i}-{s}_{i}){\left({x}_{i}-{x}_{0}\right)}^{2}({y}_{i}-{y}_{0}),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{M}_{{xyy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}}({f}_{i}-{s}_{i})({x}_{i}-{x}_{0}){\left({y}_{i}-{y}_{0}\right)}^{2},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{M}_{{yyy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}}({f}_{i}-{s}_{i}){\left({y}_{i}-{y}_{0}\right)}^{3},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{M}_{{xxxx}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}^{2}}({f}_{i}-{s}_{i}){\left({x}_{i}-{x}_{0}\right)}^{4},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{M}_{{xxxy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}^{2}}({f}_{i}-{s}_{i}){\left({x}_{i}-{x}_{0}\right)}^{3}({y}_{i}-{y}_{0}),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{M}_{{xxyy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}^{2}}({f}_{i}-{s}_{i}){\left({x}_{i}-{x}_{0}\right)}^{2}{\left({y}_{i}-{y}_{0}\right)}^{2},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{M}_{{xyyy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}^{2}}({f}_{i}-{s}_{i})({y}_{i}-{y}_{0}){\left({y}_{i}-{y}_{0}\right)}^{3},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{M}_{{yyyy}}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}\displaystyle \frac{{w}_{i}}{{r}_{i}^{2}}({f}_{i}-{s}_{i}){\left({y}_{i}-{y}_{0}\right)}^{4},\end{eqnarray} \tag{ 17 }$$

where f_i is the flux in a pixel; s_i is the local sky value for that pixel; w_i is the value of the window function for the pixel; $S={\sum }_{i}({f}_{i}-{s}_{i}){w}_{i}$ is the window-weighted sum of the source flux, used to renormalize the moments; and r_i is the radius of a pixel, $\sqrt{{\left({x}_{i}-{x}_{0}\right)}^{2}+{\left({y}_{i}-{y}_{0}\right)}^{2}}$. The sums are performed over all (unmasked) pixels in the aperture. For the centroid calculation (${x}_{0},{y}_{0}$), the peak coordinate (see Section 4.4.1) is used to define the aperture and the window function; for higher-order moments, the centroid is used to center the window function.

The motivation for measuring these higher-order moments was to select exposures with image quality problems. For example, trefoil caused by errors in the collimation and alignment can in principle be detected with the third-order moments. In our experience, these statistics can be used to select some images with such problems, but we have not been able to use these values to exclude poor images from the data processing. If we were to reject images based on these moments, we would reject too many images with image quality issues that are not so poor as to preclude a useful analysis. A future machine-learning-based analysis starting with these moments might potentially provide a better rejection statistic, but such work is beyond the scope of this article.

For sources with peak flux above the saturation limit, the moments are generally poorly measured if the aperture defined by ${\sigma }_{w}$ is used. For these sources, the quality of the measurement is compromised by the saturation. However, it is still useful to estimate the first and second moments of the source in order to allow a crude measurement of the brightness from the wings of the source. In this case, a larger aperture, three times the standard aperture, is used to make a crude estimate. For such sources, the flag bit PM_SOURCE_MODE_BIG_RADIUS is set, and the source is ignored in all analyses below except for the analysis applied to very bright stars (Section 4.6.1).

If the measured centroid coordinates (${x}_{0},{y}_{0}$) differ from the peak coordinates by a large amount (1.5${\sigma }_{w}$), then the peak is identified as being of poor quality and is skipped in further analyses; the flag bit PM_SOURCE_MOMENTS_FAILURE is set for such sources. In such a case, it is likely that the "peak" was identified in a region of flat flux distribution or many saturated or edge pixels. During the analysis of the moments, the background ("sky") model is also examined for the location of each source. The value of the background and the variance of the background are recorded for each source. In some cases, the sky model or the variance is not well defined at the location of a specific sources (e.g., due to an extrapolation failure). In these cases, the flag bits PM_SOURCE_SKY_FAILURE or PM_SOURCE_SKYVAR_FAILURE are set as appropriate, and the measurement of the moments is skipped.

In addition to the moments above, the first and half-radial moments, M_r and M_h as defined below, are calculated:

$$\begin{eqnarray}&&{M}_{r}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}({f}_{i}-{s}_{i}){r}_{i},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{M}_{h}=\displaystyle \frac{1}{S}\displaystyle \sum _{i}({f}_{i}-{s}_{i})\sqrt{{r}_{i}}.\end{eqnarray} \tag{ 19 }$$

Note that the window function is not applied in the calculation of these moments.

With the first radial moment, we can calculate a preliminary Kron radius and magnitude. Kron magnitudes are provided as an option for galaxy photometry. In addition, the comparison of Kron and PSF magnitudes is useful as a star–galaxy separator. The Kron radius (Kron 1980) is defined the be 2.5× the first radial moment. The Kron flux is the sum of (sky-subtracted) pixel fluxes within the Kron radius. We also calculate the flux in two related annular apertures: the Kron inner flux is the sum of pixel values for the annulus ${R}_{1}\lt r\lt 2.5{R}_{1}$, while the Kron outer flux is the sum of pixel values for $2.5{R}_{1}\lt r\lt 4{R}_{1}$. The first radial moment is limited at the low and high ends by ${R}_{\min }\lt {M}_{r}\lt {R}_{\max }$, where ${R}_{\min }$ is the first radial moment of the PSF stars, or $0.75{\sigma }_{w}$ if that cannot be determined. ${R}_{\max }$ is set to the size of the moment's aperture, $4{\sigma }_{w}$. These Kron measurements are performed for all sources with a valid set of moments. At this stage, the measurement of the Kron parameters are preliminary as the aperture has been chosen as a fixed size relative to the size of the PSF. At a later stage, higher-quality Kron parameters appropriate to galaxies are measured with more care paid to the exact aperture used (Section 4.6.4).

4.5.1. PSF Model versus Source Model

The PSF of an image describes the shape of all unresolved sources in the image. In a typical wide-field image, the shape of unresolved sources varies as a function of position in the image. The full PSF thus needs to include a model with parameters that vary across the image.

The PSF used by psphot consists of an analytical function combined with a pixelized representation of the residual differences between the analytical model and the true PSF. Both the shape parameters of the analytical model and the pixelized residual differences are allowed to vary in two dimensions across the images.

Within psphot, several analytical models may be used to describe the smooth portion of the PSF, but all share a few common characteristics. As an example, a simple model consists of a 2D elliptical Gaussian:

$$\begin{eqnarray}&&f(x,y)={I}_{o}{e}^{-z}+S,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&z=\displaystyle \frac{{x}^{2}}{2{\sigma }_{x}^{2}}+\displaystyle \frac{{y}^{2}}{2{\sigma }_{y}^{2}}+{\sigma }_{\mathrm{xy}}{xy},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&x={x}_{\mathrm{ccd}}-{x}_{o},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&y={y}_{\mathrm{ccd}}-{y}_{o}.\end{eqnarray} \tag{ 23 }$$

Here, the model parameters consist of the centroid coordinates (${x}_{o},{y}_{o}$), the elliptical shape parameters (${\sigma }_{x},{\sigma }_{y},{\sigma }_{\mathrm{xy}}$), the model normalization (I_o ), and the local value of the background (S).

A specific source will have a particular set of values for the model parameters, some of which depend on the PSF model and the position of the source in the image, while the rest are unique to the individual source. For the case of the elliptical Gaussian model, the PSF parameters would be the shape terms (${\sigma }_{x},{\sigma }_{y},{\sigma }_{\mathrm{xy}}$) while the independent parameters would be the centroid, normalization and local sky values (${x}_{o},{y}_{o},{I}_{o},S$), though as noted below (Section 4.6.6), in practice we do not allow the sky to be fitted independently as we subtract the background model. Thus, the shape parameters are each a function of the source centroid coordinates:

$$\begin{eqnarray}&&{\sigma }_{x}={f}_{1}({x}_{\mathrm{ccd}},{y}_{\mathrm{ccd}}),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\sigma }_{y}={f}_{2}({x}_{\mathrm{ccd}},{y}_{\mathrm{ccd}}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\sigma }_{{xy}}={f}_{3}({x}_{\mathrm{ccd}},{y}_{\mathrm{ccd}}).\end{eqnarray} \tag{ 26 }$$

psphot represents the variation in the PSF parameters as a function of position in the image in two possible ways, specified by the configuration. The first option is to use a 2D polynomial that is fitted to the measured parameter values across the image. The second option is to use a grid of values that are measured for sources within a subregion of the image. In the latter case, the value at a specific coordinate in the image is determined via bilinear interpolation between the nearest grid points. The order of the polynomial or the sampling size of the grid is dynamically determined depending on the number of available of PSF stars. In the case of the PV3 analysis, the grid of values was used, with a maximum of 6 × 6 samples per GPC1 chip image (grid cells of size ∼34). For the earlier PV2 analysis, the maximum grid sampling was 3 × 3 per GPC1 chip image (grid cells of size ∼69). For the PV1 analysis, the polynomial representation was used, with up to third-order terms. The higher order representation was used for PV3 in order to follow some of the observed PSF variations in the images.

Figure 3 illustrates the 2D variations in the PSF shapes seen in PS1 data. This figure shows the FWHM and e₁ and e₂ polarizations of the stars as a function of position in four exposures. For images with good image quality, variations of the PSF shape due to the optical aberrations can be seen. The optical aberrations vary as the active collimation and alignment are adjusted and as the focus changes. These aberrations are coupled to the piston of the chips, which have been adjusted to crudely follow the focal surface (Chambers et al. 2016). During regular operations, images with large PSFs are usually caused by the atmosphere (seeing) or by telescope tracking errors, both of which result in common shapes across the field of the camera. In the figure, the top panel shows the circularization of the PSF due to the atmosphere washing out the lower-level variations caused by the optics.

Zoom In Zoom Out Reset image size

Several analytical functions that are likely candidates to describe the smooth portion of the PSF are available in psphot:

• 1.  
  Gaussian : $f={I}_{0}{e}^{-z}$
• 2.  
  Pseudo-Gaussian : $f={I}_{0}{\left(1+z+\tfrac{1}{2}{z}^{2}+\tfrac{1}{6}{z}^{3}\right)}^{-1}$ [PGAUSS]
• 3.  
  Variable Power Law : $f={I}_{0}{\left(1+z+{z}^{\alpha }\right)}^{-1}$ [RGAUSS], $\alpha \gt 1.25$
• 4.  
  Steep Power Law : $f={I}_{0}{\left(1+\kappa z+{z}^{2.25}\right)}^{-1}$ [QGAUSS]
• 5.  
  PS1 Power Law : $f={I}_{0}{\left(1+\kappa z+{z}^{1.67}\right)}^{-1}$ [PS1_V1]

The Pseudo-Gaussian is a Taylor expansion of the Gaussian and is used by Dophot (Schechter et al. 1993). The latter profiles are similar to the Moffat profile form (Moffat 1969; Buonanno et al. 1983), with small differences. For these PSF models, the functions are evaluated at the pixel center. Unlike some galaxy model representations (see Section 5.3 ), the first derivatives of these functions approach zero as the radius approaches zero, so subpixel integration is not necessary. A user may choose to try more than one analytical function for a given image. As discussed below (Section 4.5.3), psphot can automatically choose the best model based on the quality of the PSF fits.

For the PS1 GPC1 analysis, we used the PS1_V1 model, which we found by experimentation to match well to the observed profiles generated by PS1. Figure 4 shows example radial profiles for moderately bright stars in fairly good (09) and poor (22) seeing. Using a fixed power-law exponent results in somewhat faster profile fitting compared to the variable power-law exponent model.

Zoom In Zoom Out Reset image size

The analytical models in psphot are written with a high degree of code abstraction, making it relatively easy to add different analytical models to the software. The same portion of code used to describe the analytical portion of the PSF sources is also used for galaxy models.

Once the smooth component of the PSF has been fitted with an analytical model, a pixel representation of the residuals is generated. This representation is constructed as an image of the expected residuals for any position in the image. The value of each pixel in the image model is determined from 2D fits to the measured residuals of the PSF stars.

The residual model is calculated using the residuals for all PSF stars. The residuals (and their errors) for each star are renormalized by the flux of the star to put them on a consistent flux scale. For each PSF star, all pixels within a user-specified radius (PSF.RESIDUALS.RADIUS=9) are selected for the measurement. For a given pixel in the model, the value is calculated from the four closest pixels in the PSF stars via bilinear interpolation. Pixels may be used in this analysis if their S/N exceeds a user-defined limit. For the PV3 3π analysis, we allowed all pixels within the user-specified radius, not limiting on the basis of the S/N.

Pixels for a given star that are more than a number of sigmas (PSF.RESIDUALS.NSIGMA=3.0) deviant from the median value of the pixels from all stars are rejected.

If no spatial variation is allowed, the mean or median value is calculated for the model pixel based on the user-specified mean statistic (PSF.RESIDUALS.STATISTIC=ROBUST_MEDIAN).

If spatial variation is requested, then the pixel values are fitted to a linear model:

$$\begin{eqnarray*}&&\begin{array}{l}R[({x}_{\mathrm{mod}},{y}_{\mathrm{mod}})][({x}_{\mathrm{ccd}},{y}_{\mathrm{ccd}})]={R}_{o}[({x}_{\mathrm{mod}},{y}_{\mathrm{mod}})]\\ \qquad +\,{R}_{x}[({x}_{\mathrm{mod}},{y}_{\mathrm{mod}})]{x}_{\mathrm{ccd}}+\,{R}_{y}[({x}_{\mathrm{mod}},{y}_{\mathrm{mod}})]{y}_{\mathrm{ccd}},\end{array}\end{eqnarray*}$$

where $R[({x}_{\mathrm{mod}},{y}_{\mathrm{mod}})][({x}_{\mathrm{ccd}},{y}_{\mathrm{ccd}})]$ is the value of the residual for model pixel $({x}_{\mathrm{mod}},{y}_{\mathrm{mod}})$ for a star with centroid at image pixel $({x}_{\mathrm{ccd}},{y}_{\mathrm{ccd}})$. The parameters ${R}_{o},{R}_{x},{R}_{y}$ are the elements of the 2D linear fit for each pixel $({x}_{\mathrm{mod}},{y}_{\mathrm{mod}})$ in the model.

4.5.2. Candidate PSF Source Selection

The first stage of determining the PSF model for an image is to identify a collection of sources in the image which are likely to be unresolved (i.e., stars). psphot uses the source sizes as estimated from the second moments to make the initial guess at a collection of unresolved sources. At this point, the program has measured the second-order moments for all sources identified by their peaks, as well as an approximate S/N, above the bright threshold. All sources with an S/N greater than a user-defined parameter (PSF_SN_LIM = 20.0 for PS1 PV3) are selected by psphot, though sources that have more than a certain number of saturated pixels are excluded at this stage. The program then examines the 2D plane of ${M}_{x,x},{M}_{y,y}$ in search of a concentrated clump of sources (see Figure 5). To do this, it constructs an artificial image with pixels representing the value of ${M}_{x,x},{M}_{y,y}$, using $0.1{\sigma }_{w}^{2}$ as the size of a pixel in this artificial image. The binned ${M}_{x,x},{M}_{y,y}$ plane is then examined to find a significant peak. Unless the image is extremely sparse, such a peak will be well defined and should represent the sources which are all very similar in shape. Other sources in the image will tend to land in very different locations, failing to produce a single peak. To avoid detecting a peak from the unresolved cosmic rays, sources that have second moments very close to 0 are ignored. For these sources, the flag bit PM_SOURCE_MODE_DEFECT is set.

Zoom In Zoom Out Reset image size

Once a peak has been detected in this plane, the centroid and second moments of this peak are measured. All sources that land within 2 pixels of this centroid are selected as candidate PSF sources in the image.

When the second moments are measured, psphot also counts the number of saturated pixels within the analysis aperture. If more than a single saturated pixel is found, and if the second moments of that object are more than one standard deviation larger than the clump identified above, this source is identified as a highly saturated star and marked with the flag bit PM_SOURCE_MODE_SATSTAR. Sources that have more than a single saturated pixel, but for which the second moments do not exceed the above limits, are marked as likely saturated regions (e.g., bleed trails). These sources are skipped in most additional analyses and are marked with the flag bit PM_SOURCE_MODE_SATURATED.

4.5.3. Candidate PSF Source Model Fits

All candidate PSF sources are then fitted with the selected source model, allowing all of the parameters (PSF and independent) to vary in the fit. The software uses the Levenberg–Marquardt minimization technique (e.g., Press et al. 1992; Madsen et al. 2004) for the nonlinear fitting. Nonlinear fitting can be very computationally intensive, particularly if the starting parameters are far from the minimization values. The first and second moments are used to make a good guess for the centroid and shape parameters for the PSF models. Any sources that fail to converge in the fit are flagged as invalid.

To generate the initial guess, the second moments are converted to the equivalent sigma values for a 2D elliptical Gaussian contour using the following transformations inspired by Bertin & Arnouts (1996) and Stobie (1980). First, we calculate the sigma values in the major (${\sigma }_{a}$) and minor (${\sigma }_{b}$) axis directions, along with the position angle θ from the moments using

$$\begin{eqnarray}&&\theta =\displaystyle \frac{1}{2}\,\arctan 2(2{M}_{{xy}},{g}_{2}),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\sigma }_{a}=\sqrt{\displaystyle \frac{{g}_{1}+{g}_{3}}{2}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\sigma }_{b}=\sqrt{\displaystyle \frac{{g}_{1}-{g}_{3}}{2}},\end{eqnarray} \tag{ 29 }$$

where the function arctan2 (y, x) returns the arctangent in the proper quadrant (e.g,. as implemented by the atan2(y,x) function in C) and the intermediate values g₁, g₂, g₃ are given by

$$\begin{eqnarray}&&{g}_{1}={M}_{{xx}}+{M}_{{yy}},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{g}_{2}={M}_{{xx}}-{M}_{{yy}},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{g}_{3}=\sqrt{{g}_{2}^{2}+4{M}_{{xy}}^{2}}.\end{eqnarray} \tag{ 32 }$$

Because the moments may be noisy, the calculated value of ${\sigma }_{b}$ can be numerically invalid if ${g}_{3}\gt {g}_{1}$, a situation that is especially likely for highly elongated sources. We avoid this situation by limiting the axial ratio to a maximum of 20 (setting ${\sigma }_{b}$ to ${\sigma }_{a}/20$ if the expected axial ratio would be greater than this limit). The selected value of 20 is somewhat ad hoc, chosen based on failures in real images. A more careful examination of the trade-off space would be worthwhile in the future.With ${\sigma }_{a}$, ${\sigma }_{b}$, and θ in hand, we can now transform these values to the parameters of our fits, ${\sigma }_{x}$, ${\sigma }_{y}$, and ${\sigma }_{\mathrm{xy}}$ (Equation 20 above). This transformation can be determined by rotating the 2D Gaussian equation, yielding

$$\begin{eqnarray}&&{\sigma }_{x}^{-2}={\sigma }_{a}^{-2}{\rm{co}}{{\rm{s}}}^{2}\,\theta \,+\,{\sigma }_{{\rm{b}}}^{-2}{\rm{si}}{{\rm{n}}}^{2}\,\theta ,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\sigma }_{y}^{-2}={\sigma }_{b}^{-2}{\cos }^{2}\theta +{\sigma }_{a}^{-2}{\sin }^{2}\theta \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{\sigma }_{\mathrm{xy}}=\displaystyle \frac{1}{2}\sin (2\theta )({\sigma }_{b}^{-2}-{\sigma }_{a}^{-2}).\end{eqnarray} \tag{ 35 }$$

In fact, because the calculated second moments have been measured with a window function applied (see discussion in Section 4.4.3), we instead use the measured value of M_r (Equation 18), the first radial moment as the major-axis size for the Gaussian (${\sigma }_{a}$), retaining the position angle and axial ratio from the calculation above. We use these guess parameters for all versions of the PSF analytical models, despite the fact that for the versions which are not approximations of Gaussians these guess values will be systematically incorrect. It would be worthwhile in the future to tweak the guesses for the different model version to speed up the convergence.

For the resulting collection of source model parameters, the PSF-dependent parameters of the models are all fitted as a function of position using either the 2D polynomial or the gridded representation described above. The maximum order of these fits depends on the number of PSF sources (see Table 6). The fitting process for these polynomials is iterative and rejects the 3σ outliers in each of three passes. This fitting technique results in a robust measurement of the variation of the PSF model parameters as a function of position without being excessively biased by individual sources that are not well described by the PSF model (e.g., galaxies that snuck into the sample). Sources whose model parameters are rejected by this iterative fitting technique are also marked as invalid PSF sources and ignored in the later PSF model fitting stages. Sources that are actually used to define the PSF model for a given image have the flag bit PM_SOURCE_MODE_PSFSTAR set.

Table 6. Minimum Number of Stars Required for a Given Order of the PSF 2D Variations, or for the Given Number of Grid Cells

Minimum	Order	Number of	Cell Size
# of Stars		Grid Cells	(arcmin)
16	1	4	10.3
54	2	9	6.9
128	3	16	5.1
300	4	25	4.1
576	5	36	3.4

Download table as:  ASCIITypeset image

The order of the fit or number of grid samples is modified if the number of stars available for the fit is insufficient to justify the highest value. Regardless of the requested order, if the number of stars is below the following limits, the order is limited as shown in Table 6. Note that the number of grid cells in one dimension is one greater than the equivalent polynomial order.

All of the PSF candidate sources are then refitted using the PSF model to specify the PSF-dependent model parameter values for each source. For example, in the case of the elliptical Gaussian model, the shape parameters (${\sigma }_{x},{\sigma }_{y},{\sigma }_{{xy}}$) for each source are set by the coordinates of the source centroid and fixed (not allowed to vary) in the fitting procedure. The resulting fitted models are then used to determine a metric that tests the quality of the PSF model for this particular image.

psphot allows a collection of PSF model functions to be tried on all PSF candidate sources. The number of models to be tested is specified by the configuration keyword PSF_MODEL_N. The configuration variables PSF_MODEL_0, PSF_MODEL_1, through PSF_MODEL_N---1 specify the names of the models that should be tested. The metric used by psphot to assess the PSF model is the scatter in the differences between the aperture and fit magnitudes for the PSF sources. This difference is a critical parameter for any PSF modeling software as it is a measurement of how well the PSF model captures the flux of the star. Aperture photometry is measured for a circular aperture with a radius of PSF_APERTURE_SCALE (4.5 for PV3) times ${\sigma }_{w}$ (Section 4.4.3). The average aperture correction (${m}_{\mathrm{AP}}-{m}_{\mathrm{PSF}}$) is measured and, if multiple PSF model types are selected, the PSF model with the minimum clipped scatter in this statistic is chosen for the image. For the PV3 analysis, however, only the PS1_V1 model function was used.

An approximate aperture correction is measured at this stage, with a more detailed correction measured after all source analysis is performed (see Section 4.8). Sources for which the aperture magnitude is measured have the flag bit PM_SOURCE_MODE_AP_MAGS set. These aperture magnitudes are stored in the DVO field Measure.Map and exported to the PSPS as a flux in Janskies in the field Detection.apFlux. The radius (in arcseconds) of the aperture used for each exposure is reported in PSPS as Detection.apRadius, while the unmasked fraction of the aperture is reported in PSPS as Detection.apFillF.

As noted above (Section 4.3), we do not attempt to replace or interpolate masked pixel values. Aperture photometry measurements of objects that include masked pixels are thus inaccurate. For a stellar object, the amount of error is a function of how close the masked pixels are to the core of the PSF. To provide guidance, when the PSF and aperture photometry for a source are measured, two additional quantities that are useful to assess the impact of masking are measured. First, the mask image is examined and the number of unmasked pixels is summed, weighted by the normalized PSF model. The resulting quantity, PSF_QF, has a value between 0.0 (totally masked) and 1.0 (totally unmasked). Elsewhere in the IPP system, we use this value to filter out detections that are unreliable due to the masking. For a generous cut, leaning toward completeness at the cost of some lower-quality measurements, PSF_QF $\gt 0.85$ is used in some contexts; in other cases, we require PSF_QF $\gt 0.95$ to ensure a high-quality measurement (see, for example, the calculation of average photometry in Magnier et al. 2020b). The second quantity is related to the first: PSF_QF_PERFECT uses all mask values to assess the quality factor, while PSF_QF uses only the "bad" mask bit values (see Section 4.3).

Several flag bits are raised based on statistics that are similar to the PSF_QF measurement. First, psphot calculates the normalized, PSF-weighted fraction of pixels that are masked due to one of the following four causes: a diffraction spike (SPIKE), the core of a saturated star (CORE), burntool-subtracted region (BURNTOOL), or a pixel for which, due to interpolation or convolution, a significant fraction of the pixel flux comes from a masked pixel. These masking conditions are all treated as "suspect" by psphot, which means they are included in the analysis of the source pixels. However, because they may potentially affect the photometry (or astrometry), it is useful to note a source that has a nontrivial fraction of these poor mask pixels. If the normalized PSF-weighted fraction of pixels masked due to any of these four conditions exceeds 25%, then one of the following bits is raised for the corresponding condition: PM_SOURCE_MODE2_ON_SPIKE, PM_SOURCE_MODE2_ON_STARCORE, PM_SOURCE_MODE2_ON_BURNTOOL, PM_SOURCE_MODE2_ON_CONVPOOR. In addition, the following flag bits may also be raised if the central pixel of a source lands on a pixel masked for a diffraction spike (PM_SOURCE_MODE_ON_SPIKE), an optical ghost (PM_SOURCE_MODE_ON_GHOST), or off the active pixels of the CCD (PM_SOURCE_MODE_OFF_CHIP).

Once a PSF model has been determined, the brighter sources in the image may be analyzed in detail. The goals in this stage are (1) to determine the fluxes and positions of the bright stellar sources with high precision appropriate to their high S/N and (2) to characterize the bright source flux profiles sufficiently well that they may be subtracted from the image to allow for the clean detection of the fainter sources. Note that as the analysis proceeds, there are several stages in which the 2D flux models for all sources are subtracted from the image, and individual sources are replaced in the image for a particular analysis step and then removed again. The flux limit for this analysis stage is user defined as an S/N value. In the PV3 analysis of the 3π Survey data, this limit was set to an S/N of 20.0.

In order to allow for multiple threads to process a single image, the pixels in an image are divided into a grid of superpixels (note that these superpixels are not the same as those used for either the background model or the PSF parameter variations). The superpixels are assigned to one of four groups so that each superpixel in a group is well separated from the other superpixels of that group. The analysis of the image proceeds in four steps, one for each of these groups. Each of the superpixels in the first group is assigned to a single thread until all threads are assigned. A single thread is responsible for the analysis of sources that land within their current superpixel, as determined by the centroid coordinates. Because the superpixels in a given thread group are not contiguous by construction, sources near the edge of a superpixel can be analyzed by considering the nearby pixels from a neighboring superpixel (guaranteed not to be in the current thread group).

As the threads complete their analysis, they are assigned the next unfinished superpixel in the active group. When all superpixels in one group have been processed, then the superpixels in the next group can start. This strategy allows the threading to process sources that may be extended without the danger that two threads are actively touching the same pixels. For the PV3 analysis, four threads were used for most processing tasks.

4.6.1. Very Bright Stars

4.6.2. Fast Ensemble PSF Fitting

Before the detailed analysis of the sources is performed, it is convenient to subtract off all of the sources, at least as well as possible at this stage. We make the assumption that all sources are PSF like. If the centroid of the source has been determined, we use this value for its position; otherwise, we use the interpolated position of the peak. A single linear fit is used to simultaneously measure all source fluxes. Because the local sky has been subtracted, this measurement assumes the local sky is zero. We can write a single ${\chi }^{2}$ equation for this image:

$$\begin{eqnarray*}&&{\chi }^{2}=\displaystyle \sum _{\mathrm{pixels}}{\left({F}_{x,y}-\displaystyle \sum _{\mathrm{sources}}{A}_{i}P[{x}_{0},{y}_{0}]\right)}^{2},\end{eqnarray*}$$

where ${F}_{x,y}$ is image flux for each pixel, $P[{x}_{0},{y}_{0}]$ is the PSF model realized at the position of source i, and A_i is the normalization for the source.

Minimizing this equation with respect to each of the A_i values results in a matrix equation:

$$\begin{eqnarray*}&&{M}_{i,j}\bar{{A}_{i}}\,=\,\bar{{F}_{j}},\end{eqnarray*}$$

where $\bar{{A}_{i}}$ is the vector of A_i values, the elements of ${M}_{i,j}$ consist of the dot products of the unit-flux PSF for source i and source j, and $\bar{{F}_{j}}$ is the dot product of the unit-flux PSF for source j with the pixels corresponding to source j. The dot products are calculated only using pixels within the source apertures. Because most sources have no overlap with most other sources, this matrix equation results in a very sparse, mostly diagonal square matrix. The dimension is the number of sources, likely to be 1000 s or 10,000 s. Direct inversion of the matrix would be computationally very slow. However, an iterative solution quickly yields a result with sufficient accuracy. In the iterative solution, a guess at the solution $\bar{A}$ is made assuming ${M}_{i,j}$ is purely diagonal; the guess is multiplied by ${M}_{i,j}$, and the result is compared with the observed vector $\bar{{F}_{j}}$. The difference is used to modify the initial guess. This process is repeated several times to achieve convergence. Convergence is quick (a few iterations) because of the highly diagonal matrix with small off-diagonal terms: the dot product of source i and source j is 1, where i = j and much less than 1 where $i\ne j$.

Once a solution set for A_i is found, all of the sources are subtracted from the image by applying these values to the unit-flux PSF. Sources for which a PSF model has been fitted (whether or not this is retained as the best model in the end) has the flag field PM_SOURCE_MODE_PSFMODEL set. All sources that are included in the ensemble linear fit have the flag bit PM_SOURCE_MODE_LINEAR_FIT set, including those for which the model is not the PSF.

4.6.3. Radial Profile Wings

4.6.4. Kron Magnitudes

4.6.5. Source-size Assessment

After the PSF model has been fitted to all sources and the Kron flux has been measured for all sources, psphot uses these two measurements, along with some additional pixel-level analysis, for classification based on source size. Sources identified as extended will be fitted with a galaxy model (or possibly another type of extended source model in special cases). If the source is small compared to a PSF, it is considered to be a cosmic ray and masked.

Extended sources are identified as those for which the Kron magnitude is significantly brighter than the PSF magnitude when compared to a PSF star. The value $\delta {M}_{\mathrm{KP}}\,=\,{m}_{\mathrm{Kron}}\,-{m}_{\mathrm{PSF}}$, the difference between the PSF and Kron magnitudes, is calculated for each source. The median of $\delta {M}_{\mathrm{KP}}$ is calculated for the PSF stars. This median is subtracted from $\delta {M}_{\mathrm{KP}}$ for each star. The result is divided by the quadrature error of the PSF and Kron magnitudes and called extNsigma. If extNsigma is larger than the configuration value PSPHOT.EXT.NSIGMA.LIMIT (3.0 for PV3), the source is considered to be extended, and the flag bit PM_SOURCE_MODE_EXT_LIMIT is set for the source.

We decided to use the $\delta {M}_{\mathrm{KP}}$ metric for this assessment after we tested several possible star–galaxy separation statistics. We found that the Kron PSF comparison was more reliable than second-moment and first-radial-moment based measurements. In addition, because we needed a statistic that could be calculated relatively quickly on every detected source, we rejected using a galaxy model fit for the star–galaxy separator.

Cosmic rays are identified by a combination of the Kron magnitude and the second-moment width of the source in the minor axis direction. The second moment in the minor axis direction is calculated from ${M}_{{xx}},{M}_{{xy}},{M}_{{yy}}$ as follows:

$$\begin{eqnarray*}&&{M}_{\mathrm{minor}}\,=\,\displaystyle \frac{1}{2}({M}_{{xx}}+{M}_{{yy}})-\displaystyle \frac{1}{2}\sqrt{{\left({M}_{{xx}}-{M}_{{yy}}\right)}^{2}+4{M}_{{xy}}^{2}}.\end{eqnarray*}$$

If ${M}_{\mathrm{minor}}\lt 0.8$ pixels² and the S/N of the flux measured in the moments analysis $\gt 7$, then the source is identified as a cosmic ray and the associated pixels are masked. These values are tuned empirically for the PV3 analysis based on cosmic rays identified in the GPC1 images. Sources that are determined to be a cosmic ray in this manner have the flag bit PM_SOURCE_MODE_DEFECT set.

The pixels of any suspected cosmic ray identified above are examined in additional detail to make a final judgment. The Laplacian edge detection algorithm based on van Dokkum (2001) is used to check for sharp edges in the flux distribution. If the sharpness exceeds a defined limit, then the pixels are masked and the flag bit PM_SOURCE_MODE_CR_LIMIT is set for the source.

4.6.6. Full PSF Model Fitting

4.6.7. Double and Blended Sources

4.6.8. Non-PSF Sources

After a first pass through the image, in which the brighter sources above a high threshold level have been detected, measured, and subtracted, psphot optionally begins a second pass at the image. Some of the steps described in the previous sections are repeated in this analysis, though with some modifications as discussed below.

The source detection steps described in Section 4.4 are repeated. To start, the sources detected in the previous steps are subtracted from the image, and the variance enhanced by adding the variance predicted by the model to the variance image, doubling the effective variance at the location of previously detected sources. As in Section 4.4.1, the image is smoothed, but in this pass, it is convolved with the PSF determined above, not a placeholder Gaussian. The new peaks are detected on the smoothed image. The peak-detection process again uses the variance image to test the significance of the individual peaks. All peaks with a significance greater than a user-defined minimum threshold are accepted as sources of potential interest. Footprints are again generated as in Section 4.4.2.

Next, moments are measured as in Section 4.4.3. In this pass, however, the size of the window function applied for the measurement of the moments is fixed at the value determined from the bright source analysis. All sources, including those measured in the bright-source analysis (which are readded to the image and their variance reset), are then simultaneously fit for their flux normalizations as in Section 4.6.2. In this step, the "best" model is used for each source, either a PSF model or the unconvolved extended source model determined in Section 4.6.8. For the newly detected sources, the PSF model is used, with the position set by the centroids.

After the flux normalization is calculated, the radial profile is measured (Section 4.6.3) and the moments are used to calculate the preliminary Kron radius and flux (see Section 4.6.4). These are in turn used to assess the source sizes as in Section 4.6.5. However, the nonlinear fitting steps for the PSF model fits (Section 4.6.6), and the extended source model fits (Section 4.6.8) are not performed for these faint sources. These steps are skipped for two reasons: first, the nonlinear fitting steps are costly in terms of computation time and the faint sources usually far outnumber the brighter sources. Second, because these are faint sources, they do not have the S/N to constrain models with many additional parameters. In addition, the positions (for PSF sources) are not much improved using the nonlinear fitting compared with the nonparametric centroid measurement for these faint sources.

The PV3 threshold for the bright source analysis is an S/N of 20. The flag bit PM_SOURCE_MODE2_PASS1_SRC is raised for sources detected in this initial analysis stage. The lower limit cutoff for the faint source analysis in PV3 is an S/N of 5.0.

In the psphotStack version of the code, the five filter images are processed together. In this case, any source that is detected in at least two of the five filters are then also measured on the other filter images in which it was not detected above the S/N limit. The position in the other stack images is fixed based on the pixel coordinates in the images in which the source was detected. Detection in two filters is required in order to avoid excessive forced photometry of spurious detections. There is an interesting class of astronomical objects that are extremely red (e.g., brown dwarfs and high-redshift quasars). Such sources are expected to be detected only in the reddest filter (${y}_{{\rm{P}}1}$). For the 3π PV3 processing, we therefore also force the photometry in all filters for sources that are only detected in ${y}_{{\rm{P}}1}$. All sources that are forced on one image based on detections in other images have the flag bit PM_SOURCE_MODE2_MATCHED set.

A PSF model will always fail to describe the flux of the stellar sources at some level. For high-precision photometry, we need to be able to correct for the difference between the PSF model fluxes and the total flux of the sources. In the end, all astronomical photometry is in some sense a relative measurement between two images. Whether the goal is calibration of a science image taken at one location to a standard star image at another location, or the goal is simply the repetitive photometry of the same star at the same location in the image, it is always necessary to compare the photometry between two images. If this measurement is to be consistent, then the measurement must represent the flux of the stars in the same way regardless of the conditions under which the images were taken, at least within some range of normal image conditions. So, for example, two images with different image quality, or with different tracking and focus errors, will have different PSF models. To the extent the PSF model is inaccurate, the measured flux of the same source in the two images will be different (even assuming all other atmospheric and instrumental effects have been corrected). The amplitude of the error will by determined by how inconsistently the models represent the actual source flux.

Aperture photometry attempts to avoid these problems but introduces other difficulties. In aperture photometry, if a large-enough aperture is chosen, the amount of flux that is lost will be a small fraction of the total source flux. Even more importantly, as the image conditions change, the amount lost will change by an even smaller fraction, at least for a large aperture. Aperture photometry can then be used to correct the PSF photometry.

The difficulty for aperture photometry is the need to make an accurate measurement of the local background for each source. As the aperture grows, errors in the measurement of the sky flux start to become dominant. If the aperture is too small, then variations in the image quality are dominant. The brighter the source, the smaller the error introduced by the large size of the aperture. However, the number of very bright stars is limited in any image, and of course, the brighter stars are more likely to suffer from nonlinearity or saturation.

In order to thread the needle between these effects, psphot measures the aperture photometry on a modest-sized aperture and then uses the PSF model to extrapolate to a large aperture. When the PSF fluxes are calculated, the aperture flux for the modest-sized aperture is also determined. The aperture is a circular aperture with radius set to a fixed multiple (PSF_APERTURE_SCALE) of ${\sigma }_{w}$, the width of the Gaussian window function determined based on the analysis of the second moments (see Section 4.4.3). For the PV3 3π analysis, the aperture window radius is $4.5\times {\sigma }_{w}$, while the large reference aperture radius is set to 25 pixels (PSF_REF_RADIUS = 64). These corrected aperture magnitudes are saved in the output catalogs as AP_MAG, the uncorrected aperture magnitudes are saved as AP_MAG_RAW, and the radius used to measure the raw aperture flux is saved as AP_MAG_RADIUS. The corresponding flux and the flux error are saved as AP_FLUX and AP_FLUX_SIG.

With these aperture magnitudes in hand, it is now possible to make an average correction to the PSF magnitudes to bring the PSF and aperture magnitudes to the same system. This correction is measured using the same stars from which the PSF model is measured, as long as the PSF magnitude error for the star is less than 0.03 mag. The correction is calculated using the weighted average of the values ${m}_{\mathrm{AP}}-{m}_{\mathrm{PSF}}$. Because the PSF may vary across the image, the correction is determined as a function of position in the image. Like the PSF model, the 2D variations of the aperture correction may be modeled as a polynomial or via interpolation in a grid. For the 3π PV3 analysis, a grid with a maximum of 6 × 6 samples per GPC1 chip image was used. The reported PSF magnitudes for all objects have this aperture correction applied. Note that an initial aperture correction was measured during the initial steps of the analysis before the PSF model was chosen. However, because the sources in the image were not yet measured and subtracted, that aperture could be contaminated by neighbors. The analysis here is performed one fairly bright star at a time with all other sources subtracted in order to minimize such contamination.

At the end of the psphot analysis of the sources in the image, an analysis is performed to test the detection efficiency. A number of fake PSF sources are injected into the image and the peak-detection analysis (Section 4.4.1) is used to determine if these sources would have been recovered. The PSF model fluxes are measured for the source that are detected. For a given image, the detection threshold is predicted based on the median image variance and the seeing. A series of brightness bins straddling the threshold are defined and a number of sources are injected with magnitudes corresponding to each of these bin values. The psphot recipe value EFF.NUM specifies the number of sources in each brightness bin (500 the PV3), and the value @EFF.MAG specifies the bins as magnitudes above and below the threshold. For PV3, the 13 mag offsets were (−2.0, −1.0, −0.5, −0.25, −0.1, −0.05, 0.0, 0.05, 0.1, 0.25, 0.5, 1.0, 2.0), providing fine sampling near the limit, but coarser coverage farther away. Poisson noise appropriate to the photon counts of the injected sources is included in the image. Injected sources are uniformly distributed across the image in X and Y pixel coordinates without any consideration of the masked regions. This last point means the recovered fraction in the bright bins can be used to test the masking fraction.

As the stellar density increases, the completeness suffers due to crowding and confusion. Because the injection and recovery analysis of the fake sources operates on the source-subtracted image and does not attempt to fully discover the sources, this analysis overestimates the completeness in crowded fields. To explore the completeness in crowded field images, we generate a series of simulated images using a Gaussian PSF with FWHM = 1'' for a range of stellar densities. We generate fake stars with fluxes as faint as one-fifth of the flux as the low-density detection limit, with densities ranging from ∼14,000 stars per square degree at low-density detection limit to ∼4.8 million stars per square degree at the low-density detection limit. The latter is comparable to observed densities in the Galactic plane. We run the psphot analysis on these simulated images and compare the detected stars to those injected to calculate the completeness for each image as a function of the true magnitude of the stars. Figure 6 shows the measured completeness for each of the six simulated images, labeled by the logarithm of their faint-end stellar density. The red dashed line shows the expected detection limit based on the background and seeing, while the red curve shows the completeness curve calculated automatically by psphot using the injection and recovery analysis.

Zoom In Zoom Out Reset image size

For low-density fields, the completeness function determined by injection and recovery is similar to that measured by the simulation, with the 50% completeness threshold roughly 0.3 mag too faint. As the stellar density increases, the true 50% completeness magnitude rises relative to the value estimated by injection and recovery.

Ideally, all sources detected by psphot would correspond to real astrophysical objects. In reality, many sources that do not correspond to real sources in the sky are detected in the images. In the very simplified simulations discussed above, which do not include realistic detector artifacts, we find that the fraction of bogus detections is extremely low, even at the very faint end. In real data, bogus detections are due to a variety of typical instrumental features including cosmic rays, diffraction spikes, satellite tracks, glows, non-Gaussian noise, variance misestimation, etc. See Paper III for an extensive discussion of instrumental artifacts in the Pan-STARRS images.

Figure 7 illustrates the completeness and bogus detection fraction for a set of four real PS1 exposures from the 3π Survey. This figure uses ${i}_{{\rm{P}}1}$-band exposures with Galactic longitude roughly 200° and latitudes of 0°, 10°, 30°, and 90°. We identify the real astrophysical sources in these fields by comparing with the deeper stack exposures and counting as real any source detected in both ${r}_{{\rm{P}}1}$ and ${i}_{{\rm{P}}1}$. We correct for the masking fraction in the exposures (which is roughly 80%) in the case of GPC1 and plot the completeness fraction for all detections in 0.5 mag wide bins from the saturation limit to below the detection limit. We also show the bogus fraction, calculated as $1-{f}_{\mathrm{pure}}$, where ${f}_{\mathrm{pure}}$ is the ratio of real detections to all detections for the given sample. We then apply three cuts to remove certain kinds of bogus sources. First, we exclude cosmic rays identified by psphot by rejecting sources with the flag bit PM_SOURCE_MODE_CR_LIMIT (see Section 4.6.5). Next, we also remove detections with PSF_QF less than 0.95. Because this cut removes detections with heavy masking, it exclude a number of bogus detections due to glows and edge defects. Finally, we also exclude detections with PSF_QF_PERFECT less than 0.95. This cut removes detections on residual persistent glows and diffraction spikes.

Zoom In Zoom Out Reset image size

For the exposures at high Galactic latitude, with a relatively low density of sources, the cosmic rays represent a significant contamination, as seen in the excess of bogus sources with ${i}_{{\rm{P}}1}$-band magnitudes in the range 17–19. These are efficiently removed with the cosmic-ray cut listed above without noticeable impact on the completeness. The other two cuts remove significant numbers of bogus detections, especially at the faint end, but at a significant cost in completeness at even brighter magnitudes. The completeness impact of these cuts is more significant at low Galactic latitude, likely because the chance of having a source lie on the diffraction spikes or persistence glows is greatly increased at higher stellar densities. The impact of the crowding on the completeness is also clear in this data set.

To illustrate the quality of the stellar photometry as measured with PSF and aperture magnitudes, we show the results of an analysis of a set of 18 images obtained by PS1 on 19 February 2010. These images were obtained for the stellar transit survey "Pan-Planets" (Obermeier et al. 2016) and thus target a relatively dense Galactic plane field. The observations were obtained with approximately consistent pointing, reducing our sensitivity to small-scale variations in the flat-field structures.

Figures 8 and 9 show comparisons of the PSF and aperture photometry measured for these 18 images. In these figures, the photometry has been measured using the configuration for psphot as used for the full PV3 chip analysis. The first image of the sequence is compared to the remaining 17 images. A relative zero-point correction is applied, measured as the median of the photometry difference for stars with S/N greater than 50. The combined error is reported and used to generate the histograms shown in the figures. From these two figures, one can observe the trade-off between PSF and aperture photometry. For the brightest instrumental magnitudes, corresponding to S/Ns greater than roughly 300, aperture magnitudes provide a more consistent measurement of the stellar fluxes, while the PSF magnitudes are more reliable for fainter sources. Catastrophic failures or extreme outliers are also reduced for the PSF photometry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We largely attribute the behavior on the bright end to systematic errors in the photometry due to our inability to perfectly represent the shape of the PSF. The PSF of stars at the bright end will depend on the brightness because of the "Brighter-Fatter" effect (Antilogus et al. 2014; Gruen et al. 2015), in which the charge already present in the pixels will force the newly arriving photoelectrons to be systematically pushed away from the accumulating stellar image, but we do not include a brightness term in our PSF model. Detector or electronic nonlinearity may also affect the PSF shape and thus the PSF photometry, though nonlinearity will affect the reported photometry for both PSF and aperture magnitudes.

We believe the observed behavior at the faint end is primarily a result of the increased crowding. Aperture photometry is more adversely affected by close neighbors than PSF photometry. Compared to the formal errors, the faint PSF photometry is the most reliable, with the aperture photometry degrading rapidly as the flux of the star decreases.

The figures above show the relative photometric accuracy for observations at a consistent pointing compared to the photon-counting statistics. A related question is to ask how consistent is the photometry of the very brightest stars in terms of magnitudes. Figure 10 shows the accuracy of the brightest stars in these images for both PSF and aperture magnitudes. The relative zero-point between the first image in the sequence and each of the remaining images was calculated and the standard deviations were measured using stars 7 to 8 mag brighter than the detection threshold, for which the photon noise is less than 1 mmag. Significant zero-point differences between the images are observed, largely due to the atmospheric transparency variations. Even so, the relative zero-points calculated from the aperture magnitudes have standard deviations of 2.4–7.4 mmag with a median of 3.5 mmag, while for PSF magnitudes, the standard deviations are in the range 6.7–14.2 mmag, with a median of 9.2.

Zoom In Zoom Out Reset image size

Our ultimate ability to accurately measure the brightness of individual sources depends on a few factors: the accuracy of the flat-field response, the consistency of the flux measurement across the image (either due to the accuracy of the PSF model or the accuracy of the aperture correction), and the accuracy of our correction for any zero-point changes. Our ability to accurately measure the zero-point of each exposure depends in part on the characteristics of the observing site. In hazy conditions, the transparency of the atmosphere may vary substantially in time but be relatively stable across the field of view of the camera, as is shown in Figure 10. Conversely, thin patchy clouds can result in small average transparency changes but substantial localized variations. If the site experiences more patchy clouds than smooth haze, photometric calibration will be difficult. A large fraction of time with cloudless conditions will benefit the calibration.

To examine the Pan-STARRS site characteristics, we extracted ${i}_{{\rm{P}}1}$ zero-points for the lifetime of the observatory (2009 June–2020 April), shown in Figure 11. These zero-points were measured as part of the PV3 analysis of the 3π Survey, and from the nightly data analysis after the end of the 3π Survey, in both cases using the Pan-STARRS-based reference catalog. The zero-points vary from night to night and over long periods. Over the 11 years of PS1 operation, the observed ${i}_{{\rm{P}}1}$-band zero-point (for data in good weather, extrapolated to the zenith) has varied over 0.175 mag (see Figure 11). The long-term variations are believed to be due mostly to dust accumulation on the primary mirror and occasional cleaning, though the effect of the atmosphere cannot be ruled out.

Zoom In Zoom Out Reset image size

Figure 12 shows a log-scale histogram of the ${i}_{{\rm{P}}1}$-band zero-points after subtracting a smoothly varying spline fit to the median of groups of 10 nights. A Gaussian fit to this distribution has $\sigma \,=\,26.6$ mmag. If we alternatively subtract a median zero-point for each night, the standard deviation is reduced to 17.6 mmag. These values can be compared to the results of Schlafly et al. (2012), in which only photometric nights were included, yielding a standard deviation of 9.0 mmag. On short timescales, weather (e.g., clouds & haze) causes the deviations to lower zero-point values. A small fraction of positive deviations also seen in Figure 12, which are not expected from the normal effects of weather. We believe these are largely due to aperture correction errors.

Zoom In Zoom Out Reset image size

This section is focused on the basic analysis of the image for point-source detection and measurement. This analysis is applied as described to the individual exposures in the chip-stage analysis and the measurements are exposed in the public release PSPS database in the Detection table. The same analysis is applied to the individual skycells in the stack-stage analysis and the resulting values are presented in the PSPS StackObjectThin and StackObjectAttribute tables, with the latter presenting values in instrumental units and the former giving calibrated values. The detection efficiency information determined from the injection and recovery analysis is stored in the ImageDetEffMeta and StackDetEffMeta tables for the chip- and stack-stage analysis.

Galaxies with regular profiles, such as elliptical galaxies and regular spiral galaxies, may be described as primarily a radial surface brightness profile, with additional structure acting as a perturbation on that profile. For many galaxies, the azimuthal shape at a given isophotal level may be described as an elliptical contour. To first order, a galaxy may be well described with a single elliptical contour and radial profile.

In order to facilitate the Petrosian photometry analysis below, psphot generates a radial profile for each suspected galaxy. This analysis starts by generating a radial profile in 24 azimuthal segments. Near the center of the galaxy, the profile is defined for radial coordinates in steps of 1 pixel, with the closest pixel values interpolated to that radial position. Further from the center, the profile is defined using the median of the pixels landing in an annular segment of size $\delta R=r\sin \theta$, rounded up to the nearest integer pixel value. The median of all pixels within a rectangular approximation to the radial wedge is used.

The resulting 24 radial profiles are subject to contamination from neighboring sources or to NAN values from masked pixels. To clean the profiles, pairs of radial profiles from opposite sides of the source are compared. Any masked values are replaced by the corresponding value in the other profile. The minimum of both profiles is then kept for both profiles. The result of this analysis is a set of profiles of the form ${f}_{i}({r}_{i})$. In this case, f_i is effectively the surface brightness for each radius in instrumental counts per pixel. If fewer than four radial surface brightness values are available for the analysis, the source is skipped and the flag bit PM_SOURCE_MODE2_ECONTOUR_FEW_PTS is set. Some apparently extended sources are in fact bright stars with a central saturation. These sources show up in this analysis as having many NAN-valued pixels in the central regions. During the radial profile analysis, such sources are flagged with the bit PM_SOURCE_MODE2_RADBIN_NAN_CENTER and are skipped from the rest of the analysis.

The surface brightness profiles are then used to define the azimuthal contour at a specific isophotal level. This contour will be used to rescale the radial profiles into a single set of profiles normalized by the elliptical contour. This contour is defined by determining the median radius for profile bins with surface brightness in the range ${F}_{\min }+0.1{F}_{\mathrm{range}}$ to ${F}_{\min }\,+0.5{F}_{\mathrm{range}}$. The result of this analysis is a value for the radius as a function of the angle for a well-defined surface brightness regime. We then determine the elliptical shape parameters for this elliptical contour: ${R}_{\mathrm{major}},{R}_{\mathrm{minor}},\theta$. This ellipse is then used to redefine a single radial profile normalized by the elliptical contour:

$$\begin{eqnarray*}&&\rho =\sqrt{\displaystyle \frac{{x}^{2}}{{S}_{{xx}}^{2}}+\displaystyle \frac{{y}^{2}}{{S}_{{yy}}^{2}}+{{xyS}}_{{xy}}}.\end{eqnarray*}$$

The surface brightness values are sampled at a number of radial annuli, with the radii defined in the configuration (RADIAL.ANNULAR.BINS.LOWER & RADIAL.ANNULAR.BINS.UPPER). For each source, the resulting surface brightness profile is saved in the output FITS table as a vector (PROF_SB). The flux at each radial position and the fill factor (fraction of pixels used to the total possible) are also saved as equal-length vectors in the FITS table (PROF_FLUX and PROF_FILL). The values of the radial bins are saved in the output file FITS header (RMIN_NN, RMAX_NN). These measurements are saved in the catalog FITS files generated by psphot, but they are not currently exported to the PSPS database for easy access.

Petrosian (1976) defined an adaptive aperture based on a ratio of surface brightnesses. The motivation is to define an aperture that can be determined for galaxies without significant biases as a function of distance from the observer. As the surface brightness profile in a resolved source is conserved as a function of distance, using a ratio of surface brightness to define a spatial scale results in a spatial scale that is constant regardless of galaxy distance.

To measure the Petrosian radius and flux, we start by defining a series of radial apertures with power-law spacing: ${r}_{i+1}\,=\,1.25{r}_{i}$. We calculate the surface brightness for the annulus from ${r}_{i}-{r}_{i+1}$ by calculating the median of the values in the range ${r}_{i}/\sqrt{1.25}$ to ${r}_{i+1}\sqrt{1.25}$ and dividing the effective area of the annulus corresponding to ${r}_{i}-{r}_{i+1}$.

For any annulus i spanning the radii ${r}_{\min }$ to ${r}_{\max }=\beta {r}_{\min }$, the Petrosian ratio for that annulus is defined as the ratio of the surface brightness in the annulus to the average surface brightness within ${r}_{\max }$. The Petrosian radius is defined to be r_max for the annulus for which the Petrosian ratio = 0.2, i.e., the point on the galaxy radial profile at which the surface brightness is 20% of the average surface brightness at that point. If the profile falls below the Petrosian ratio for the first radial bin, the flag bit PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN is set to note that the Petrosian radius may be poorly determined.

We determine the Petrosian radius for the galaxy by quadratic interpolation between the last two of the fixed annuli with Petrosian ratio $\gt 0.2$ and the first annulus with Petrosian ratio $\lt 0.2$. In general, the Petrosian ratio for a galaxy with a regular morphology (spiral or elliptical) is falling monotonically, so this interpolation is unambiguous. However, irregular galaxy morphologies, noise, and/or significant masking can cause the Petrosian ratio to have rises as well as drops. We track the Petrosian ratio until the value is no longer significant (${\sigma }_{\mathrm{Ratio}}\lt 2\mathrm{Ratio}$). If the Petrosian ratio drops below 0.2 for more than one radius, we choose the largest such radius. If the Petrosian ratio does not fall below 0.2 for any of the measured radii, we choose the annulus for which the ratio falls to the lowest (yet still significant) value. In such a case, the flag bit PM_SOURCE_MODE2_PETRO_INSIG_RATIO is set.

Once the Petrosian radius has been determined, we can now measure the Petrosian flux: this is defined to be the total flux within an aperture corresponding to 2× the Petrosian radius. Using the Petrosian flux, we can calculate two other interesting radii: R₅₀ and R₉₀, the radii inside which 50% and 90% of the total Petrosian flux is contained. Sources for which the Petrosian parameters are successfully measured have the flag bit PM_SOURCE_MODE_EXTENDED_STATS set. Sources for which the Petrosian parameters were attempted but for which the radial profile analysis failed have the flag bit PM_SOURCE_MODE2_PETRO_NO_PROFILE set. These measurements are available from the PSPS StackPetrosian table.

Our implementation of the Petrosian apertures and fluxes is designed to match the SDSS implementation (Stoughton et al. 2002), and therefore, the measured parameters should be quite comparable between the two surveys. Figure 14 compare the Petrosian magnitudes and radii as measured by psphot on the 3π Survey observations and the values measured by SDSS for the same objects. Objects identified by SDSS as galaxies (probPSF_r $\lt 0.5$) near the Galactic north pole (α = 180° to 190°s, δ = 25° to 35°s) are selected from the PS1 3π Survey data set base on positional coincidence. The figure shows the difference in the r-band Petrosian magnitudes as a function of the Petrosian magnitude and as a function of the difference in the measured Petrosian radii. Differences in the measured magnitudes are driven by differences in the size estimates from the two data sets and analysis methods. The PS1 analysis tends to find larger radii than the SDSS analysis for the same objects, with a mean difference of 03. The larger aperture results in more flux captured in the aperture and thus brighter magnitudes for the same object: the mean difference is –0.23 magnitude in the sense of larger fluxes for the PS1 measurements.

Zoom In Zoom Out Reset image size

In the galaxy model-fitting stage, sources that meet certain criteria are fitted with analytical models for galaxies. Three traditional analytical galaxy models are implemented in psphot and used in the PV3 analysis:

• 1.  
  Exponential profile : $f={I}_{0}{e}^{-\rho },$
• 2.  
  de Vaucouleur (1948) profile: $f={I}_{0}{e}^{-{\rho }^{1/4}},$
• 3.  
  Sérsic (1963) : $f={I}_{0}{e}^{-{\rho }^{1/n}},$

where ρ is a normalized radial term: $\rho =\sqrt{\tfrac{{x}^{2}}{{R}_{{xx}}^{2}}+\tfrac{{y}^{2}}{{R}_{{yy}}^{2}}+{{xyR}}_{{xy}}}$. The terms (R_xx , R_yy , R_xy ) describe the elliptical contour and the profile scale in all three models and the coordinates x and y are determined relative to the centroids ($x,y={X}_{\mathrm{chip}}\,-{x}_{0},{Y}_{\mathrm{chip}}-{y}_{0}$). Including the normalization (I₀) and a local sky value, the Exponential and de Vaucouleur profiles have seven free parameters and the Sérsic profile has the additional free parameter of the Sérsic index n. In this stage, the galaxy model is convolved with an approximation to our best guess for the PSF model at the location of the galaxy.

Sources that passed the extended source restrictions described above were fitted with all three galaxy models, unless (a) the morphological test identified the source as a likely cosmic ray (Section 4.6.5) or (b) the peak of the PSF profile was above the saturation limit for the chip (see the discussion in Waters et al. 2020, regarding the masking of saturated pixels). For all sources for which the extended source model fits were attempted, the flag bit PM_SOURCE_MODE2_EXT_FITS_RUN is set. If any of the attempted model fits failed, then the flag bit PM_SOURCE_MODE2_EXT_FITS_FAIL is set. If all model fits failed, then the flag bit PM_SOURCE_MODE2_EXT_FITS_NONE is set.

Before the nonlinear fitting may be performed, it is necessary to determine initial values for the parameters to be fitted. For each of the three model types, the position determined from the PSF fitting analysis is used as the initial centroid ${x}_{0},{y}_{0}$. A guess for the terms (R_xx , R_yy , R_xy ) is generated based on the second moments. The guess does not attempt to use the PSF model to adjust the (R_xx , R_yy , R_xy ) values; it was found that such a guess tended to be too small and resulted in more iterations rather than fewer. The first radial moment (see Section 4.4.3) is used to estimate the effective radius of the model based on the results of (Graham & Driver 2005, Table 1). They quantify the relationships between the first radial moment used to calculated a Kron magnitude and the effective radius for different Sérsic index values, n. Because the Exponential and de Vaucouleur models are equivalent to Sérsic models with n = 1 and 4, respectively, this work can be used to generate the initial effective radius values for all three model types. Once the effective radius is chosen, the second moments are used to define the aspect ratio and position angle of the elliptical contour, as described for PSF sources in Section 4.5.3. The Kron flux is used to generate a guess for the normalization, applying an appropriate scale factor based on the (R_xx , R_yy , R_xy ) values, generated by integrating normalized Sérsic models and determining the relationship between the central intensity and the integrated flux as a function of the Sérsic index.

The PSF-convolved galaxy model fitting analysis uses the Levenberg–Marquardt minimization method to determine the best fit. In this process, the ${\chi }^{2}$ value to be minimized is

$$\begin{eqnarray*}&&{\chi }^{2}(\bar{a})=\displaystyle \sum _{p}\displaystyle \frac{1}{{\sigma }_{p}^{2}}{\left[{I}_{p}-{M}_{p}(\bar{a})\otimes \mathrm{PSF}\right]}^{2},\end{eqnarray*}$$

where I_p represents the pixel values in the image (within some aperture) and ${M}_{p}(\bar{a})$ represents the unconvolved galaxy model, a function of a number of parameters $\bar{a}$, which is then convolved with the PSF model.

To determine the minimization, we need the gradient and Laplacian of ${\chi }^{2}$ with respect to the model parameters, a_m :

$$\begin{eqnarray}&&{\chi }^{2}(\bar{a})=\displaystyle \sum _{p}{f}_{p}^{2},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&2{\rm{\nabla }}{\chi }^{2}=\displaystyle \sum _{p}{f}_{p}\displaystyle \frac{\partial {f}_{p}}{\partial {a}_{m}},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\chi }^{2}\sim {H}_{m,n},\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&2{H}_{m,n}=\displaystyle \sum _{p}\displaystyle \frac{\partial {f}_{p}}{\partial {a}_{m}}\displaystyle \frac{\partial {f}_{p}}{\partial {a}_{n}},\end{eqnarray} \tag{ 39 }$$

where we define

$$\begin{eqnarray}&&{f}_{p}({a}_{m})=\displaystyle \frac{1}{{\sigma }_{p}}({I}_{p}-{M}_{p}\otimes \mathrm{PSF}),\end{eqnarray} \tag{ 40 }$$

and we have approximated the Laplacian with the Hessian matrix, ${H}_{m,n}$, by dropping the second derivatives (which are assumed to be a small perturbation).

Because

$$\begin{eqnarray*}&&\displaystyle \frac{\partial {f}_{p}}{\partial {a}_{m}}\,=\,-\displaystyle \frac{1}{{\sigma }_{p}}\displaystyle \frac{\partial {M}_{p}\otimes \mathrm{PSF}}{\partial {a}_{m}}\end{eqnarray*}$$

and because the order of the derivative and convolution may be exchanged, we can write these in terms of the convolved image of the model and the convolved images of the derivatives of the model M_p with respect to the model parameters, a_m :

$$\begin{eqnarray}&&{{ \mathcal M }}_{p}={M}_{p}\otimes \mathrm{PSF},\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{{ \mathcal M }}_{p,m}^{{\prime} }=\displaystyle \frac{\partial {M}_{p}}{\partial {a}_{m}}\otimes \mathrm{PSF},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&2{\rm{\nabla }}{\chi }^{2}=-\displaystyle \sum _{p}\displaystyle \frac{{I}_{p}-{{ \mathcal M }}_{p}}{{\sigma }_{p}}{{ \mathcal M }}_{p,m}^{{\prime} },\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&2{H}_{m,n}=\displaystyle \sum _{p}\displaystyle \frac{1}{{\sigma }_{p}^{2}}{{ \mathcal M }}_{p,m}^{{\prime} }{{ \mathcal M }}_{p,n}^{{\prime} }.\end{eqnarray} \tag{ 44 }$$

The gradient vector and Hessian matrix are used in the Levenberg–Marquardt minimization analysis using the standard technique of determining a step from the current set of model parameters to a new set by solving the matrix equation:

$$\begin{eqnarray*}&&(1+{\lambda }_{m,n}){H}_{m,n}\,=\,\delta {\rm{\nabla }}{\chi }^{2},\end{eqnarray*}$$

where ${\lambda }_{m,n}$ is zero for $m\ne n$ and for m = n set to be large when the last iteration produced a large change in the parameters compared to the local linear expectation and small when the last change was small. The iteration ends when the change in the parameters is small or the change in the ${\chi }^{2}$ value is small.

In the analysis, convolved galaxy fits, the galaxy model image, and the model derivative images must be convolved with the PSF at each iteration step. To save computation time, this convolution is performed using a circularly symmetric approximation of the PSF model, with the PSF model scale size set to the average of the major and minor axis direction scale size of the full PSF model, with the same radial profile term as the PSF model. The convolution is performed directly using the circular symmetry to reduce the number of multiplications performed: all points in the 2D circularly symmetric PSF model that have the same radial pixel coordinate can be evaluated in the convolution by summing up the corresponding pixels in the (galaxy model) image to be convolved before multiplying by the PSF model profile at that radial coordinate. This approximation reduces the number of multiplications by a factor of ∼8 for larger radii. For the small size of the PSF model used to convolve the galaxy model images, it was found that this direct convolution was faster than using an FFT-based convolution.

For the Exponential and de Vaucouleur fits, all parameters are fitted in the nonlinear minimization stage. For the Sérsic model, we do not fit the index within the Levenberg–Marquardt analysis. Instead, we start with a coarse grid search over a range of possible index values ($n=0.5,1.0,1.5,2.0,3.0,4.0,5.0,6.0$) and a range of possible values for ${R}_{\mathrm{eff}}$ based on the value of R₁, the first radial moment. For a given value of the Sérsic index, the ${R}_{\mathrm{eff}}$ is related to the first radial moment by the scale factor specified by Graham and Driver. We use the observed value of the 1st radial moment and try ${R}_{\mathrm{eff}}$ values of a factor of (0.8, 0.9, 1.0, 1.12, 1.25) times the value predicted by the Graham and Driver equation. For each of these steps, the aspect ratio and position angle are held constant and the normalization is determined to minimize the ${\chi }^{2}$.

We next perform three Levenberg–Marquardt minimization fits allowing the shape parameters (R_xx , R_yy , R_xy ) and the normalization to be fitted, holding the centroid (${x}_{0},{y}_{0}$), Sérsic index n and sky constant. In these fits, the index n is set to the minimum value previously calculated as well as values halfway to the next, and previous, values in the grid above. For example, if the minimum fitted index value is 3.0, then the LMM fits are performed using n = 2.5, 3.0, 3.5. The resulting ${\chi }^{2}$ values are then used to perform quadratic interpolation to find the index n which produces the locally minimum ${\chi }^{2}$ value. Finally, this best-fit index value is held constant while Levenberg–Marquardt minimization is used to find the best fit values of all other parameters. Sources for which a convolved galaxy model fit was successful have the flag bit PM_SOURCE_MODE_EXTENDED_FIT set.

The central pixel of the Sérsic, de Vaucouleur, and Exponential models requires special handling. When comparing an analytical model to the pixelized image recorded by a CCD, one normally treats the value in a pixel as equivalent to the value of the model at the center of the pixel. However, in reality, the number of counts observed in a pixel represents the integral of the surface brightness across the area of the pixel. This average will be equal to the central surface brightness times the area of a pixel as long as the second and higher derivatives of the analytical model are zero. However, if the first and second derivatives are nonzero, the curvature of the function within the pixel will make the integral differ from the central surface brightness times a fixed pixel area. If the curvature of the model function is sufficiently large, this difference will have a significant impact on the evaluation of the model. This situation is particularly true for the central portion of the Sérsic-like model functions.

In order to accurately compare the observed galaxy flux distribution to a model, it is necessary to integrate the pixel flux for a given set of model parameter values. In the psphot implementation, we currently use a brute-force numerical evaluation of the integral, dividing the central pixel into a grid of subpixels, with the sampling set by the Sérsic index of the model being evaluated as ${N}_{{\rm{sub}}}\,=\,2{\rm{Integer}}(6n/{R}_{min})$, where ${N}_{\mathrm{sub}}$ is the subpixel scale n is the Sérsic index and ${R}_{\min }$ is the size of the minor axis in pixel units. The value of ${N}_{\mathrm{sub}}$ is constrained to be in the range 11 to 121, so the number of subpixel evaluations ranges from 121 to ${121}^{2}=14,641$. Faster approximations to this analysis were explored, but they resulted in unsatisfactory results. This is definitely an area where psphot could benefit from some of the lessons in the literature (e.g., Hogg & Lang 2013).

The convolved galaxy model fit results are available in one of three PSPS database tables: StackModelFitExp, StackModelFitDeV, StackModelFitSer for the Exponential, de Vaucouleur, and Sérsic models, respectively.

For some science goals, a well-measured color of a galaxy is more important than an accurate total magnitude. In the case of PS1, the image quality variations for stacks of different filters presents a serious challenge for the determination of precise colors. psphot determines a set of PSF-matched radial aperture flux measurements in order to minimize the impact of the stack image quality variations.

In psphotStack, the stack analysis version of psphot, the five filter images are processed together. After the PSF models have been fitted and a best set of galaxy models have been determined, three sets of fixed circular apertures are measured. In the first set, the fluxes in the apertures are measured using the raw stack images. The centers of the apertures for each source across the five filters are fixed so that the pixels represent the equivalent portions of the same galaxy for all five filters. In this analysis, the best model for each source is subtracted from the image pixels for all sources excluding the source in consideration. The "best model" is determined based on the minimum ${\chi }^{2}$ value for the model fits.

In addition to the raw fixed circular apertures, the stack images are each convolved with a circular Gaussian with σ chosen to yield an image with a typical FWHM of 6 pixels (15). The full set of circular apertures is again measured on these convolved images. Again, the best source models are subtracted from the image for sources not being measured. This subtraction includes the convolution to smooth the model to the effective FWHM of the convolved image. The entire procedure is then repeated with a target FWHM of 8 pixels (2''). Note that we do not attempt to restrict the stack image quality to match these convolution targets. If the stack has an effective FWHM larger than either the 6 or 8 pixel targets, the convolution does not take place and the resulting analysis is performed on the raw stack. In such cases, the smallest apertures are sensitive to seeing variations and should be avoided. Of the individual images, the fraction with FWHM larger than 8 pixels is (${g}_{{\rm{P}}1}$, ${r}_{{\rm{P}}1}$, ${i}_{{\rm{P}}1}$, ${z}_{{\rm{P}}1}$, ${y}_{{\rm{P}}1}$) = (9.8%, 5.1%, 4.9%, 3.4%, 3.7%), and the fraction of stacks with an effective FWHM larger than this limit will be much smaller.

For the PV3 analysis of the 3π Survey data, the fluxes are measured for a set of up to nine circular apertures with sizes chosen to match the similar circular apertures measured by the SDSS analysis. These apertures have radii of (4.16, 7.04, 12.0, 18.56, 29.76, 45.68, 72.80, 112.80, 176.88) pixels = (104, 176, 300, 464, 744, 1142, 1820, 2820, 4422). If the object is too faint, the larger apertures will be largely noise, and the computation is wasteful. We only calculate the circular apertures out to the second aperture larger than the "sky radius" (defined in Section 4.6.3), but we calculate photometry for at least the smallest four apertures. Sources for which photometry in these fixed aperture is calculated have the flag bit PM_SOURCE_MODE_RADIAL_FLUX set. Although these apertures are chosen to match the SDSS apertures, the SDSS images are measured on unconvolved images. Because the median seeing for the SDSS images is ∼14 in the r band (Adelman-McCarthy et al. 2007), our 15 aperture photometry should generally compare well to the SDSS aperture magnitudes.

The measurements described in this subsection are presented in the PSPS database (Paper VI) in the StackApFlxExGalUnc, StackApFlxExGalCon6, StackApFlxExGalCon8, and StackApFlx tables. The first three tables present measurements for all apertures from the unconvolved, and 6 and 8 pixel FWHM convolved images (respectively) while the last table presents a subset of the radii from all three sets of measurements joined together for ease of access.

To test the galaxy model analysis, we have generated a series of simulated images containing both stars and galaxies on which we have run the psphot PSF-convolved galaxy model fitting analysis. The images generated for this analysis have dimensions of 4000 × 4000 pixels with a spatial scale of 025 per pixel. The images are generated using an effective exposure time of 30 s, with zero-points matching the PS1 ${r}_{{\rm{P}}1}$ filter and a realistic sky brightness of 20.86 mag per arcsec². The stars are injected into these images with fluxes drawn from a realistic stellar luminosity function and random spatial locations. For each image, the same underlying simulated stellar population was used. Galaxies are injected into the image with positions on a regularly spaced grid with a separation of 120 pixels. A total of 1089 galaxies are injected in each image, with all galaxies in an image having the same total magnitude. Galaxies were injected for 51 magnitude bins, ranging from 17.0 to 22.0. Given the parameters of the image, these values span the range from high S/N ($\gt 100$) to undetectable (S/N < 1). Note that we are not attempting to represent a realistic astronomical distribution of galaxies but rather attempting to understand our ability to recover galaxies of a certain type, size, and shape in a given image.

The galaxies are injected using Exponential and de Vaucouleur profiles in separate simulation runs. The major-axis values are randomly distributed between 1 and 10 pixels (025–25) while the aspect ratios are randomly chosen in a range from 0.25 to 1.0. The position angles are set by the sequence in the image and allowed to vary from 0° to 180°. The images are then convolved with a circular PSF model using the PS1_V1 profile (FWHM = 10, $\kappa \,=\,0.2$), and noise is added using Poisson statistics for the detected photons.

For the figures below, we present results as a function of the (input) instrumental magnitude of the galaxy minus the instrumental magnitude corresponding to the stellar 5σ detection limit. We make the simplifying assumption that the stellar detection threshold encapsulates enough information about the sensitivity of the images that this magnitude difference may be used to compare the results shown here to images with other depths. Thus, this and subsequent figures may be compared with the reported detection limits from the PS1 3π Survey. Note for reference that the typical stellar detection limits in the PS1 3π stack images (Paper I) are (${g}_{{\rm{P}}1}$, ${r}_{{\rm{P}}1}$, ${i}_{{\rm{P}}1}$, ${z}_{{\rm{P}}1}$, ${y}_{{\rm{P}}1}$) = (23.3, 23.2, 23.1, 22.3, 21.4). The minimum Kron magnitudes for which galaxy model fits were performed for the PV3 analysis (Section 5) thus correspond to −1.6 to −1.8 in these plots.

Figure 15 shows completeness for the detection of the Exponential and de Vaucouleur model galaxies. This analysis does not indicate if the galaxy was detected as a galaxy (i.e., was the extended nature of the source sufficiently clear), only if the source was detected by the peak-finding algorithm. As expected, the more compact galaxies are more likely to be detected; Exponential profile galaxies, with a broader light distribution for the same effective radius, are less likely to be detected for the same magnitude than de Vaucouleur profile galaxies. This completeness should be compared to our earlier work (Metcalfe et al. 2013) in which we injected a realistic population of simulated galaxies into real PS1 images. That work found that the 50% completeness for the typical galaxy was roughly 0.5 mag brighter than the 50% stellar completeness limit, somewhat fainter than the completeness shown in Figure 15. However, that previous work did not explore the dependency of the completeness on the galaxy size or profile. The difference suggests that the galaxies in the earlier work were generally compact.

Zoom In Zoom Out Reset image size

Figures 16 and 17 demonstrate the recovery of galaxy parameters in these simulations for galaxy using the Exponential and de Vaucouleur models, respectively. Both figures show the reliability of the measured magnitudes, major and minor axis sizes, and ellipticities, which we represent as $\tfrac{{R}_{\mathrm{Major}}-{R}_{\mathrm{Minor}}}{{R}_{\mathrm{Major}}+{R}_{\mathrm{Minor}}}$. Galaxies are grouped in 0.1 mag bins based on their injected magnitudes. For all recovered parameters, the standard deviation of the difference between the measured parameter and the truth value is shown for all galaxies in each magnitude bin, as well as for subsets in major-axis ranges. The mean of the difference, illustrating any biases, is also given for all galaxies. The comparison for the major and minor axis sizes are shown as absolute measurements (in pixels) and as a fraction of the axis in question.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some overall patterns can be observed. First, the parameters measured for the exponential profile galaxies are consistently more reliable than those measured for the de Vaucouleur profiles. Second, the errors in the estimated magnitudes are primarily driven by errors in the size measurements: if the sizes are overestimated, the fluxes are also overestimated. Finally, the magnitudes and major axes are more accurate for the smaller galaxies, but the ellipticities are more accurate for the larger galaxies.

PSF photometry is measured for all objects detected in the stack using the positions determined by the stack photometry analysis. Candidate PSF stars are pre-identified as those stars used to generate the PSF model in the stack photometry analysis. A PSF model is generated for each input warp image based on those stars; PSF stars that are excessively masked on a particular image are not used to model the PSF. The normalization of the PSF model is fitted to all of the known source positions in the warp images to determine the PSF fluxes. This measurement is performed simultaneously for all sources in the image at once using the method described above (Section 4.6.2).

Aperture fluxes, Kron fluxes, and moments are also measured at this stage for each warp. For the Kron fluxes, the radii are fixed to the value determined in the analysis of the stack. Fluxes are also measured in three of the fixed apertures discussed in Section 5.4: those with 300, 464, and 744 radii. Note that the flux measurement for a faint, but significant, source from the stack image may be at a low significance (less than the 5σ criterion used when the photometry is not run in this forced mode) in any individual warp image; the measured flux may even be negative due to statistical fluctuations. When combined together, these low-significance measurements result in a significant measurement as the S/N increases with the combination of more data.

Individual warp images are processed independently with separate executions of the psphotFullForce program. Sources that are loaded by psphotFullForce for analysis are marked with the flag bit PM_SOURCE_MODE_EXTERNAL. This bit is also used to mark user-supplied sources loaded for analysis by the regular version of psphot. Once all of the forced photometry measurements (for point sources as well as galaxies, discussed below) are completed for all of the warps that contributed to a stack image, the measurements are combined together by other portions of the IPP system. The PSF photometry measurements are combined in the context of the DVO database system (Magnier et al. 2020a), including recalibration of the zero-points for the individual warps. These measurements for each warp are available from the PSPS database ForcedWarpMeasurement and ForcedWarpExtended tables, the latter containing the three fixed-aperture fluxes. The average values calculated over the warps are found in the ForcedMeanObject tables.

With the inclusion of the forced-warp photometry, we have three distinct methods for measuring the PSF photometry of stars in the Pan-STARRS survey data: the average of the chip-stage photometry from the individual exposures, the measurement from the stacks, and the average of the forced-warp photometry described here. It is worth considering which of these should be used in which circumstance. Figure 18 shows a comparison of these three different methods to deeper data from the Medium-Deep Survey observations (MD06 field). Our conclusion from this and other analysis is that the average chip-stage photometry is the best (most accurate) measurement for brighter objects, where the signal-to-noise is roughly 10 or more. This is the photometry source that was used for the global photometry solution discussed by Schlafly et al. 2012 and used in the overall calibration (see Paper V).

Zoom In Zoom Out Reset image size

As can be clearly seen in the figure, the average from the forced-warp photometry is slightly worse than the chip photometry, while the stack PSF photometry is significantly degraded. We attribute the latter effect to the highly textured PSF observed in the stack images due to the combination of variable PSFs in each exposure and significant masking fraction in the PS1 camera. At the faint end, the chip photometry is significantly worse than both average warp and stack photometry. First, in order to have a measurement, a source must be detected above the detection threshold in at least one of the exposures, limiting the depth possible of the average chip photometry. Second, at the faint end, only bright fluctuations will be detected, resulting in a bright bias, a form of Eddington bias (Eddington 1913). This latter effect is clearly seen in Figure 18 as the average chip magnitudes diverge from the deeper Medium-Deep photometry measurements. As has been noted elsewhere (Best et al. 2018), the warp and stack photometry is also degraded for objects that have significant proper motion over the course of the 3π Survey because the position is held constant for all epochs, while the average chip photometry is calculated on detections that are cross-matched in the database. Thus, warp and stack photometry should be avoided for sources with proper motion greater than roughly 100 mas per year.

The convolved galaxy models are also remeasured on the warp images by the psphotFullForce analysis. In this analysis, the galaxy models determined from the stack image analysis are used to seed the analysis in the individual warp images. The motivation of this analysis is the same as the forced PSF photometry: the PSF of the stack image is poorly determined due to the masking and PSF variations in the inputs. Without a good PSF model, the PSF-convolved galaxy models are of limited accuracy.

In the forced galaxy model analysis, we assume that the galaxy position and position angle, along with the Sérsic index if appropriate, have been sufficiently well determined in the analysis of the stack image. In this case, the goal is to determine the best values for the major and minor axes of the elliptical contour and at the same time the best normalization corresponding to the best elliptical shape, and thus the best galaxy magnitude value. This technique is similar to the joint fitting of multiple exposures performed by the CFHT Lensing Survey team (Miller et al. 2013).

For each warp image, the stack values for the major and minor axes are used as the center of a grid search of the major and minor axis parameter values. The grid spacing is defined as a function of the S/N of the galaxy in the stack image so that bright galaxies are measured with a much finer grid spacing than faint galaxies. For the PV3 3π analysis, a 5 × 5 grid was used; values in both the major and minor axis directions of $(1-\tfrac{3.0}{S/N}$, $1-\tfrac{1.5}{S/N}$, 1.0, $1+\tfrac{1.5}{S/N}$, $1+\tfrac{3.0}{S/N})$ times the dimension are tested. For each grid point, the major and minor axis values at that point are used to generate the model. The model is then convolved with the PSF model for the warp image at that point. The resulting convolved model is then compared to the warp pixel data values and the best-fit normalization value is determined. The integrated flux, flux error, and the ${\chi }^{2}$ value for each grid point are recorded. This analysis is performed on the warp images one object at a time with the other objects in the image subtracted according to the best model currently known.

For a given galaxy, the result is a collection of ${\chi }^{2}$ values, fluxes, and flux errors for each of the grid points spanning all warp images. A single ${\chi }^{2}$ grid can then be made by combining each grid point across the inputs. The combined ${\chi }^{2}$ for a single grid point is simply the sum of all ${\chi }^{2}$ values at that point. If, for a single warp image, the galaxy model is excessively masked, then that image will be dropped for all grid points for that galaxy. The reduced ${\chi }^{2}$ values can be determined by tracking the total number of pixels used across all inputs to generate the combined ${\chi }^{2}$ values. From the combined grid of ${\chi }^{2}$ values, the point in the grid with the minimum ${\chi }^{2}$ is found. Quadratic interpolation is used to determine the major, minor axis values for the interpolated minimum ${\chi }^{2}$ value. The errors on these two parameters are then found by determining the contour at which the ${\chi }^{2}$ increases by 1.

In this way, the forced galaxy model analysis uses the PSF information from each warp image to determine a best set of convolved galaxy models for each galaxy model measured for the stack image. The results of these galaxy model fits are available from the PSPS database ForcedGalaxyShape table.

Weak-lensing studies frequently use nonparametric measurements of the ellipticities of galaxies to quantify the strength of gravitational lensing, and thus directly measure mass distributions in the universe. The classic approach was originally described by Kaiser et al. (1995, hereafter KSB), and applied to a set of deep HST observations. The details of the technique were further refined by Hoekstra et al. (1998, hereafter HFK); in the discussion below, we primarily use their notation, though we explicitly cast their integrals as sums over discrete pixels.

The KSB-style analysis of object ellipticities has also been used by several authors to search for marginally resolved binary stars in wide-field imaging data. The use of the lensing statistics for this application was described by Hoekstra et al. (2005) in the context of vetting planet transit events in data from the Optical Gravitational Lensing Experiment (OGLE). Terziev et al. (2013) applied the technique to PTF data to search for binary stars and Deacon et al. (2017) used the same technique to search for binary companions to known ultracool dwarfs using Pan-STARRS 3π data. The work by Deacon et al. (2017) used images and their own analysis of the pixels with the program Sextractor (Bertin & Arnouts 1996).

For the Pan-STARRS 3π PV3 analysis, we have measured the full set of KSB lensing parameters for all objects with measured second moments (i.e., excluding saturated stars, suspected cosmic rays, and other likely defects) of the data to enable both lensing studies and binary/multiple star searches. Here, we describe the measurements as performed within psphot, reviewing the mathematical framework as described by KSB and HFK. Just like the forced PSF photometry and the constrained galaxy models above, this analysis is performed by measuring the KSB lensing parameters on the individual warp images and averaging these measurements for each object.

The goal of the KSB technique is to measure the intrinsic ellipticity of objects (i.e., galaxies, in the case of weak lensing studies) as would be observed on the sky on the without instrumental effects and to determine the impact weak gravitational lensing would have on the observed shapes, after correcting for the instrumental effects. The analysis starts with the observed ellipticity of objects as represented by the two polarization components derived from the second moments (see Section 4.4.3):

$$\begin{eqnarray}&&{e}_{1}=\displaystyle \frac{{M}_{{xx}}-{M}_{{yy}}}{{M}_{{xx}}+{M}_{{yy}}},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{e}_{2}=\displaystyle \frac{2{M}_{{xy}}}{{M}_{{xx}}+{M}_{{yy}}}.\end{eqnarray} \tag{ 46 }$$

These two quantities have values that range from −1 to +1, with a circularly symmetric object having $({e}_{1},{e}_{2})=0.0,0.0$. The two polarization components vary sinusoidally with twice the position angle of the object: an elongated object aligned with the x-pixel axis will have positive values of e₁ and e₂ values near zero, while the same object aligned with the y-pixel axis will negative e₁ values. An object with a position angle on the 45° lines between the pixel axes will have large positive or negative values of e₂ and low absolute values of e₁.

Note that in our analysis of the second moments, we are applying a Gaussian window function to down-weight the noise contributions from pixels at high radii and low flux (see Section 4.4.3). This type of window function is also assumed in the KSB formalism and is represented in the equations below as W.

The measured ellipticity of an object observed in a real instrument will be affected by the PSF of the instrument. To first order, the effect on the polarization components can be described as a combination of the circularly symmetric seeing disk, which smears the observed shapes (driving ${e}_{1},{e}_{2}$ to low absolute values) and the shearing effect of the anisotropic component of the PSF, in which the observed shape is stretched in one direction relative to the others (driving ${e}_{1},{e}_{2}$ to larger absolute values).

KSB and HFK quantify the change in the observed polarization due to the smearing effect of the PSF with

$$\begin{eqnarray}&&\delta {e}_{\alpha }^{\mathrm{sm}}\,=\,{P}_{\alpha ,\beta }^{\mathrm{sm}}{p}_{\beta }.\end{eqnarray} \tag{ 47 }$$

p_β is a measurement of the anisotropy of the PSF (see below), and ${P}_{\alpha ,\beta }^{\mathrm{sm}}$ is the "smear polarizability" of the object, defined as

$$\begin{eqnarray}&&{P}_{\alpha \beta }^{\mathrm{sm}}\,=\,{X}_{\alpha \beta }^{\mathrm{sm}}-{e}_{\alpha }{e}_{\beta }^{\mathrm{sm}},\end{eqnarray} \tag{ 48 }$$

where

$$\begin{eqnarray}&&{X}_{1,1}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\sum f\left[W+2{W}^{{\prime} }{r}^{2}+{W}^{{\prime} ^{\prime} }{\left({x}^{2}-{y}^{2}\right)}^{2}\right],\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{X}_{1,2}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\sum f\left[2{W}^{{\prime} ^{\prime} }({x}^{2}-{y}^{2}){xy}\right],\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{X}_{2,2}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\sum f\left[W+2{W}^{{\prime} }{r}^{2}+4{W}^{{\prime} ^{\prime} }{x}^{2}{y}^{2}\right],\end{eqnarray} \tag{ 51 }$$

and

$$\begin{eqnarray}&&{e}_{1}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\sum f\left[2{W}^{{\prime} }+{W}^{{\prime} ^{\prime} }({x}^{2}+{y}^{2})\right]({x}^{2}-{y}^{2}),\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&{e}_{2}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\sum f\left[2{W}^{{\prime} }+{W}^{{\prime} ^{\prime} }({x}^{2}+{y}^{2})\right]2{xy}.\end{eqnarray} \tag{ 53 }$$

In these equations, $T={M}_{{xx}}+{M}_{{yy}}$ and W is the window function applied when measuring the second moments. The terms ${W}^{{\prime} }$ and ${W}^{{\prime} ^{\prime} }$ are the derivatives of the window function with respect to ${r}^{2}={x}^{2}+{y}^{2}$. Because the window function is a circularly symmetric Gaussian with width ${\sigma }_{w}$, the derivatives are simply ${W}^{{\prime} }=-\tfrac{1}{2{\sigma }_{w}^{2}}W$ and ${W}^{{\prime} ^{\prime} }\,=\,\tfrac{1}{4{\sigma }_{w}^{4}}W$.

The elements of the equations above can be written in terms of the second- and higher-order moments calculated in Section 4.4.3:

$$\begin{eqnarray}&&{X}_{1,1}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\left[1-\displaystyle \frac{{R}_{2}}{{\sigma }^{2}}+\displaystyle \frac{({M}_{{xxxx}}-2{M}_{{xxyy}}+{M}_{{yyyy}})}{4{\sigma }^{4}}\right],\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&{X}_{1,2}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\left[\displaystyle \frac{({M}_{{xyyy}}-{M}_{{xxxy}})}{2{\sigma }^{4}}\right],\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&{X}_{2,2}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\left[1-\displaystyle \frac{{R}_{2}}{{\sigma }^{2}}+\displaystyle \frac{{M}_{{xxyy}}}{{\sigma }^{4}}\right],\end{eqnarray} \tag{ 56 }$$

and

$$\begin{eqnarray}&&{e}_{1}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\left[\displaystyle \frac{{M}_{{xx}}-{M}_{{yy}}}{{\sigma }^{2}}+\displaystyle \frac{{M}_{{xxxx}}-{M}_{{yyyy}}}{4{\sigma }^{4}}\right],\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&{e}_{2}^{\mathrm{sm}}=\displaystyle \frac{1}{T}\left[\displaystyle \frac{({M}_{{xxxy}}+{M}_{{xyyy}})}{2{\sigma }^{4}}-\displaystyle \frac{2{M}_{{xy}}}{{\sigma }^{2}}\right],\end{eqnarray} \tag{ 58 }$$

where ${R}_{2}={M}_{{xx}}+{M}_{{yy}}$.

KSB and HFK use the observed ellipticities of stars and the smear polarizability of the stars to estimate the anisotropy due to the PSF:

$$\begin{eqnarray}&&{p}_{\alpha }=\displaystyle \frac{{e}_{\alpha }^{* }}{{P}_{\alpha \alpha }^{\mathrm{sm},* }},\end{eqnarray} \tag{ 59 }$$

where the terms with the ∗ represent parameters measured on stars.

Similarly, the impact of shear can be quantified by the "shear polarizability" in a similar fashion:

$$\begin{eqnarray}&&\delta {e}_{\alpha }^{\mathrm{sh}}={P}_{\alpha ,\beta }^{\mathrm{sh}}{p}_{\beta },\end{eqnarray} \tag{ 60 }$$

where now the shear polarizability ${P}_{\alpha \beta }^{\mathrm{sh}}$ is defined as

$$\begin{eqnarray}&&{P}_{\alpha \beta }^{\mathrm{sh}}={X}_{\alpha \beta }^{\mathrm{sh}}-{e}_{\alpha }{e}_{\beta }^{\mathrm{sh}},\end{eqnarray} \tag{ 61 }$$

where

$$\begin{eqnarray}&&{X}_{1,1}^{\mathrm{sh}}=\displaystyle \frac{1}{T}\sum f\left[2W({x}^{2}+{y}^{2})+2{W}^{{\prime} }{\left({x}^{2}-{y}^{2}\right)}^{2}\right],\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&{X}_{1,2}^{\mathrm{sh}}=\displaystyle \frac{1}{T}\sum f\left[4{W}^{{\prime} }({x}^{2}-{y}^{2}){xy}\right],\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&{X}_{2,2}^{\mathrm{sh}}=\displaystyle \frac{1}{T}\sum f\left[2W({x}^{2}+{y}^{2})+8{W}^{{\prime} }{x}^{2}{y}^{2}\right],\end{eqnarray} \tag{ 64 }$$

and

$$\begin{eqnarray}&&{e}_{1}^{\mathrm{sh}}=2{e}_{1}+\displaystyle \frac{2}{T}\sum {{fW}}^{{\prime} }({x}^{2}+{y}^{2})({x}^{2}-{y}^{2}),\end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray}&&{e}_{2}^{\mathrm{sh}}=2{e}_{2}+\displaystyle \frac{2}{T}\sum {{fW}}^{{\prime} }({x}^{2}+{y}^{2})2{xy}.\end{eqnarray} \tag{ 66 }$$

Rewriting in terms of the second- and higher-order moments calculated in Section 4.4.3, we find

$$\begin{eqnarray}&&{X}_{1,1}^{\mathrm{sh}}=\displaystyle \frac{1}{T}\left[2{R}_{2}-\displaystyle \frac{({M}_{{xxxx}}-2{M}_{{xxyy}}+{M}_{{yyyy}})}{{\sigma }^{2}}\right],\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray}&&{X}_{1,2}^{\mathrm{sh}}=\displaystyle \frac{1}{T}\left[\displaystyle \frac{2({M}_{{xyyy}}-{M}_{{xxxy}})}{{\sigma }^{2}}\right],\end{eqnarray} \tag{ 68 }$$

$$\begin{eqnarray}&&{X}_{2,2}^{\mathrm{sh}}=\displaystyle \frac{1}{T}\left[2{R}_{2}-\displaystyle \frac{4{M}_{{xxyy}}}{{\sigma }^{2}}\right],\end{eqnarray} \tag{ 69 }$$

and

$$\begin{eqnarray}&&{e}_{1}^{\mathrm{sh}}=\displaystyle \frac{1}{T}\left[2({M}_{{xx}}-{M}_{{yy}})+\displaystyle \frac{({M}_{{yyyy}}-{M}_{{xxxx}})}{{\sigma }^{2}}\right],\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{e}_{2}^{\mathrm{sh}}=\displaystyle \frac{1}{T}\left[4{M}_{{xy}}-\displaystyle \frac{2({M}_{{xxxy}}+{M}_{{xyyy}})}{{\sigma }^{2}}\right].\end{eqnarray} \tag{ 71 }$$

In the Pan-STARRS PV3 analysis, we have measured the elements of the smear polarizability (${X}_{\alpha \beta }^{\mathrm{sm}}$, ${e}_{\alpha }^{\mathrm{sm}}$) and the shear polarizability (${X}_{\alpha \beta }^{\mathrm{sh}}$, ${e}_{\alpha }^{\mathrm{sh}}$) for all objects on each of the warp images. We have also selected only the PSF stars from the images and interpolated a smoothed version of these parameters from the PSF stars to the locations of all objects, using the grid described above to interpolate the parameters. We then determine the ellipticities of the stars (${e}_{1}^{* },{e}_{2}^{* }$) from the equivalent smooth grid. Thus, for every object in the 3π Survey, we are able to report the PSF and object elements of the KSB analysis. These lensing parameters are measured for each of the warps, and then averaged over all warps for each of the filters. The average values are calculated by including only measurements from the same warp detection used in the average photometry (nominally, the primary skycell; see Paper V, Section 5.4.4) and excluding any measurements for which the PSF_QF or PSF_QF_PERFECT is less than 0.85.

The lensing parameters measured for individual warps are available from the PSPS database ForcedWarpLensing table while the average values calculated over the warps is found in the ForcedMeanLensing tables. Although the software used here was not involved in any of the GRavitational lEnsing Accuracy Testing (GREAT) challenges, it is similar to the code of the EPFL_KSB team (Mandelbaum et al. 2015) and likely to perform similarly.