PEARLS: Near Infrared Photometry in the JWST North Ecliptic Pole Time Domain Field

We present Near-Infrared (NIR) ground-based Y, J, H, and K imaging obtained in the James Webb Space Telescope North Ecliptic Pole Time Domain Field (TDF) using the MMT-Magellan Infrared Imager and Spectrometer (MMIRS) on the MMT.These new observations cover a field of approximately 230 arcmin^2 in Y, H, and K and 313 arcmin^2 in J. Using Monte Carlo simulations we estimate a 1 sigma depth relative to the background sky of (Y, J, H, K}) = (23.80, 23.53, 23.13, 23.28) in AB magnitudes for point sources at a 95% completeness level. These observations are part of the ground-based effort to characterize this region of the sky, supplementing space-based data obtained with Chandra, NuSTAR, XMM, AstroSat, HST, and JWST. This paper describes the observations and reduction of the NIR imaging and combines these NIR data with archival imaging in the visible, obtained with the Subaru Hyper-Suprime-Cam, to produce a merged catalog of 57,501 sources. The new observations reported here, plus the corresponding multi-wavelength catalog, will provide a baseline for time-domain studies of bright sources in the TDF.


INTRODUCTION
One of the greatest discoveries in the last few decades is the existence of Dark Energy, which was inferred from the observations of Type Ia supernovae ranging from the local universe to redshifts of 1 and above.This discovery has prompted several large scale surveys of galaxies, e.g., SDSS-IV (Albareti et al. 2017), Palomar Transient Factory (Rau et al. 2009), J-PAS (Bonoli et al. 2021), and DESI (Dey et al. 2019) among others, to better characterize Dark Energy and understand its origin.Because the detection and follow-up of very distant supernovae from the ground becomes very challenging at redshifts of z ∼ 2 and above, the use of space-based facilities becomes essential.However, a major limitation in using the space observatories is the ability to monitor regions of the sky on a somewhat regular time cadence, that allows identifying and following the changes transient sources undergo.
The combination of mirror size, sensitivity, and smaller Galactic attenuation in the infrared, allow the James Webb Space Telescope (JWST) to reach much fainter flux levels than ever before, and make it an ideal instrument to detect sources presenting changes, be it through proper motions (e.g., outer Solar System bodies and nearby stellar and substellar objects), or in light output (e.g., supernovae, variable stars, AGN).Jansen & Windhorst (2018) proposed taking advantage of a region of the sky that can be observed by JWST any time of the year to detect transient, variable, and moving objects much fainter than is feasible from the ground even with near-future facilities like the Rubin Observatory.Examining the distribution of stars and mid-infrared emitters (bright stars, galactic cirrus), Jansen & Windhorst (2018) selected a patch of sky close to the North Ecliptic Pole (NEP) within JWST's northern continuous viewing zone (CVZ), which is a position on the sky used for spacecraft house keeping and will be frequently visited by the telescope.The position chosen in the northern CVZ is centered at (R.A., Dec.) J2000 = (17:22:47.896, +65:49:21.54),and has a particularly low number density of bright (m K ≤ 15.5 mag) sources and low Galactic extinction (E(B−V ) ≲ 0.03 mag) where such time-monitoring survey work will be feasible and efficient.This area was dubbed the JWST NEP Time Domain Field (NEP TDF) and will be a prime field for time-domain measurements during the lifetime of the telescope.It is also one of two CVZ targets of the Prime Extragalactic Areas for Reionization and Lensing Science (PEARLS) JWST/GTO program 2738 (PIs: Windhorst & Hammel) (Windhorst et al. 2023), the other being the IRAC Dark Field (Yan et al. 2023) a few degrees away.
In this paper, we report the results of a NIR imaging survey in the NEP TDF.These observations increase the depth relative to 2MASS by more than 5 magnitudes and provide baseline observations for brighter sources detected by JWST.The NIR observations cover a wider field than targeted by JWST and provide NIR-counterparts for sources outside the JWST footprint but which are contained within the fields-of-view of the other space observatories.We supplement the NIR infrared catalogs with visible data obtained from our own reduction of archival Subaru Hyper-Suprime-Cam (HSC) data taken in this part of the sky.In § 2, we describe the observations, in § 3, the reduction and calibrations, in § 4, the analyses of MMIRS and the archival HSC data, and in § 5, our conclusions.All magnitudes are quoted in the AB system (Oke & Gunn 1983) and, following Willmer (2018), we adopt the Vega to AB transformations (MMIRS-Y, MMIRS-J, MMIRS-H, MMIRS-K, MMIRS-K spec ) = (0.574, 0.891,1.333,1.836, 1.840) mag.These filters will be noted throughout the text as Y, J, H, K and K spec .The NIRCam observations as executed use a field center offset to the W, such that the full NEP TDF is covered with MMIRS only in J.

