PEARLS: Near-infrared Photometry in the JWST North Ecliptic Pole Time Domain Field^*

The reduction of the MMIRS imaging followed a procedure similar to those used by Labbé et al. (2003), McCracken et al. (2012), and Pelló et al. (2018), taking into account some of the characteristics of the H2RG detector. The first step in the reduction was subtracting the dark current from each science image plane using the average of the corresponding plane in the set of dark images. Next, image counts were corrected for nonlinearity using lookup tables compiled by the MMIRS team, cosmic rays were removed and crosstalk subtracted, after which the count-rate images ("slopes") for each ramp in units of data numbers per second (DN s⁻¹) were calculated. These reduction steps were carried out using the Chilingarian et al. (2015) IDL procedures ²¹ with adaptations for use with imaging data by one of the authors (C.N.A.W.), as briefly described in Appendix B. The procedures as modified are publicly available at https://github.com/cnaw/mmirs_imaging and in Willner et al. (2023).

The sky subtraction followed the IRAF script xdimsum, ²² which was specifically designed for the reduction of NIR data; similar algorithms were used by Thompson et al. (1999) and Huang et al. (2003). The first step in the process was removing a baseline sky value from each image and storing the reciprocal of this value. ²³ The initial sky value per slope image was obtained taking a pixel-by-pixel average using the 16 (baseline-subtracted) images taken nearest in time (both before and after a given exposure) with no masking of sources or bad pixels. In the case of images at the beginning or end of a sequence, we used the nearest images in time following or preceding a given observation (still totalling 16). The next step calculated the ratio of the image divided by the sky, from which the median ratio (after outlier rejection) was obtained. The sky image multiplied by this median ratio was then subtracted from the slope image.

Initial source catalogs per (sky-subtracted) image were calculated using SExtractor (Bertin & Arnouts 1996), and masks were created using the segmentation images, where each detection has its segmentation map extended radially by 2 pixels (042) to mask any residual source light close to the detection-level isophote. The initial catalogs were also used to identify "ghost" sources due to image persistence. This was done by examining the brightest sources and using the pixel coordinates between sequential images; the persistence images must fall ≲8 pixels away from the pixel position in the previous exposure occupied by the bright star and be fainter by more that 2 mag. In general, after two ramps, the remaining signal is negligible compared to the sky background, though in the case of very bright stars, some residuals can be detected even after ∼200 s (in the case of the Y exposures, which are the longest integrations we used). All persistence ghosts had their segmentation maps expanded radially by 5 pixels. Bad (i.e., dead or hot) pixels were also flagged, and their positions were noted in the segmentation images and added to the image masks for each slope image.

After the initial sky-subtraction step, we measured the seeing of individual images with the aim of optimizing the mosaic depth relative to the image resolution. To estimate the seeing, we used the sky-subtracted images to identify stellar sources in each image using distributions of the instrumental surface brightness of the brightest pixel (μ_max) versus the detected area (isoarea_image) to remove hot pixels and residual cosmic rays, as well as the half-light radius (flux_radius) versus the "total" magnitude (mag_auto), all calculated by SExtractor. The total magnitudes used a minimum Kron radius of 3.5 pixels (147 diameter) and a Kron factor of 2.5. The average FWHM for each image was calculated from the measurements of these candidate stars using the biweight estimator (Beers et al. 1990).

Figure 2 shows the number of images in each filter for a given FWHM. The left panel shows the distribution for the Y, J, and H filters; the right panel shows the distribution for the observations using MMIRS with the K_spec, K_0.95, and K_2.68 configurations. The majority of observations peak under ∼1'' values, and the distribution of K_spec, K_0.95, and K_2.68 reflects the varying conditions during the observations, taken over several nights between 2017 June and 2017 November (see Figure 10 in Appendix A for the variations over individual nights). In the case of K_0.95, data were taken in a single night, while those for K_spec were taken over the course of two nights. Most of the K values are ≤15 but have a somewhat broad distribution, in particular in the case of the K_2.68 observations.

Zoom In Zoom Out Reset image size