OBSERVATIONS
We imaged the JWST NEP TDF with the MMT-Magellan Infrared Imager and Spectrometer (MMIRS; McLeod et al. 2012, Chilingarian et al. 2015) on the MMT.The observations were taken in queue-mode over 31 nights between UT 2017 May 14 and UT 2019 June 21.In addition to operating MMIRS, the queue observers also recorded the sky conditions, overheads, and any anomalies for each executed group of exposures ("blocks").A log of these observations is provided in Table 10 (Appendix A).MMIRS is equipped with a single HAWAII-2RG (H2RG) detector with 2040 × 2040 light-sensitive pixels1 .The detector pixel scale of ∼0.′′ 21 pixel −1 gives a field of view of ∼6.′ 8×6.′ 8, so that 4 pointings (a 2×2 mosaic) are required to cover the 14 ′ diameter NEP TDF footprint.MMIRS allows using different gain valuesthe default low gain of 2.68 e − /DN and a high gain of 0.95 e − /DN, designed for the observations of faint sources.The detector is read out non-destructively using 32 amplifiers, generating "ramps" for each exposure ("data cubes", where the third dimension contains the individual readouts).The observations used the Y, J, H, and K filters; however, by accident, during the instrument setup early in the survey, some observations were taken using the wider spectroscopic filter K spec , others with the (short) K using the high gain value of 0.95 e − /DN rather than the default 2.68 e − /DN.Because of the different throughputs of the K and K spec filters and the different noise characteristics due to detector gain (which are discussed in greater detail in Section 4.2), we differentiate the observations according to filter and gain values; thus, the (short) K observations are represented as K 2.68 and K 0.95 .To keep the sky background and source flux within the roughly linear regime and minimize the number of saturated sources, the integration time of individual exposures ranged from ∼8 s in K to∼ 57 s in Y.The ramps for the NEP TDF observations used the LogGain/ramp 4.5s readout mode for Y and J and ramp 1.475 for H and K/K spec , generating data cubes with 14, 12, (8, 11)2 , and 7 frames in (Y, J, H, K+K spec ).For each science image, a corresponding sky-camera exposure was stored, which allows monitoring changes in the sky conditions over the duration of each observing block.At the end of each night, a series of dark frames were taken with integration times identical to those used for the science exposures.
After each exposure, the telescope was moved using a pseudo-random dither pattern to a new position at least 3 times the expected FWHM (assumed as 0. ′′ 8) away in order to mitigate both bad pixels and persistence.Table 1 shows the exposure times per dither, the total number of images in each filter combining all quadrants and the average exposure time per pixel.The accumulated total on-target science exposure time was ∼3.9, 7.1, 5.9, and 9.7 hours in Y, J, H, and K/K spec , respectively.
Figure 1 shows outlines of the MMIRS NIR observations reported here, as well as outlines of the JWST observations in the NEP TDF which are now complete, overlaid on a background g, i2, z color composite constructed from our reduction of archival HSC observations taken in 2017 for the HEROES survey (Songaila et al. 2018;Taylor et al. 2023).3. REDUCTION

Individual images
The reduction of the MMIRS imaging followed a similar procedure as 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 non-linearity using look-up tables compiled by the MMIRS team, cosmic rays were removed and cross-talk 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 procedures3 with adaptations for use with imaging data by one of the authors (CNAW) which are briefly described in Appendix B. The procedures as modified are publicly available at https://github.com/cnaw/mmirsimaging.
The sky subtraction followed the IRAF script xdimsum (Eisenhardt et al. 2001), which was specifically designed for the reduction of near-infrared 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. 4he 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 (0. ′′ 42) 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 magnitudes.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 pixels (dead pixels, hot pixels) were also flagged and their positions 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, and 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 (1.′′ 47 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 bi-weight 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 the K spec K 0.95 and K 2.68 reflect the varying conditions during the observations, taken over several nights between 2017 June and 2017 November (see Figure 10 for the variations over individual nights).In the case of K 0.95 , data were taken in a single night, for K spec over the course of two nights.Most of the K values are ≤1.′′ 5 but have a somewhat broad distribution, particularly in the case of the K 2.68 observations..68,Kspec) = (0.′′ 90, 0. ′′ 96, 0. ′′ 85, 0. ′′ 80, 1. ′′ 06, 0. ′′ 86).The distribution of K2.68 values reflects the variation of sky conditions over the several nights the data were collected.When making mosaics for the whole field, a maximum seeing cut at 2. ′′ 1 was adopted as a compromise between attaining depth and maintaining resolution; this cut mainly affects the K2.68 measurements.
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 astrometry.netsoftware of Lang et al. (2010) and the GAIA DR3 catalogue (Gaia Collaboration 2022) as reference.The improved astrometry for individual images has typical RMS residuals ∼0.′′ 03.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: Masked pixels were assigned inverse-variance values of 0 so swarp (Bertin et al. 2002) ignores these pixels when assembling mosaics, while the lowest weight values were set to 10 −9 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 swarp5 .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. 2006) for J, H, K and K spec and 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.

Mosaics
To create the full-field mosaics, we followed a similar procedure as Ashcraft et al. (2018), Ashcraft et al. (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 these 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 ∼ 2. ′′ 0, the deeper mosaic reaching ∼ 1 magnitude 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 2. ′′ 1 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 which includes the contribution from individual image weights.For the output mosaics, we adopted the pixel size of 0. ′′ 168 and center position (17:22:48.1298, +65:50:14.853,J2000) used in the Subaru/HSC imaging discussed in §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 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 on the GAIA-DR3 system, and using GAIA-DR3 stars we find that the MMIRS positions have uncertainties ≲ 0. ′′ 070.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, 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 χ 2 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, Pan-STARRS and 2MASS and only used sources that were classified as stars in GAIA and which 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), particularly 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 of 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.

Catalog completeness
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 S/N 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); 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 modelled 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 • ], [0.8 ′′ , 3.3 ′′ ] 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 b Sextractor detection threshold in terms of standard deviations above the background level. 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.

Building a combined K sample
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 2. ′′ 1 apertures in each of the K band mosaics.The K 2.68 measurements (orange distribution) show the smallest dispersion, though there is large amount of overlap between these and the measurements using the wider K spec filter (grey) 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 given quadrant.The regions where each of these settings overlap were used to estimate average offsets between the different imaging samples (Table 4 and Figure 5).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

Kspec K095 K268
Figure 4. Distribution of the dispersion of measured background fluxes using 100,000 random apertures of 12.5 pixels (2.′′ 1) calculated for the final science-grade mosaics in K2.68 (orange), K0.95 (green) and Kspec (grey).The number of apertures contributing to the histogram for each filter is proportional to area with exposure times greater than 600s and comprise 78%, 53% and 31% of the random variates for K2.68, K0.95 and Kspec respectively.The higher noise level for the K0.95 is immediately apparent; the different peaks are associated to the large contiguous regions of the mosaics with very similar exposure times.Shorter exposures shift the peaks towards higher noise values as expected; the peaks seen for K2.68 (from left to right) correspond to total exposures ∼ 9100 s, 8600 s and 3800 s, while for K0.95 the peaks correspond to ∼ 4000 s and 1000 s; for Kspec the counts peak at ∼ 6000 s.
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 of the order of ∼0.02 magnitudes.
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.