A second calculation of the sky background was then performed using the object masks, after which we refined the World Coordinate System (WCS) of individual images using the http://astrometry.net software of Lang et al. (2010) and the Gaia-DR3 catalog (Gaia Collaboration et al. 2023) as reference. The improved astrometry for individual images has typical rms residuals ∼003. During this phase of sky-subtraction, weights for individual images, described by the inverse variance due to the total number of photons in an exposure, the readout noise (3.14 e⁻), and the gain, have also aggregated a contribution due to the seeing:

$$\begin{eqnarray}&&\begin{array}{l}{\mathrm{weight}}_{i,j}\\ \quad =\,\displaystyle \frac{1}{({\mathrm{count}\_\mathrm{rate}}_{i,j}\times \mathrm{exptime}+\mathrm{readnoise}/{\mathrm{gain}}^{2})\times {\mathrm{FWHM}}^{2}}.\end{array}\end{eqnarray} \tag{ 1 }$$

Masked pixels were assigned inverse-variance values of zero, so swarp (Bertin et al. 2002; Bertin 2010) ignores these pixels when assembling mosaics, while the lowest-weight values were set to 10⁻⁹ for the SExtractor parameter weight_thresh; pixels with associated weights lower than this value are ignored by SExtractor. With the image masks and weights in place, another round of object detection was carried out. These new catalogs were used by scamp (Bertin 2006) to include higher-order terms in the astrometry, providing WCS keywords compatible with swarp. ²⁴ The astrometric calibration of scamp also used the Gaia-DR3 catalog matched to sources classified as stellar on individual slope images. The zero-point magnitude offsets between individual slope images as well as initial ensemble zero-point offsets were calculated using sources matched to 2MASS (Skrutskie et al. 2003, 2006) for J, H, K, and K_spec, and to Pan-STARRS (Tonry et al. 2012) for the Y band. These corrections by scamp were written into ASCII headers subsequently read by swarp when creating the full-field mosaics.

To create the full-field mosaics, we followed a procedure similar to those of Ashcraft et al. (2018, 2023) and McCabe et al. (2023), where images in a given filter are combined in two different ways: one to optimize the spatial resolution at the cost of depth, and the other to optimize the depth at the cost of resolution. In those works, the optimal resolution images result from stacking the best-seeing images, comprising about 10% of the total number of images (Ashcraft et al. 2018), while for the optimal depth, the 5%–10% worst-seeing images are excluded (Ashcraft et al. 2018; McCabe et al. 2023); typically, the cutoff is made at ∼20, with the deeper mosaic reaching ∼1 mag fainter than the higher-resolution mosaic (McCabe et al. 2023).

Because the aim of the NIR imaging was to provide a first epoch of variability studies, we opted to reach as faint as possible, and following these works, we cut at 21 when stacking the images. The greatest impact of this cut is on the K_2.68 images, where 213 out of 2571 (∼8%) were excluded. These numbers are smaller for the remaining setups—26/796 (∼3%) for K_spec, 7/563 for K_0.95, 3/1502 for H, and 1/248 for Y. No J images were excluded.

In addition to mosaics of individual images, stacked images combining the Y + J, K + K_spec, and H + K + K_spec were created to have deeper (or in the case of K + K_spec, full-field) coverage. In addition to the image mosaics, swarp calculated a stacked weight image that includes the contribution from individual image weights. For the output mosaics, we adopted the pixel size of 0168 and center position (17:22:48.1298, +65:50:14.853, J2000) used in the Subaru/HSC imaging discussed in Section 4.3. The final mosaics were combined using swarp parameters combine_type=clipped with clip_sigma=3.0.