Archival Subaru Hyper-Suprime-Cam data
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 SMOKA6 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 The full sample of galaxies, with no magnitude cut, is plotted as black diamonds; green dots represent sources that remain after a 3.5σ cut.The orange diamonds represent the sub-sample of galaxies after the magnitude cuts and black plusses the same sample after a 3.5σ cut.The dashed lines represent the mean calculated using the resistant mean procedure in IDL and the horizontal dot-dashed lines the ±1 σ range for the full sample (green) and the magnitude-limited one (blue).The vertical dashed lines represent the bright and faint magnitude limits used to calculate the average offsets.
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) (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 contains data for 56786 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 is tied to the GAIA DR1 astrometric reference frame (Gaia Collaboration et al. 2016).The HSC images are calibrated into maggies7 , which are in the AB system (Oke & Gunn 1983), and the flux value corresponding to magnitude 0 is specified by the FLUXMAG0 image header keyword8 .The conversion from maggies into Jy or erg cm −2 s −1 Hz −1 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: magerr = 2.5 ln(10) area × variance(background) + flux/gain flux (2) assumes flux measurements in DN/s (data numbers per second), and 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 converting the AB magnitudes into photon numbers per cm 2 per second using where f ν 0 is the flux measured in erg cm −2 Hz −1 s −1 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 use the pivot wavelength (e.g., Bessell & Murphy (2012), eq.A13) as the characteristic wavelength.Assuming Poisson errors both for background and source photons the expression above becomes: magerr = 2.5 ln (10) [area × variance(background) + flux] × photon numbers × mirror area flux × photon numbers × mirror area where area is the total number of pixels covered by the source, and flux and 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 9 , this expression produces estimated uncertainties of 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), discussed in the next paragraph.Both flux and flux uncertainties are converted into nJy by multiplying each value by 10 31.4/2.5 .
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,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) ≤ Dec. ≤ 66.125 (67:07:30) -and contains 150,216 sources after setting the flags and for each filter {g,i2,z,NB816,NB921} base PixelFlags flag edge = f alse, (6) and {g,i2,z,NB816,NB921} base PixelFlags bad = f alse, 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 for all wide band filters and only considering 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, e.g., the HSC pipeline uses pixel radii to define the aperture photometry (SExtractor uses diameters) and the HSC pipeline adds an aperture correction when producing the final magnitudes for stars and galaxies, but 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 4.03" 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.

Merged Visible-NIR catalog
To create a merged catalog joining the PEARLS HSC and MMIRS photometry, sources were matched using a search radius of 0. ′′ 8 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 -0.′′ 070 and 0. ′′ 030 respectively -the larger search radius allows matching the larger and brighter objects.Cases where there were multiple matches were visually inspected and resolved in the cases of obvious mis-matches.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 57501 (including 56752 PEARLS HSC sources and respectively (Y, J, H, K) = (8128, 10558, 7612, 6361) of which 354 were found to be spurious.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 which 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 57467, of which 715 are NIR-detected without counterparts in the PEARLS HSC catalog.
The number of sources without PEARLS HSC detections10 but present in all NIR (Y, J, H, and K) catalogs ("zdropouts") is 118; those with H and K ("J-dropouts") is 323.The latter comprise candidate objects close to bright sources which 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 0.50" matching radius we identify 1858 counterparts, 1419 of which have MMIRS measurements in at least one filter.Color-color plot of sources matching the PEARLS HSC-MMIRS with a preliminary JWST-TDF catalog measured from the first two NEP TDF spokes.The magnitudes of the latter were measured using SExtractor with a fixed aperture determined from the F444W filter.The 174 plotted NIR sources have no HSC counterparts and the key identifies the different characteristics, i.e., sources that have at least one NIR measurement (grey dots; sources with measurements in all NIR filters (blue upward pointing triangles); sources with measurements in H and K (orange downward pointing triangles) and only in K (red circles).The 5 galaxies which match SCUBA2 sources are noted by open black squares and the 28 MMIRS-detected galaxies with no HEROES counterparts are represented by green open diamonds.
8 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 grey 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 as 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, private communication).
We used the HEROES (Taylor et al. 2023) subcatalog discussed in §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 0.5 arcseconds while 153 have no HEROES match.Of these, 12 are detected in all MMIRS bands, 93 in H, K and 20 in K only; 5 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).By 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.