The source catalogs were generated in a three-step process. The first step used SExtractor with a high detection threshold (5σ) per pixel relative to the average background, optimized to detect stars; these candidate sources were then processed by PSFEX (Bertin 2011) to calculate the average point-spread function (PSF) for each mosaic. The PSFs were used in a second pass of SExtractor to calculate the photometry and image classification using the parameters listed in Table 2. All runs used the dual-image mode, where the detection images were mosaics combining J + Y for J and Y, H + K + K_spec for H, and combined K + K_spec for the K and K_spec images. The positions of these catalogs are in the Gaia-DR3 system, and using Gaia-DR3 stars we find that the MMIRS positions have uncertainties ≲0070. The photometric zero-points reported in Table 2 were calculated following Almeida-Fernandes et al. (2022), who used external photometric measurements coming from wide-area surveys (e.g., Gaia and SDSS) with the stellar population models of Coelho (2014) to calculate the unknown zero-point offsets of uncalibrated filters. This procedure creates a grid of flux measurements by folding the stellar population models with the filter curves of the reference data and considering several values of Galactic extinction. By means of χ² minimization, the best-fitting model is found that also provides an estimate of the Galactic (and atmospheric) extinction. From the ensemble of fitted stars, the best-fitting zero-point is then estimated.

In this calibration, we matched the MMIRS-derived catalogs with Gaia, SDSS (Abdurro'uf et al. 2022), Pan-STARRS, and 2MASS, and only used sources that were classified as stars in Gaia and also had measurements in all other catalogs. While the number of sources used is small (ranging from 32 in the case of K_spec to 136 for J), in particular compared to the thousands used by Almeida-Fernandes et al. (2022), the average values for the zero-point estimates show uncertainties of only ∼1%–4%. These are much smaller than the estimates calculated using scamp, which are on the order of 10%–30% (per image) when using 2MASS measurements alone. These final zero-point values are shown in the third column of Table 1.

To create the catalogs for each MMIRS filter, we used exposure maps to remove sources contained in regions covered by fewer than (5, 4, 7, 9, 9, 9) images for Y, J, H, K_0.95 , K_2.68 , K_spec because of the low signal-to-noise ratio and greater likelihood of including spurious detections. These regions are those outlined in Figure 1.

The completeness of individual bands and for the K + K_spec mosaic was estimated following, e.g., Caldwell (2006) and Finkelstein et al. (2015), by adding simulated sources to cutouts of the science images. We used sections of 1000 × 1000 pixels—with the exception of the K + K_spec mosaic, for which a 2000 × 2000 region was used—centered on the position where all K imaging setups overlap. Using a grid of 0.1 mag steps, 100 simulations per magnitude bin were run, where in each simulation 100 sources were added at random positions and catalogs created using the same procedure as for the science catalogs, i.e., using a minimum area of 10 pixels with a detection threshold at 1σ above the mean background and using the SExtractor mag_auto estimator for the total magnitude and the detection parameters presented in Table 2. Table 3 lists the completeness limits in each band estimated for point sources modeled by the PSF constructed by PSFEX from the mosaics and for sources with a Sérsic n = 1 profile, using GalSim (Rowe et al. 2015) models convolved with these PSFs. The Sérsic models were sampled from a uniform distribution of position angle, half-light radius, and axial ratio ranging between [0°, 360°], [08, 33], and [0.3, 1), respectively. Table 3 shows the completeness limits (in AB magnitudes) at the 95%, 80%, and 50% levels for each of the filters (and gain value in the case of K) and for the combined K + K_spec mosaic. The completeness curves from which these limits were computed are presented in Figure 3.

Zoom In Zoom Out Reset image size

Table 3. Completeness Levels

Filter	mAB(95%)	mAB(80%)	mAB(50%)	Source Type, Minarea a , Threshold b
Y	23.80	24.22	24.59	point source, 10, 1
Y	22.32	22.82	23.23	Sérsic n = 1, 10, 1
J	23.53	24.03	24.38	point source, 10, 1
J	21.80	22.57	23.01	Sérsic n = 1, 10, 1
H	23.13	23.52	23.88	point source, 10, 1
H	21.63	22.13	22.55	Sérsic n = 1, 10, 1
K0.95	22.71	22.97	23.32	point source, 10, 1
K0.95	21.26	21.59	22.00	Sérsic n = 1, 10, 1
K2.68	22.91	23.34	23.68	point source, 10, 1
K2.68	21.39	21.91	22.35	Sérsic n = 1, 10, 1
Kspec	23.33	23.58	23.93	point source, 10, 1
Kspec	21.83	22.19	22.58	Sérsic n = 1, 10, 1
K + Kspec	23.28	23.62	24.02	point source, 10, 1
K + Kspec	21.79	22.25	22.69	Sérsic n = 1, 10, 1

Notes.

^a SExtractor minimum detection area (in pixels² of 0168 per pixel). ^b Sextractor detection threshold in terms of standard deviations above the background level.

Download table as:  ASCIITypeset image

The different filter throughputs and detector gain values have a small but detectable effect on the data quality of the mosaics. We show in Figure 4 the distribution of the standard deviation of the background flux measured from 100,000 randomly placed 21 apertures in each of the K-band mosaics. The K_2.68 measurements (orange distribution) show the smallest dispersion, though there is a large amount of overlap between these and the measurements using the wider K_spec filter (gray) and the K_0.95 (green) setups. Because of this, we created catalogs for the individual instrument settings. The different peaks seen for the K_2.68 and K_0.95 setups are due to a large number of pixels that have very similar exposure times and are concentrated at the center of the mosaics in a given quadrant.

Zoom In Zoom Out Reset image size

The regions where each of these settings overlap (see Figure 5) were used to estimate average offsets between the different imaging samples (Table 4 and Figure 6). In this calculation, we used the IDL task resistant_mean with a cutoff at 3.5σ and restricted the sample to sources with Sextractor flags = 0. We considered two cases: in the first one, we used sources with magnitudes in the range 18 ≤ mag_auto ≤ 23 mag to exclude bright saturated sources and faint sources where the photometry becomes uncertain; the second case repeated the same calculation but without magnitude cuts. From Table 4, using the bright sample, we found an offset of −0.031 ± 0.007 mag between K_2.68 and K_0.95. The offset between K_0.95 − K_spec is 0.122 ± 0.010 mag. Combining the K_2.68 − K_0.95 and K_0.95 − K_spec differences implies that K_2.68 − K_spec ∼ 0.091 mag. The direct measurement of K_2.68 − K_spec = 0.073 mag suggests that constructing a single-magnitude sample can have systematic uncertainties on the order of ∼0.02 mag.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  K Magnitude Differences

Filter and Gain Setup	Mean Difference	Uncertainty a	σ	Matching Sources	Bright (B) or Full (F) Sample	Quadrants
K2.68 − K0.95	−0.031	0.007	0.282	1613	B	SE, SW
K2.68 − Kspec	0.073	0.017	0.418	621	B	SE
K0.95 − Kspec	0.122	0.010	0.340	1284	B	SE
K2.68 − K0.95	−0.010	0.007	0.307	1921	F	SE, SW
K2.68 − Kspec b	0.130	0.018	0.513	776	F	SE
K0.95 − Kspec	0.212	0.011	0.442	1686	F	SE

Notes.

^a The uncertainty is calculated by resistant_mean as $\sqrt{\sigma }/(N-1)$. ^b Both sets use gain of 2.68 e⁻ DN⁻¹.

Download table as:  ASCIITypeset image

To create a uniform sample, we used K_2.68 as reference magnitudes, subtracted 0.031 mag from K_0.98, and added 0.073 mag to the K_spec measurements. The final magnitude for each source is an inverse-variance-weighted average using the corrected fluxes in each modality with weights derived from the estimated flux uncertainty measured by SExtractor.

To expand the wavelength coverage and provide visible counterparts to the NIR sources, we used archival Subaru Hyper-Suprime-Cam (HSC) data that were available in early 2020 from the SMOKA ²⁵ archive maintained by the National Astronomical Observatory of Japan. These images were obtained as part of HEROES (Songaila et al. 2018; Taylor et al. 2023) and comprise data of a single field: NEP-wide-A05. The observation log is shown in Table 5, which lists the filter name, total exposure time, the limiting magnitude estimated from the peak of the magnitude distribution in 0.5 mag bins, filter pivot wavelength and bandwidth (e.g., Willmer 2018), the flux corresponding to a count of 1 photon per second, and dates of observation.

The data were reduced using the HSC pipeline (Bosch et al. 2018), and for the NEP TDF analyses, a region of ∼3080 × 3359 centered at (17:22:48.1298, +65:50:14.853, J2000) was excised, which contains the footprints of JWST (PI: R. Windhorst & H. Hammel) and HST (PI: R. Jansen & N. Grogin) imaging and most of the other ancillary data, e.g., JCMT/SCUBA2 and VLA, (Hyun et al. 2023); Chandra (PI: W. Maksym); NuSTAR (PI: F. Civano; Zhao et al. 2021), and XMM-Newton (PI: F. Civano). The source catalog was created using SExtractor (Bertin & Arnouts 1996) in dual-imaging mode with a detection image stacking the i2, z, NB0816, and NB0921 images following Szalay et al. (1999). The detection was run in two steps, "cold" and "hot," where some of the detection parameters are changed to minimize the fragmentation of large and bright galaxies; these parameters are presented in Table 6 and those that change according to the detection mode are noted. The final PEARLS HSC catalog is a single list concatenating the "cold" and "hot" output from the five individual bands, and it contains data for 56,786 sources. The astrometric and photometric calibration used the same process described by Bosch et al. (2018) and Aihara et al. (2022) with the Pan-STARRS DR2 as reference, and it is tied to the Gaia DR1 astrometric reference frame (Gaia Collaboration et al. 2016a, 2016b). The HSC images are calibrated into maggies, ²⁶ which are in the AB system (Oke & Gunn 1983), and the flux value corresponding to a magnitude of zero is specified by the FLUXMAG0 image header keyword. ²⁷ The conversion from maggies into Jy or erg cm⁻² s⁻¹ Hz⁻¹ requires multiplying the image values by 10^8.9/2.5 or 10^−48.6/2.5, respectively. Because of the flux units adopted by the HSC pipeline, the SExtractor magnitudes are calculated directly in the AB system. However, the expression used by SExtractor to calculate the magnitude uncertainties,

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{magerr}=\displaystyle \frac{2.5}{\mathrm{ln}(10)}\\ \quad \times \,\displaystyle \frac{\sqrt{\mathrm{area}\times \mathrm{variance}(\mathrm{background})+\mathrm{flux}/\mathrm{gain}}}{\mathrm{flux}},\end{array}\end{eqnarray} \tag{ 2 }$$

assumes flux measurements in DN s⁻¹ (data numbers per second), and it will produce incorrect uncertainties (by several orders of magnitude) because the statistical uncertainty in photon counts has to be separated from the conversion of photon counts into flux units. The magnitude uncertainties we report are obtained by converting the AB magnitudes into photon numbers per cm² per second using

$$\begin{eqnarray}&&\mathrm{photon}\_\mathrm{numbers}=\displaystyle \frac{{f}_{\nu 0}\times d\lambda }{h\times \lambda },\end{eqnarray} \tag{ 3 }$$

where f_ν 0 is the flux measured in erg cm⁻² Hz⁻¹ s⁻¹ corresponding to one event per second in a given filter, λ and dλ are the filter's characteristic wavelength and bandwidth, and h is Planck's constant. We make the assumption that the instrument gain is already contained in the calibration, and we use the pivot wavelength (e.g., Bessell & Murphy 2012, Equation (A13)) as the characteristic wavelength. Assuming Poisson errors both for background and source photons, the expression above becomes

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{magerr}=\displaystyle \frac{2.5}{\mathrm{ln}(10)}\\ \quad \times \,\displaystyle \frac{\sqrt{[\mathrm{area}\times \mathrm{variance}(\mathrm{background})+\mathrm{flux}]\times \mathrm{photon}\_\mathrm{numbers}\times \mathrm{mirror}\,\mathrm{area}}}{\mathrm{flux}\times \mathrm{photon}\_\mathrm{numbers}\times \mathrm{mirror}\,\mathrm{area}},\end{array}\end{eqnarray} \tag{ 4 }$$

where area is the total number of pixels covered by the source, and the flux as well as the background variance are values measured by SExtractor. Using the filter parameters in Table 5 and an effective diameter of 820 cm for the Subaru telescope, ²⁸ this expression produces estimated uncertainties that are on the same order of magnitude as those available in the SDSS database for objects in common and are slightly more conservative than the uncertainty estimates in the full HEROES catalog (Taylor et al. 2023), as discussed in the next paragraph. Both flux and flux uncertainties are converted into nJy by multiplying each value by 10^31.4/2.5.

Table 6. HSC SExtractor Detection Parameters

Parameter	Value	Note
DETECT_THRESH	5.0	σ above mean background cold mode
ANALYSIS_THRESH	5.0	σ above mean background cold mode
DETECT_MINAREA	10	pixels
DEBLEND_NTHRESH	32
DEBLEND_MINCONT	0.001	cold mode
DETECT_THRESH	1.10	σ above mean background hot mode
ANALYSIS_THRESH	1.10	σ above mean background hot mode
DETECT_MINAREA	10	pixels
DEBLEND_NTHRESH	32
DEBLEND_MINCONT	0.0001	hot mode
FILTER_NAME	Gauss_1.5_3 × 3.conv
PHOT_FLUXFRAC	0.5
PHOT_AUTOPARAMS	2.5,3.5	Kron factor, minimum radius
PHOT_AUTOAPERS	0.0,0.0
CLEAN_PARAM	0.3
WEIGHT_TYPE	BACKGROUND
RESCALE_WEIGHTS	N
WEIGHT_THRESH	1.e-08
SEEING	07
SATUR_LEVEL	8.0e-09	maggies

Download table as:  ASCIITypeset image

A catalog for the full HEROES survey was published by Taylor et al. (2023), who also used the pipeline described by Bosch et al. (2018) to produce a band-merged aperture-matched catalog for the full HEROES region. The HEROES catalog divides the survey into several patches, four of which (17,9, 17,10, 18,9, and 18,10) cover the NEP TDF.

The catalog of unique sources (HEROES_Full_Catalog.fits) was used to recover some of the photometric measurements only available in the patch-level catalogs, e.g., aperture magnitudes. We restricted this catalog to the NEP TDF region: 260.0583 (17:20:13.992) ≤ R.A. ≤ 261.3350 (17:25:20.400) and 65.550 (65:33:00) ≤ decl. ≤ 66.125 (67:07:30). It contains 150,216 sources after setting the flags

$$\begin{eqnarray}&&{\mathtt{is}}\_{\mathtt{primary}}=1,\end{eqnarray} \tag{ 5 }$$

and for each filter

$$\begin{eqnarray}&&\begin{array}{l}\{g,i2,z,{NB}816,{NB}921\}\_{\mathtt{base}}\_{\mathtt{PixelFlags}}\_{\mathtt{flag}}\_{\mathtt{edge}}\\ \,\,=\,{false},\end{array}\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}\{g,i2,z,{NB}816,{NB}921\}\_{\mathtt{base}}\_{\mathtt{PixelFlags}}\_{\mathtt{bad}}\\ \,\,=\,{false},\end{array}\end{eqnarray} \tag{ 7 }$$

following Taylor et al. (2023). Matching the latter with the PEARLS catalog of HSC sources produces a list of 55,796 objects in common. To assess the quality of the PEARLS photometry, we selected stellar sources using the distribution of HEROES MAG_PSF magnitudes versus the MAG_PSF - MAG_CMODEL magnitude difference following Bosch et al. (2018). In the comparison that follows, we restricted sources to having

$$\begin{eqnarray}&&\parallel {\mathtt{MAG}}\_{\mathtt{PSF}}-{\mathtt{MAG}}\_{\mathtt{CMODEL}}\parallel \leqslant 0.01\end{eqnarray} \tag{ 8 }$$

for all wide-band filters and only considered sources with 17.5 ≤ i2 ≤ 23.5.

This excised catalog was used to verify the SExtractor photometry after accounting for some of the procedure differences. For example, the HSC pipeline uses pixel radii to define the aperture photometry, whereas SExtractor uses diameters. The HSC pipeline also adds an aperture correction when producing the final magnitudes for stars and galaxies, but it is not used in the aperture magnitudes (Bosch et al. 2018). The photometry comparison used the aperture magnitude measured at a radius of 12 pixels, which corresponds to a 403 aperture diameter, as this showed the best compromise between the amount of light being measured and the dispersion between both sets of measurements. The results from this comparison are shown in Table 7, where the mean value and estimated uncertainty of the differences in magnitude are calculated using the IDL resistant_mean procedure. The last column shows the multiplicative factor that corrects the fluxes measured by SExtractor to be consistent with the HEROES values.