Star/Galaxy separation and Number Counts
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 Figure 8. Star-galaxy separation using the total magnitude (mag auto) versus half-light radius measured from the HSC z-band photometry.Green plusses represent sources classified as stars; small blue squares represent galaxies; brown asterisks represent saturated stars or merged images where one or more sources is a star.For sources without an HSC match within 0. ′′ 80, the NIR half-light radius is used, and these are represented by black diamonds.The vertical lines represent the completeness for extended sources for the 95% level (dashed line), 80% (dot-dashed line), and 50% level (dotted line).
of MMIRS sources without PEARLS HSC half-light radius measurements, we used the corresponding NIR-band halflight 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, dot-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.
sample.The uncertainties in Table 8 assume Poisson statistics δN c (m) = ( N (m)×area−N (m)×∆area)/(∆m×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.5mag −1 per unit solid angle, where converted into N gal mag −1 deg −2 .The areas and limiting magnitudes come from the faintest magnitudes in the original references11 .Given the small area covered by the NEP TDF MMIRS imaging-(0.0635,0.0870, 0.0643, 0.0634) deg 2 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.
a Italics refer to counts affected by incompleteness estimated from the Monte Carlos simulations of Section 4.1.
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 for the 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 multi-field number counts of Keenan et al. (2010) throughout the common magnitude range.

APPENDIX
A. MMIRS OBSERVATION LOG Table 10 is the MMIRS observation log, which identifies the date of observation (UT), the modified Julian Date, the mean seeing value and dispersion measured from all images taken during the night, the quadrant(s) and filter(s), and in parenthesis the gain values.When no gain value is shown, the default 2.68 e − /DN gain was used.

B. IDL IMAGING WRAPPER
A wrapper to run the initial imaging reduction scripts is available at https://github.com/cnaw/mmirsimaging (also available at Willmer (2023)), though currently there is limited documentation on how to use them.Most of the procedures are those of Chilingarian et al. (2015) but with small modifications in the initial reduction that converts data cubes into slopes.Also in the repository is a perl script, pre reduction.pl,that will read the contents of the directory where raw images are stored, sort images according to the DARKTIME keyword value, and check whether they are science or calibration data.This script attempts to emulate the perl script described in the Chilingarian et al. (2015) paper.The pre-reduction script sorts imaging from spectroscopic data (this feature was not tested extensively) and creates a batch file that should be run under IDL.The final product of this step is a set of dark-subtracted, linearity-corrected slope images.
Figure 10.Distribution of the average seeing FWHM and dispersion measured from all images during a night presented in Table 10.The colors represent different filters as noted in the key.

Figure 1 .
Figure 1.Footprint of the imaging for different MMIRS filters in the NEP TDF using a ∼23.′ 9×22.′ 6 Subaru HSC g, i2, z color composite as backdrop, with North up and East to the left with the grid showing coordinates in sexagesimal units.The NEP TDF is centered at (R.A., Decl.)J2000 = (17:22:47.896,+65:49:21.54).The white solid line outlines the Y imaging, the green dashed line outlines J, the thick blue dash-dotted line H, and the orange dots K.Most outlines overlap each other, with the exception of J.The K outline combines observations taken with the (short) K and Kspec filters; the latter occupy the SE corner quadrant.The figure also includes the outline of the completed JWST/NIRCam observations as a thin red solid line.When the MMIRS observations were initiated, a preliminary NEP TDF pointing and set of NIRCam position angles was used.The NIRCam observations as executed use a field center offset to the W, such that the full NEP TDF is covered with MMIRS only in J.

Figure 2 .
Figure 2. Distribution of the seeing FWHM in arcsec measured for sources classified as stars in each individual image.The left panel shows the distributions for Y, J and H, and the right panel the distributions for the low and high gain K2.68, K0.95 short K filter and Kspec as indicated in the key.The median seeing values for the different configurations are (Y, J, H, K0.95, K2.68, Kspec) = (0.′′ 90, 0. ′′ 96, 0. ′′ 85, 0. ′′ 80, 1. ′′ 06, 0. ′′ 86).The distribution of K2.68 values reflects the variation of sky conditions over the several nights the data were collected.When making mosaics for the whole field, a maximum seeing cut at 2. ′′ 1 was adopted as a compromise between attaining depth and maintaining resolution; this cut mainly affects the K2.68 measurements.

Figure 3 .
Figure 3. Completeness levels for the individual bands and for the combined K+Kspec mosaic (g) for point sources (orange points and lines) and for n = 1 Sérsic profiles (blue).The vertical lines indicate (from left to right) the 95%, 80%, and 50% completeness levels.

Figure 5 .Figure 6 .
Figure 5. Distribution in right ascension and declination of sources in the regions where the three modes of K band imaging overlap.Pink circles show sources in common between all modalities (K2.68,K0.95, and Kspec).Blue points represent sources with K2.68 and K0.95 imaging; green points sources with K0.95 and Kspec.Sources detected in Kspec and K0.95only are represented by black "+'s" and green "x's" respectively.Grey diamonds represent sources with Kspec and K2.68 measurements.
Figure7.Color-color plot of sources matching the PEARLS HSC-MMIRS with a preliminary JWST-TDF catalog measured from the first two NEP TDF spokes.The magnitudes of the latter were measured using SExtractor with a fixed aperture determined from the F444W filter.The 174 plotted NIR sources have no HSC counterparts and the key identifies the different characteristics, i.e., sources that have at least one NIR measurement (grey dots; sources with measurements in all NIR filters (blue upward pointing triangles); sources with measurements in H and K (orange downward pointing triangles) and only in K (red circles).The 5 galaxies which match SCUBA2 sources are noted by open black squares and the 28 MMIRS-detected galaxies with no HEROES counterparts are represented by green open diamonds.

Table 1 .
Summary of MMT/MMIRS Observations of the NEP TDF NE a,b , NW a , SE b , SW b
Filter m

Table 4 .
K magnitude differencesFilter and gain setup mean difference uncertainty a σ

Table 7 .
Magnitude differences between HEROES and this work.

Table 9 .
Depth and area of literature number counts.These observations provide the first epoch deep NIR photometry in this part of the sky and will continue being used in the search for transient and variable objects in the NEP TDF.This region of the sky has already amassed in very few years an extensive coverage in wavelength space, ranging from X-rays to radio frequencies.The mosaics and catalogs are publicly available at https://sites.google.com/view/jwstpearls/data-productsand on Zenodo under an open-source Creative Commons Attribution license: doi:10.5281/zenodo.7934393Wewouldlike to dedicate this paper to the memory of our PEARLS colleague Mario Nonino and acknowledge his suggestions in preparing images, catalogs and documentation for the public release.We respectfully acknowledge that the University of Arizona and Arizona State University are on the land and territories of Indigenous peoples.Today, Arizona is home to 22 federally recognized tribes, with Tucson being home to the O'odham and the Yaqui and Tempe home to the Akimel O'odham (Pima) and Pee Posh (Maricopa).Committed to diversity and inclusion, the University of Arizona strives to build sustainable relationships with sovereign Native Nations and Indigenous communities through education offerings, partnerships, and community service.We wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the Kānaka Maoli community.We are most fortunate to have the opportunity to conduct observations from this mountain.We thank Brian McLeod (MMIRS PI) and Igor Chilingarian for their assistance and insight in preparing and reducing the MMIRS data and for providing the filter throughputs used in the photometric calibration.We also thank Dallan Porter and Sean Moran for their assistance in implementing more MMTO data retrieval options and Catherine Grier for proving a script that enabled calculating exposure maps using swarp.We thank the anonymous referee for comments and suggestions that allowed improving the presentation.CNAW and MJR acknowledge support from the NIRCam Development Contract NAS5-02105 from NASA Goddard Space Flight Center to the University of Arizona.CNAW and RAJ gratefully acknowledge funding from the HST-GO-15278.008-Agrant.RAW, SHC, and RAJ acknowledge support from NASA JWST Interdisciplinary Scientist grants NAG5-12460, NNX14AN10G and 80NSSC18K0200 from GSFC.JFB was supported by NSF grant No. PHY-2012955.IS acknowledges support from STFC (ST/T000244/1).MI acknowledges the support from the National Research Foundation of Korea (NRF) grants, No. 2020R1A2C3011091, and No. 2021M3F7A1084525, funded by the Korean government (MSIT).M.H. acknowledges the support from the Korea Astronomy and Space Science Institute grant funded by the Korean government (MSIT) (No. 2022183005) and the support from the Global Ph.D. Fellowship Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2013H1A2A1033110).We acknowledge the SMOKA archive of the National Astronomical Observatory of Japan from which Subaru Hyper-Suprime-Cam data were retrieved.This research has made use of the NASA/IPAC Extragalactic Database, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.