To create a merged catalog joining the PEARLS HSC and MMIRS photometry, sources were matched using a search radius of 08 between PEARLS HSC and the NIR sources and between the different NIR catalogs. While the position uncertainties of the MMIRS and PEARLS HSC are smaller than this—0070 and 0030, respectively—the larger search radius allows matching the larger and brighter objects. In the cases where there were multiple matches, these were visually inspected and resolved in the cases of obvious mismatches. The final list of unique sources includes those with matches as well as sources with detections in single NIR bands. The total number of sources in the final merged catalog is 57,501 (including 56,752 PEARLS HSC sources and respectively (Y, J, H, K) = (8128, 10,558, 7612, 6361). Of the sources identified only in the MMIRS imaging, 354 were found to be spurious and removed from the sample. The number of sources with a detection in a single NIR band is (50, 146, 151, 73) for Y, J, H, K, respectively. As the number of single detections is small, these were visually inspected to remove spurious sources due to hot pixels, residual cosmic rays, background fluctuations, and stellar diffraction spikes. The final tally is (1, 32, 82, 35) sources in each of the NIR filters, but we note that several sources have PEARLS HSC counterparts that are below the catalog detection threshold. The larger number of single detections in H is due to the use of the combined H and K images in the source detection and the better image quality of the observations (e.g., Figure 2). Once the spurious detections are removed, the remaining number of sources is (Y, J, H, K) = (7990, 10355, 7535, 6316), and the total number of bona fide sources is 57,467, of which 715 are NIR-detected without counterparts in the PEARLS HSC catalog.

The number of sources without PEARLS HSC detections ²⁹ but present in all NIR (Y, J, H, K) catalogs ("z-dropouts") is 118; the number of those with H and K ("J-dropouts") is 323. The latter comprise candidate objects close to bright sources that were not deblended correctly in the PEARLS HSC catalog (20 sources, ∼6% of the unmatched subsample), sources with PEARLS HSC counterparts below the detection level (96 sources, ∼29.5%), and legitimate higher-redshift sources (210 galaxies, ∼64.5%). Given the relatively bright limits of these samples, these "dropouts" are most likely caused by the Balmer break shifting through the HSC and MMIRS filters. To better understand the properties of these sources without PEARLS HSC matches, we matched the PEARLS HSC-MMIRS catalog with a preliminary catalog of the NEP TDF comprising the first two spokes, which contains 24,119 sources. Using a 050 matching radius, we identify 1858 counterparts, 1419 of which have MMIRS measurements in at least one filter. Figure 7 shows the 174 sources with no PEARLS HSC counterparts but detected in at least one MMIRS filter. Of these, eight have only K measurements (red circles), 66 have H and K (orange triangles), and the remaining 29 are detected in all NIR filters (blue triangles). The gray circles represent the remaining 69 sources that have no HSC detection and have one or two NIR detections but not necessarily in contiguous filters; in one case, this is a source contaminated by a diffraction spike. The visual inspection of these objects suggests the majority are higher-redshift (z ≳ 1.5) or dusty galaxies—five of these galaxies are also SCUBA2 sources discovered by Hyun et al. (2023). One source is a quiescent galaxy at z ∼ 2 (N. Adams 2023, private communication).

Zoom In Zoom Out Reset image size

We used the HEROES (Taylor et al. 2023) subcatalog discussed in Section 4.3 to have an independent assessment of the 715 galaxies without PEARLS HSC matches. HEROES includes additional photometry specifically obtained in the NEP TDF with significantly deeper imaging in the z band and additional data in the HSC-r2 and HSC-y bands. We find 562 sources with HEROES counterparts within 05, while 153 have no HEROES match. Of these, 12 are detected in all MMIRS bands, 93 in H and K, and 20 in K only; five are detected in J, where there are no MMIRS observations in the other filters. There are 17 objects detected in Y and J only, and most of these are due to imperfect matches, because of merged images or sources close to bright stars. We also matched this catalog to the JWST catalog of spokes 1 and 2, and we find 28 objects in common, which are identified in Figure 7 as green diamonds. Four of these are detected in K only, and the majority (16) are detected in both H and K.

The matched PEARLS HSC-MMIRS-JWST catalog was also used to verify if any sources show evidence of variability, comparing measurements in the closest pair of filters (e.g., F090W − z, F115W − J, F150 − H, and F200W − K). When considering only the unsaturated bright sources, there is no evidence of variability. The sources presenting large magnitude differences are all due to contamination either by diffraction spikes or very close bright objects.

The PEARLS HSC+MMIRS merged catalog was used to classify sources detected in the NIR bands as stars or galaxies (or more correctly, compact or extended sources). For this, we used the total magnitude versus the half-light radius, e.g., Kron (1980), measured from the HSC z band and the total magnitude (mag_auto) in each NIR band. We used the radii measured for the PEARLS HSC because of the greater depth and overall better image seeing. In the case of MMIRS sources without PEARLS HSC half-light radius measurements, we used the corresponding NIR-band half-light radius values. The results from this classification are shown in Figure 8, where sources of different types are distinguished by symbol shapes and colors as noted in the key of each panel. We also show the 95%, 80%, and 50% completeness levels for extended sources as dashed, dotted–dashed, and dotted lines. The plots show that, below the 80% completeness level, the star/galaxy classification starts breaking down, though this changes with the filter being used.

Zoom In Zoom Out Reset image size

In addition to characterizing the large-scale properties of a sample (e.g., Koo & Kron 1992; Smail et al. 1995), the number counts (shown in Table 8) provide an independent check on the quality of the photometry, the effectiveness of the star/galaxy separation (e.g., Conselice et al. 2008), and the completeness of the catalog (e.g., Windhorst et al. 2023). Figure 9 shows the number counts ${N}_{c}(m)=\tfrac{N(m)}{{\rm{\Delta }}m\times \mathrm{area}}$ measured in the NEP TDF for the Y, J, H, and combined K samples. The uncertainties in Table 8 assume  Poisson  statistics  $\delta {N}_{c}{(m)=(\sqrt{N(m)}\times \mathrm{area}\,-N(m)\times {\rm{\Delta }}\mathrm{area})/({\rm{\Delta }}m\times \mathrm{area})}^{2}$ and are underestimated at larger numbers of galaxies because no cosmic variance is included. Vertical lines representing the completeness levels for extended sources from Table 3 are also shown and indicate that the turnover due to incompleteness is well characterized by the Monte Carlo simulations. The figures also show a selection of published measurements noted in Table 9, and when required, the latter were converted into AB magnitudes by adding the Vega to AB offsets of (J, H, K) = (0.870, 1.344, 1.814) mag derived for the 2MASS (Skrutskie et al. 2006) filters by Willmer (2018). The published counts in N_gal 0.5 mag⁻¹ per unit solid angle were converted into N_gal mag⁻¹ deg⁻². The areas and limiting magnitudes come from the faintest magnitudes in the original references. ³⁰ Given the small area covered by the NEP TDF MMIRS imaging—(0.0635, 0.0870, 0.0643, 0.0634) deg² for (Y, J, H, K), respectively)—some differences relative to other measurements are expected, though the agreement is good for magnitudes of 18 and fainter in all four bands.

Zoom In Zoom Out Reset image size

The deepest data in Y, J, and H come from the HST imaging in GOODS-S by the Early Science Release using the Wide-Field Camera 3 of Windhorst et al. (2011). In K, the deepest coverage comes from the ground-based observations of Fontana et al. (2014) in the UDS and GOODS-S fields. Because of the large number of samples in K, it is also the band showing the greatest variation between samples. The greatest differences are seen for the K number counts at bright magnitudes by Bielby et al. (2012) and Fontana et al. (2014), because of small number statistics. Fainter than K ∼ 18, the largest differences are seen for the Iovino et al. (2005) and Conselice et al. (2008) counts. However, the J counts of both references are good agreement with other works. For J, H, and K, the MMIRS counts show good agreement with the multifield number counts of Keenan et al. (2010) throughout the common magnitude range.