JWST PEARLS. Prime Extragalactic Areas for Reionization and Lensing Science: Project Overview and First Results

PEARLS targets for First Light and Reionization studies include high-redshift Lyα galaxies and protoclusters. In light of observations with HST WFC3 over the past 13 yr, PEARLS will also image several rich galaxy clusters that boost the signal of faint, high-redshift objects via strong gravitational lensing. Blank-field surveys will contribute to the First Light theme via number counts. To study the Assembly of Galaxies, we will observe galaxies up to the highest redshifts, and lowest masses and luminosities, in different environments. We will also investigate SMBH growth by observing high-redshift galaxies having an active nucleus: quasars and radio galaxies. The blank fields at high Ecliptic latitude will contribute time-domain information. PEARLS will also study VV 191, a nearby, overlapping galaxy pair, to provide a benchmark dust-attenuation profile for studying high-redshift, dusty environments. Table 1 summarizes the PEARLS fields observed thus far (as of 2022 July 31), and Table 2 the PEARLS fields to be observed subsequently. The 112.3 calendar hours allocated to PEARLS include 2.3 hr from H. Hammel. In all, PEARLS will image 16 NIRCam and four NIRISS fields in up to eight filters to AB ≲ 28.5–29 mag and will cover ∼165.66 arcmin² or 0.046 deg², equivalent to ∼34 HUDF/XDF fields (e.g., Beckwith et al. 2006; Koekemoer et al. 2013).

Table 1. PEARLS Targets with NIRCam Images Taken as of 2022 July: Depth from ETC, SourceExtractor and Galaxy Counts

Instr.+Filters	R.A. (J2000)	Decl.	Obs. Date	Visit	Area	SCeff	Net	Net texp (s)
Target	(h m s.sss)	(° ' '')	YYYY-MM-DD	No.	(') × (')	(%)	(hr)	5σ Point-source AB Limit
NIRCam Broadband:								F090W	F115W	F150W	F200W	F277W	F356W	F410M	F444W
VV191-Backlit	13 48 22.0990	+25 40 40.01	2022-07-02	341.1	2.15 × 4.30	32.0	0.52	0934	⋯	0934	⋯	⋯	0934	⋯	0934
PSF-FWHM (")								0.066	⋯	0.068	⋯	⋯	0.164	⋯	0.163
ETC 5$\sigma$ AB lim								27.62	⋯	28.01	⋯	⋯	28.00	⋯	27.59
Cat 5$\sigma$ AB lim								27.88	⋯	28.24	⋯	⋯	29.01	⋯	28.81
Counts 80% compl								27.3	⋯	27.6	⋯	⋯	28.5	⋯	28.3
${\rm{\Delta }}$ AB lim(80%-ETC)								−0.3	⋯	−0.4	⋯	⋯	+0.5	⋯	+0.7
IRAC-Dark-ep1	17 40 08.5352	+68 58 27.00	2022-07-08	121.1	2.15 × 4.30	61.2	1.76	⋯	⋯	3157	3157	⋯	3157	⋯	3157
PSF-FWHM (")								⋯	⋯	0.063	0.075	⋯	0.166	⋯	0.164
ETC 5$\sigma$ AB lim								⋯	⋯	28.96	29.13	⋯	28.81	⋯	28.41
Cat 5$\sigma$ AB lim								⋯	⋯	28.75	28.93	⋯	29.67	⋯	29.43
Counts 80% compl								⋯	⋯	28.1	28.2	⋯	29.0	⋯	29.0
${\rm{\Delta }}$AB lim(80%-ETC)								⋯	⋯	−0.9	−0.9	⋯	+0.2	⋯	+0.6
El-Gordo	01 02 55.4000	−49 15 38.00	2022-07-29	241.1	2.15 × 4.30	57.3	2.50	2491	2491	1890	2104	2104	1890	2491	2491
PSF-FWHM (")								0.062	0.057	0.062	0.074	0.119	0.171	0.153	0.160
ETC 5$\sigma$ AB lim								28.43	28.61	28.60	28.88	28.58	28.55	28.03	28.32
Cat 5$\sigma$ AB lim								28.57	28.69	28.57	28.87	29.43	29.55	29.06	29.30
Counts 80% compl								27.9	27.9	27.7	28.1	28.8	28.9	28.1	28.9
${\rm{\Delta }}$AB lim(80%-ETC)								−0.5	−0.7	−0.9	−0.8	+0.2	+0.4	+0.1	+0.6
NIRCam Medium-band:								F115W	F150W	F182M	F210M	F300M	F335M	F360M	F444W
TNJ1338-1942	13 38 26.1000	−19 42 28.00	2022-07-01	361.1	2.15 × 4.30	37.9	0.86	⋯	1031	1031	1031	1031	1031	1031	⋯
PSF-FWHM (")								⋯	0.064	0.071	0.079	0.125	0.169	0.160	⋯
ETC 5$\sigma$ AB lim								⋯	27.86	27.51	27.30	27.15	27.25	27.28	⋯
Cat 5$\sigma$ AB lim								⋯	27.7	27.4	27.2	28.35	28.25	28.16	⋯
Counts 80% compl								⋯	27.1	26.6	26.4	27.8	27.4	27.3	⋯
${\rm{\Delta }}$AB lim(80%-ETC)								⋯	−0.8	−0.9	−0.9	+0.7	+0.2	+0.0	⋯

Note. For each object, line 1 lists the J2000 (R.A., decl.) tangent point to which the images were drizzled, the observing date, the APT visit number, the area covered, the net exposure time per filter and the net total hours per visit, as well as the visit's spacecraft efficiency. Line 2 lists for each filter the stellar PSF-FWHM in arcsec as measured from unsaturated stars in the drizzled images. Line 3 lists the 5σ point source sensitivity in AB-mag predicted by the prelaunch ETC for the net integration time on the first line of each target. NIRCam ETC Parameters used were aperture radii r = 008 for SW and r = 016 for LW, and sky annuli r = 03–099 for SW and r = 06–198 for LW. Line 4 lists the 5σ detection limit derived from the AB level in Figures 4–6 where the median SourceExtractor catalog flux error is 0.20 mag. Line 5 indicates the AB level in Figures 4–8 where the galaxy counts are ∼80% complete compared to a power-law extrapolation. Line 6 indicates the difference between the 80% galaxy count and predicted ETC 5σ point-source completeness limits in AB-mag. All PEARLS NIRCam images have a zero-point of 28.0865 to convert the flux (in MJy sr⁻¹) in each drizzled 00300 pixel to AB-mag.

Download table as:  ASCIITypeset image

PEARLS will obtain data over at least 13 independent lines of sight more than three degrees apart from each other and is therefore more robust against cosmic variance (CV) at AB ≲ 28.5 mag than programs that image only a few areas (e.g., Somerville et al. 2004; Driver & Robotham 2010; Windhorst et al. 2022). Figure 1 compares the area and depth covered by PEARLS to other JWST Cycle 1 surveys. While not as deep or wide as other contiguous JWST Cycle 1 surveys, PEARLS covers more fields across the sky to decrease the effects of CV. The expected CV for PEARLS fields can be found with the calculator ⁵⁸ of Driver & Robotham (2010) based on the areas covered and sensitivity limits in Tables 1 and 2. To AB ≲ 28.5 mag, the PEARLS fields sample a typical redshift range of z ≃ 0.3–8 with a median redshift of z ≃ 1–2 (see Section 4.5 and Appendix B.2). The NIRCam field of view (FOV) covers ∼0.0026 deg² (Section 3.1 or ∼1.1 × 2.2 Mpc), over which its CV is then predicted to be 30%. For the two PEARLS fields with galaxy counts presented here in Section 4.5, CV is expected to be ≲9%. At the end of JWST Cycle 1, large JWST NIRCam parallel programs like PANORAMIC (PID 2514; C. Williams PI) may push CV of the sampled objects to ∼1%–2%.

Zoom In Zoom Out Reset image size

In four of our NIRCam pointings, coordinated NIRISS grism and imaging parallels will cover a significant portion of our NIRCam images (Table 2), while UV-optical images are available from HST WFC3+ACS. The coordinated NIRISS parallels will be used for both object characterization and redshifts, and to expand the area and time-baseline of time-domain studies. The coordinated parallel observations are critical to obtain imaging and grism data that is as homogeneous as possible, over as large an area as possible, and in the least amount of time feasible with JWST.

Two of the PEARLS blank fields and two galaxy cluster fields will be observed more than once. This time-domain component will allow us to find and study Galactic brown dwarfs via high proper motion or atmospheric variability, variable active galactic nuclei (AGN), high-redshift supernovae, and any time-varying objects seen behind lensing clusters, including possible cluster caustic transits. The PEARLS time-domain data at the North Ecliptic Pole (NEP) may also reveal some faint moving objects at high Ecliptic latitude in our outer solar system.

To encourage immediate use of JWST data by the community and follow-up proposals by JWST Cycle 2 GO proposers, we will make the first epoch of our JWST NEP Time-Domain Field (TDF) public immediately (#112.* in Table 2). The other three JWST NEP TDF epochs will be released together with the v1 data products as soon as we have these. Also public right away are the Cycle 25 (R. Jansen et al. 2022, in preparation), 28, and 29 HST WFC3/UVIS F275W and ACS/WFC F435W+F606W observations of the NEP TDF and other ancillary data across the electromagnetic spectrum, as these become available and their data reduction is completed. These include 600 ks NuSTAR 3–24 keV images (Zhao et al. 2021; F. Civano et al. 2022, in preparation), 900 ks of Chandra ACIS 0.2–10 keV images (P. Maksym et al. 2022, in preparation), 31 hr of JCMT/SCUBA-2 plus 66 hr of SMA data at 0.85 mm, as well as 70 hr of VLA 3 GHz A+B-array images (M. Hyun et al. 2022, in preparation), 147 hr of VLBA 4.5 GHz data at milliarcsec (mas) resolution to sub-micro-Jy levels, and 75 hr of LOFAR 150 MHz images including LOFAR VLBI. The presence of a S_{3 GHz} ≃ 239 mJy quasar at z = 1.4429 in the JWST NEP TDF that is unresolved at VLBI milliarcsecond resolution is used as phase calibrator to provide high resolution VLA/VLBA and LOFAR/VLBI images of very high dynamic range in the NEP TDF. The NEP TDF database also includes multi-epoch Large Binocular Telescope/LBC + Subaru/HSC Ugiz images to AB ≲ 26.0 mag, Gran Telescopio Canarias/HiPERCAM ugriz images to AB ≲ 27 mag, Multiple Mirror Telescope/MMIRS images to YJHK ≲ 24–23 mag, and MMT/Binospec and MMIRS spectra to 22–24 mag (C. Willmer et al. 2022, in preparation), plus JPAS 56-narrow-band spectrophotometry to provide confirmation of the astrometric, photometric, and spectroscopic calibration of our JWST NIRCam+NIRISS observations.

PEARLS target selection began in the early 2010s, when it became clear that JWST had a viable path toward launch and that it could perform as designed. The largest blank field is in the JWST continuous viewing zone (CVZ) near the NEP. This NEP TDF has the best combination of low foreground extinction and absence of AB ≲ 16 mag stars (Jansen & Windhorst 2018). A second blank field is within the IRAC Dark Field, which is a Spitzer/IRAC calibration field near the NEP observed repeatedly for over 15 yr. These historical light curves offer several examples of what might be high-redshift, dusty supernovae in ultraluminous infrared galaxies selected by Herschel (Yan et al. 2018). Figures 1(a) and 2(a) of Jansen & Windhorst (2018) give a layout of the JWST CVZ in the NEP, where the IRAC Dark Field (IDF) is ∼356 NE of the TDF. Our Figure 2 shows the first-epoch NIRCam observation of the JWST IDF (hereafter the JWIDF). The final blank field is in the WFC3 ERS area (Windhorst et al. 2011), which is in the northern part of the GOODS-South area.

Zoom In Zoom Out Reset image size

PEARLS gravitational-lensing clusters were selected to have high mass and central compactness or to have apparent double-cluster nature. The latter could result in higher transverse motions and therefore makes caustic transits more likely. Possible transiting sources include distant, luminous single stars, double stars, and possibly stellar-mass black hole accretion disks (e.g., Miralda-Escude 1991). All of our selected clusters show gravitationally lensed arcs, and all have lensing magnification maps produced with multiple independent lensing models, which will be refined with the JWST data. Other lensing clusters were similarly selected because of their high mass and high central compactness, and their lower IntraCluster Light (ICL) content, which could make it easier to detect caustic transits with less microlensing by foreground stars in the cluster ICL (e.g., Windhorst et al. 2018). The PEARLS lensing clusters are:

• 1.  
  The HFF (Lotz et al. 2017) cluster MACS J0416.1−2403 at z ≃ 0.397. This field will be covered by three JWST epochs about six months apart to maximize the chance of seeing caustic transits at z ≥ 6 (e.g., Windhorst et al. 2018; Welch et al. 2022a, 2022b). A number of plausible caustic transits at z ≃ 0.9–1.5 have already been observed for this cluster by HST (Dai et al. 2018; Kelly et al. 2018; Dai 2021).
• 2.  
  Abell 2744 at z ≃ 0.31 and MACS J1149.5+2223 at z ≃ 0.54. These are likewise HFF (Lotz et al. 2017) clusters. They will have additional GTO observations by Wilmott & the NIRISS GTO team (2022), and by Stiavelli (2022) to look for variable objects in or behind these clusters, and by the GLASS team (PID 1324; PI: Treu et al. 2022). This allows us to monitor potential high-redshift caustic transits on timescales longer than a year.
• 3.  
  The cluster known as El Gordo at z ≃ 0.87 (Menanteau et al. 2012; Zitrin et al. 2013; Cerny et al. 2018; Diego et al. 2020; Caputi et al. 2021). This cluster was selected because of its enormous mass (Menanteau et al. 2012, ∼2 × 10¹⁵ M_⊙), its elongation to maximize the probability of caustic transits (Windhorst et al. 2018), and its rich collection of distant lensed source candidates. The field includes a background galaxy grouping at z ≃ 4.3 (Caputi et al. 2021) that is lensed by the cluster. Figure 3 shows the module around El Gordo that did not cover the central part of the cluster (hereafter referred to as the "noncluster" module; see Section 4.5).
• 4.  
  PLCK G165.7+67.0 (G165) is a double cluster at z ≃ 0.35 selected by its FIR colors, and not by the Sunyaev–Zeldovich effect (e.g., Cañameras et al. 2015; Harrington et al. 2016; Planck Collaboration et al. 2016). Many caustic-crossing arcs are detected, which are well-suited to transient science. One example is a strongly lensed red and dusty submillimeter galaxy detected in the HST imaging, whose counter image appears in the LBT/LUCI+ARGOS laser-guided AO K-band images (e.g., Frye et al. 2019; Rabien et al. 2019). Spectral energy distributions (SEDs) fit to the optical–near-IR images yield photometric redshift estimates for some image families, which constrain the lens model (Pascale et al. 2022). The model confirms the bi-modal mass distribution of this ongoing merger that is only a low-luminosity X-ray source. The JWST observations aim to constrain the dynamical state of this cluster and detect a significant number of lensed background sources.
• 5.  
  The Clio cluster at z ≃ 0.42 from the GAMA survey (Driver & Robotham 2010; Alpaslan et al. 2012) is a massive, compact cluster selected to have significant potential for lensing background sources. A ground-based VLT image (Griffiths et al. 2018) already showed strongly lensed arcs and a lower-than-average amount of IntraCluster Light. This is attractive because low-mass IntraCluster Medium stars can significantly lengthen caustic-transit times (e.g., Diego et al. 2018; Windhorst et al. 2018) and complicate their lensing analysis.
• 6.  
  The cluster RXC J1212+27 = A1489 at z ≃ 0.35. This cluster was chosen because of strong gravitational lensing, using the automated implementation (Zitrin et al. 2012) of the light-traces-mass method (e.g., Zitrin & Broadhurst 2016; Zitrin 2017). HST images showed a significant number of lensed sources that resulted in a good lensing model (Zitrin et al. 2020).

Zoom In Zoom Out Reset image size

The first public JWST images (Pontoppidan et al. 2022) released starting 2022 July 12 have already inspired a number of further studies. Relevant for PEARLS, a possible caustic transit candidate has been suggested at z ≃ 3 in some of the public JWST images of the cluster SMACS0723 (e.g., Chen et al. 2022), for which mass models were made by (e.g., Pascale et al. 2022), and in which also a significant number of red spirals were identified (e.g., Ferreira et al. 2022; Fudamoto et al. 2022). Indeed, some very high redshift candidates were already suggested in some of the very first JWST ERS images (Adams et al. 2022; Finkelstein et al. 2022). The PEARLS high redshift protoclusters are:

• 1.  
  PHz G191.24+62.04 (G191) is a protocluster candidate at z = 2.55 with one of the highest star formation rates (SFR ≃ 23,000 M_⊙ yr⁻¹) in the parent sample of Planck high-z sources (PHz; Planck Collaboration et al. 2016). G191 hosts an overdensity of red Spitzer sources (Planck Collaboration et al. 2015; Martinache et al. 2018), containing ∼14 objects arcmin⁻² with IRAC 3.6–4.5 μm colors ≳–0.1 mag. Two of the Herschel sources have spectroscopic redshifts and a large estimated SFRs ≃1000–1500 M_⊙ yr⁻¹ (Polletta et al. 2022), i.e., high enough that they present challenges for theoretical models (Granato et al. 2015; Lim et al. 2021; Gouin et al. 2022). The JWST observations will constrain the stellar mass assembly and fueling mechanism (e.g., major mergers, cold accretion) occurring in this highly star-forming high-z structure.
• 2.  
  TNJ1338−1942 is a protocluster at z = 4.1 that was discovered with the VLT as 60 Lyα-emitters near a luminous, steep-spectrum radio source (e.g., De Breuck et al. 1999, 2000; Venemans et al. 2002; Miley et al. 2004; Intema et al. 2006; Saito et al. 2015). The radio source's AGN activity and outflow will be studied by a JWST Cycle 1 GO program (PID 1964, PI R. Overzier). PEARLS has imaged the field in the five NIRCam medium-band filters that best straddle the Balmer/4000 Å break at z = 4.1 to help delineate the ages of ∼30 of the Lyα-emitters. Figure 4 shows part of the NIRCam image around TNJ1338−1942. To maximize the scientific return on TNJ1338−1942, the analysis of the PEARLS and GO data will be coordinated.

Zoom In Zoom Out Reset image size

In addition to our above two protocluster targets, the z ≃ 4.3 group of galaxies behind El Gordo (Caputi et al. 2021) may also turn out to be a protocluster candidate. Additional PEARLS targets are:

• 1.  
  Two QSOs, QSO 1425+3254 (or NDWFS J142516.3+325409 at z = 5.85; e.g., Mechtley et al. 2012), and QSO J0005−0006 (or SDSS J000552.35−000655.6 at z = 5.86). The first has a number of possible z ≃ 6 companions (e.g., Marshall et al. 2020, 2021). PEARLS IFU observations will address whether these form a group around the QSO. QSO J0005−0006 was selected because it lacks both hot and cold dust (Wang et al. 2008; Jiang et al. 2010). It therefore represents a rare subpopulation of dust-free high-z quasars.
• 2.  
  The VV 191 system (Figure 5) consists of a foreground spiral galaxy with an unassociated elliptical galaxy behind it (e.g., Keel et al. 2013). Light from the elliptical suffers extinction from dust in the spiral. PEARLS NIRCam imaging maps the extinction and determine its wavelength dependence (Keel et al. 2022).

Zoom In Zoom Out Reset image size

2.3.1. JWST Observation Planning of PEARLS Targets

2.3.2. PEARLS' Primary NIRCam Observations and Areas

2.3.3. PEARLS' Coordinated NIRISS Parallel Observations and Areas

Calibration of PEARLS data obtained as of 2022 July used the calibration files on the STSCI JWST website as of 2022 July 12. All data were processed with the standard STScI pipeline CALWEBB, ⁶⁶ which comes in three stages: (1) detector-level corrections to the raw individual exposures to produce count-rate images from the nondestructive readouts ("ramps"); (2) photometric and astrometric calibration of the individual exposures; and (3) drizzling the calibrated and distortion-corrected images into mosaics. The NIRCam pipeline CALWEBB was used to process all our images. The Calibration Reference Data System (CRDS) provides the latest reference files that we used to calibrate our data. ⁶⁷ CRDS version 11.13.1 was used for all images.

Our NIRCam TNJ1338, VV191, and JWIDF images taken in early 2022 July were initially reduced with Pipeline version v1.6.1.dev2+g408c711 and context file jwst_0916.pmap_filters, which contained the prelaunch ZP values available then. Our more recent El Gordo images of 2022 July 29 were reduced with Pipeline version 1.6.2 in early 2022 August using context file jwst_0942.pmap_filters, which implemented Rigby et al. (2022) in-flight ZPs affecting all NIRCam filters. We refer to these calibrations and their resulting mosaics as version v0.5. When more accurate on-orbit NIRCam flat-field and ZP calibrations became available in early 2022 October, we reprocessed our PEARLS images into a version v1 with context file jwst_0995.pmap_filters and Pipeline version 1.7.2. We will make v1 of the NEP TDF available to the community as soon as its catalogs are completed and verified (R. Jansen et al. 2022, in preparation). Further details of the 2022 October calibration improvements are given in Section 3.3 and Appendix B.1.

Performance of the NIRCam detectors is relevant to depth, calibration quality, and accuracy of the sky-SB values discussed in Section 4–5. SW and LW module characteristics relevant for PEARLS are: ⁶⁸

• 1.  
  Average NIRCam read-noise values are ∼16.2 and ∼13.5 e⁻ pixel⁻¹ using correlated double-sampling.
• 2.  
  Average dark-current values in typical exposures are very low: ∼0.0019 ± 0.002 and ∼0.027 ± 0.005 e⁻ s⁻¹.
• 3.  
  Average detector gains are ∼2.05 and ∼1.82 e⁻ ADU⁻¹, respectively.
• 4.  
  Persistence of charge from a previous equal-length exposure is ≲0.01% of the original charge detected in the previous image.

The first step in mosaicing was to anchor all individual frames into the Gaia DR3 reference frame (Gaia Collaboration et al. 2022). This step used catalogs based on recent, deep ground-based or HST images already referenced to Gaia DR3. Both stars and galaxies were used for the correction with star positions corrected for proper motion from the epoch of the reference image to the JWST observing epoch (Table 1). Coordinate differences between the catalog and JWST frame were measured, and each frame's center position and position angle were adjusted to minimize the differences. Typical adjustments were ≲10 mas, and the final uncertainty in each frame's position is about 2–3 mas rms. We used AstroDrizzle ⁶⁹ (Koekemoer et al. 2013; Avila et al. 2015) to drizzle the NIRCam images as calibrated in the CALWEBB pipeline Stages 1 and 2 into two mosaics for each field. Pixel sizes used were 00300 and 00600 for SW and LW, respectively. For the LW module, we also provide 00300 pixel⁻¹ mosaics to facilitate aligned analysis. The higher resolution of the former samples the NIRCam PSFs better in the SW channel, while the bigger pixels of the latter have better SB sensitivity for the generally short PEARLS exposures. All NIRcam images were drizzled after we removed wisps as well as possible, applied a 1/f correction, flagged the snowballs before object detection, and subtracted a surface-fit to the sky-SB between all the detected objects. Details on these aspects are given below.

1. NIRCam 1/f Noise Pattern Removal: NIRCam images are read out simultaneously in four vertical 510 × 2040 pixel strips. Differing gains or zero-points in the four amplifiers causes banding or striping across the images as they are read out. Rest (2014) presented a mathematical method to separate this 1/f noise from the random read noise ⁷⁰ and derived the time dependence of the 1/f noise. Rest (2014) found that the 1/f noise strongly correlates between the amplifiers of a given detector because it is caused by a common reference voltage. Rest (2014) also found that the 1/f noise has reproducible spatial structure at the 10%–20% level down to spatial scales of tens of pixels, and this structure does not seem to change on timescales of months. Schlawin et al. (2020) described methods to reduce 1/f-noise patterns in the highly demanding NIRCam grism time-series observations of exoplanets. As a bonus, their spatial background subtraction also efficiently removes many random detector defects, including preamplifier offsets, amplifier discontinuities, and even–odd column offsets. The NIRCam 1/f noise can also be reduced by subtracting the values from background pixels or reference pixels that are read closely in time. Schlawin et al. (2020) removed the 1/f noise as a step in the pipeline after the superbias correction. ⁷¹ Hilbert et al. (2016) presented a method to subtract 1/f noise from NIRCam integrations before averaging the data to produce superbias maps. This method produced superbias images with significantly lower noise levels than images produced using the more traditional approach. Wilmott & the NIRISS GTO team (2022) provided a more recent code to remove 1/f noise that runs on the calibrated files. ⁷² Bagley et al. (2022) presented a code to remove both the detector-level offsets in the SW modules and the 1/f noise patterns. This code was produced for the CEERS project and is part of their SDR1 release. ⁷³

We used both the Willott code and the ProFound code (Robotham et al. 2017, 2018) to remove the 1/f noise patterns. Together with the low-level pedestal removal between the SW detectors below, the ProFound-package resulted in images that are mostly visually flat without major row-based artifacts. We visually verified that the 1/f noise-removal parameter settings did not introduce new artifacts in the final images. Further details are given in Appendix A, which compares the results from both 1/f-noise removal methods and verifies that these do not noticeably affect object photometry at S/N levels ≳5σ.

2. SW Detector-Level Offsets: the eight detectors in the SW modules A and B have detector-level offsets that are a combination of additive and multiplicative corrections, although most of the effect seems to be additive. Except in very crowded fields, these detector-level offsets are relatively easy to remove early in our data reduction workflow with the ProFound-based code (Robotham et al. 2018), as implemented in our previous HST WFC3/IR work (e.g., Windhorst et al.2022). In short, with ProFound we created a new fits extension SKY that removed bad pixels and real objects detected to AB ≲ 28.5 mag (Section 4) and that interpolated the local sky-SB plus its local rms noise underneath each object. With a number of images in the NIRCam broadband filters now available, we created a low-frequency supersky image "SKY_SUPER" as a clipped mean over these SKY images. The SKY images were then subtracted from each science SCI image using:

$$\begin{eqnarray}&&({SCI}-{SKY})={SCI}-(M\times {SKY}\_{SUPER}+P),\end{eqnarray} \tag{ 1 }$$

where M and P describe the linear model and pedestal, respectively, of the SKY image pixels made for each detector and filter combination from these SKY_SUPER frames. In this process, we also determined the lowest object-free sky-SB in each detector following Section 4.2 of Windhorst et al. (2022), which is used for sky-SB estimates in our 13 PEARLS filters in Section 5. By design, our medium-deep PEARLS images come from relatively short integrations (Tables 1–2), and we adjusted our signal-to-noise-ratio criteria for object detection (Section 4) to achieve uniform detection above any residual 1/f-noise patterns. Further details of the 1/f and pedestal removal procedures are given in Appendix A.

3. NIRCam SW Detector Wisps: the so-called "wisps" are caused by straylight hitting a secondary mirror support bar and then being reflected into the main light path. Only four NIRCam detectors are affected (SW A3, A4, B3, and B4) and mainly in the F090W, F115W, F150W, F182M, 200W, and F210M filters. Wisp positions are fixed on the telescope and each detector, and therefore a sky-flat made in a given filter from available images is able to subtract much of the wisp pattern locally. We used our NIRCam images to improve available wisp templates to be subtracted from the images that show visible wisps. We then subtracted the wisp template from each image that best matches the wisp amplitude. This amplitude needs to be high enough that no significant positive wisp signal is left but not so high that a negative hole is created in the local sky-SB. Some wisp residuals are left in some of the images, and we masked these areas in the faint-object-detection phase (Section 4.1) and in the sky-SB analysis (Section 5). Steady collection of NIRCam images over time is expected to improve the wisp templates to allow more accurate subtraction in future image reductions.

4. NIRCam Snowballs: NIRCam detectors can show large artifacts that resemble "snowballs" when an energetic cosmic-ray impact occurs. ⁷⁴ The vast majority of CR hits impact only a few detector pixels, but snowballs can affect several hundred pixels. They occur approximately at a rate of 50 snowballs per 1000 s of exposure time in each NIRCam detector. Our four-point dithers enable the pipeline to remove many of the snowball artifacts, but this process is not perfect. E.g., Figure 2 shows a few dim green rings left over from snowballs that were not fully removed in the drizzling process of the JWIDF F200W images. Because the drizzle weight-maps give these pixels a lower weight, we can still derive catalogs from this image. Nevertheless, we vet these catalogs carefully and mask out areas that are visibly affected by residual snowballs when doing faint object counts as in Section 4. That is, any large remaining defects that generate a local excess in the object catalogs are masked before running our final catalogs.

5. NIRCam PSFs: JWST NIRCam Point Spread Functions are detailed on the STSCI website. ⁷⁵ We generated JWST PSFs with the WebbPSF tool. ⁷⁶ Brighter star images (Section 4.1) have PSF FWHM values consistent with a diffraction limited telescope at 1.1 μm wavelength (Rigby et al. 2022), i.e., much better than the JWST Observatory requirement of a diffraction limit that was designed and kept during development at 2.0 μm (Section 4.4).

6. NIRCam Flat Fields: the accuracy of the NIRCam flat-fields is better than 7% rms (B. Sunnquist 2022, private communication) and has improved to ∼2% with the release of the flat fields captured in jwst_0952.pmap_filters and most recently jwst_0995.pmap_filters. We processed our images with this latest context file to estimate the sky-SB values in each exposure and assess their quality in Section 5. Table 3 summarizes the parameters that characterize the 0.9–4.5 μm galaxy counts and integrated galaxy light, which are needed in Tables 4–6 (see Section 5). Tables 4–5 gives the predicted and observed sky-SB for the three PEARLS targets observed in 4–8 broadband NIRCam filters, and Table 6 gives the values for TNJ1338 and its five medium-band filters.

7. NIRCam Zero-points: Measured zero-points (ZPs) for all JWST+NIRCam+Filters are on the STScI website. ⁷⁷ The JWST Mission Requirement is that the absolute ZPs of the imaging filters are known to better than 5%, ⁷⁸ and they are stated to be good to ∼4% or better for Cycle 1 (Rigby et al. 2022). The "Throughput" column in their Tables 4–5 lists the in-flight zero-points in units of DN/s/nJy.

Results based on the initial v0.5 calibrations are recorded in our first submission of this paper on https://arxiv.org/abs/2209.04119. Boyer et al. (2022) ⁷⁹ and their cited URLs analyzed new standard star observations and updated the NIRCam ZPs. Our v1 results below use this latest calibration, which more accurately corrects for ZP variations between each of the 10 NIRCam detectors. Typical ZP changes for individual detectors were ≲10%–20%. The new NIRCam F356W and F444W zero-points produce photometry in the JWIDF field consistent with the deepest Spitzer images available to within 2.6%–2.9% (Yan et al. 2018). The v1 calibration also tightened the dispersion between the bright end of the NIRCam galaxy counts and the faint end of the ground-based, HST, and Spitzer galaxy counts. The uncertainty in our estimates of the Integrated Galaxy Light in Section 4 went down, the rms variation between the sky-SB measurements decreased, and our limits on diffuse light in Section 5 improved. Further details of the calibration improvements are given in Appendix B.1.

For context, all HST ZPs were defined in units of AB-mag for a count rate of 1.000 e⁻ pixel⁻¹ s⁻¹. This definition permitted monitoring of the ZPs' wavelength and time dependence over many decades (e.g., Calamida et al. 2022; Windhorst et al. 2011, 2022). (For a summary, see Section 4 and references therein of Windhorst et al. 2022.) However, all JWST ZPs instead have been defined in units of MJy sr⁻¹. Conversion between the two sets of units can be made using the footnotes of Section 1 but also requires knowledge of the drizzled pixel scale in the case of JWST. For completeness, we therefore list both the drizzled pixel scale and the resulting equivalent ZP in AB-mag for our PEARLS NIRCam images in the footnotes of Table 1 and in Appendix B.1.

Note that given this different JWST ZP definition, the equivalent JWST ZPs in AB-mag are no longer wavelength dependent—unlike the case of HST—but only depend on the image pixel scale, which therefore should always be stated. To leave no further ambiguity, for our basic drizzled pixel scale of 00300 pixel⁻¹, the JWST ZP for 1.000 MJy sr⁻¹ converted to AB-mag would thus be the following same value for every wavelength:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ZP} & = & 8.900-2.5\mathrm{log}[{10}^{6}/{\left((360\times 3600)/(2\times \pi \times 0.0300)\right)}^{2}]\\ & = & 28.0865\,\mathrm{AB} \mbox{-} \mathrm{mag}\,\mathrm{per}\,\mathrm{pixel}.\end{array}\end{eqnarray} \tag{ 2 }$$

The constant 28.0865 in Equation (2) will be valid at all wavelengths for all our images at 00300 pixel⁻¹ if the flux calibration is correct. With the v1 calibration, this appears to be the case to within 3%–4% (see Appendix B.1.)

8. NIRCam Straylight Levels: JWST has an open-architecture Optical Telescope Element (OTE), and it will have more straylight (SL) than a closed-tube design such as HST or Spitzer. The JWST Project designed the JWST sunshield and baffles to minimize SL with expected levels ≲20%–40% (worst case) of the Zodiacal SB in a given direction. Bright near-IR sources like the Zodiacal cloud, Galactic Center and Galactic plane—the brightest NIR sources in the sky other than the Sun, Earth, and Moon—can add SL that scatters off dust accumulated on the primary mirror into the telescope FOV. Estimates of the SL levels have been made by, e.g., Lightsey (2016) and are incorporated in extensions of the JWST ETC predictions. Tables 4–6 show the predicted SL levels for our PEARLS targets. The predicted SL levels are modest at 0.9–3.5 μm, but they increase at wavelengths longer than 4.5 μm due to the increased thermal foreground from the telescope and Zodiacal belt. Further details on the adopted SL levels are given in Section 5 and Appendix C.

Once we have verified that the astrometry and zero-points of the images are robust, our v1 mosaics and catalogs will be made available via our PEARLS websites. ⁸⁰

We used both the SourceExtractor (Bertin & Arnouts1996) and ProFound (Robotham et al. 2018) packages to generate object catalogs from our processed NIRCam images. Both packages were designed to deblend close objects and find the object total fluxes. Details of these procedures are given by Windhorst et al. (2011, 2022), and we applied similar procedures to the PEARLS NIRCam images. The current paper focuses on single-filter PEARLS object catalogs, and so we use the single-filter mode of object detection with SourceExtractor. That is, we defer the SourceExtractor steps necessary to produce accurate object colors, such as dual image-mode extraction, the production of band-merged catalogs, and the application of PSF-matching and aperture corrections to future papers, which will study the colors of faint stars and galaxies (R. Ryan et al. 2022, in preparation) as well as high-redshift dropout candidates (e.g., Yan et al. 2022).

The SourceExtractor input parameters for the NIRCam images used a minimum detection threshold above sky of 1.5σ and nine connected pixels above this threshold for inclusion in the catalog. We found that using fewer connected pixels resulted in too many small spurious sources, particularly in the LW images, where the original pixel size is larger. While the more stringent nine-pixel requirement may result in missing a few real sources at the faint end, we found this to be a good compromise between reliability and completeness at all wavelengths. To detect sources, we used a 5 × 5 pixel convolution filter with Gaussian FWHM of 3.0 pixels (0090). This value is close to the median size of the faintest galaxies in Figures 6–8 (i.e., with object FWHM ∼ 01 or half-light radii r_e ≃ 005), which enables us to better detect very faint, low-SB, or clumpy galaxies. The SourceExtractor parameter DEBLEND_MINCONT was set to 0.06 to assure that real objects were not over-deblended. These parameters were chosen as a balance between extracting objects deep enough to achieve sample completeness to approximately the 5σ detection level (Section 4.2) but not so deep that a visible number of bogus objects were detected around remaining low-level image artifacts (Section 3.3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For all images and filters, the corresponding weight maps were used to account for image borders and properly characterize the photometric uncertainties and the effective areas for each drizzled mosaic. We compared our catalogs to the actual images and weight maps to look for any visible excess objects near residual image structures and snowballs and where needed applied additional masking to the images and their weight maps. Most masked regions are due to edge effects, low-exposure regions, or strips due to the dither pattern adapted for our shallow PEARLS exposures. Some bright stellar diffraction spikes (e.g., Figure 3) and residual wisp patterns also required masking. All of these masked areas were excluded when calculating galaxy number counts (Section 4.5).

The first line of Table 1 lists the J2000 tangent point to which the images in all filters of each target were drizzled, the observing date, the APT visit number, the area covered, the net exposure time per filter and total net hours, and the spacecraft efficiency of that visit. The second line for each filter lists the PSF FWHM, and the third line lists the 5σ point-source sensitivity in AB-mag predicted by the prelaunch ETC for the net integration time on the first line of each target. The fourth line in each filter indicates the achieved ∼5σ detection limits. The values were derived from Figures 6–8 as the median AB magnitude where SourceExtractor reports flux error bars of 0.20 mag. The fifth line indicates the AB level in Figures 6–8 where the galaxy counts are ∼80% complete compared to a power-law extrapolation, as derived from the figures in Section 4.5.

Separating stars from galaxies is important, especially for the PEARLS NEP fields, which are located at Galactic latitudes +316 (JWIDF) and +336 (TDF). Fields at these latitudes are expected to have more faint brown dwarfs than fields at higher Galactic latitudes (e.g., Ryan et al. 2011, 2017, 2022; Jansen & Windhorst 2018). The PEARLS star–galaxy classification procedure is based on the method described by Windhorst et al. (2011) for the ten-band WFC3 ERS images and by Windhorst et al. (2022) for the HST Archival Legacy project SKYSURF sample of ∼249,000 HST images. Figures 6–8 show the object detection, classification, and count diagnostics used for the JWIDF, as well as for the El Gordo noncluster module, respectively.

For each detected object, the left panels show the SourceExtractor magnitude error versus MAG_AUTO AB-mag. The horizontal dashed line indicates the adopted detection limit where the MAG_AUTO error is ≥0.20 mag, which approximately corresponds to a ∼5σ detection for point sources. Table 1 lists the average AB-mag values where the MAG_AUTO error reaches ≥0.20 mag, as derived from the left panels of Figures 6–8.

The middle panels shows the star–galaxy classification diagram using SourceExtractor MAG_AUTO AB-magnitudes versus image FWHM. The NIRCam diffraction limit is indicated in by the full-drawn left-most vertical line, and the FWHM of the PSF is listed in the legend of each middle panel and in Table 1. Objects with FWHM < FWHM(PSF) have been flagged and removed from this plot as spurious detections or border imperfections. In short, objects detected by SourceExtractor with sizes straddling the NIRCam diffraction limit to a certain magnitude limit are classified as stars (red dots), following Windhorst et al. (2011). Stars are generally located in a thin nearly vertical column bordered by the right-most vertical line. The remaining objects are classified as galaxies (blue dots). The green dotted–dashed lines indicate the 5σ sensitivity limits of each image, with the horizontal part showing the ∼5σ point-source limit and the slanted part showing the SB limit. For further details, see Windhorst et al.(2022).

The right panels of Figures 6–8 show the resulting star counts (red filled circles) and galaxy counts (blue filled circles). At most wavelengths ≳0.9 μm and at intermediate to high Galactic latitudes, the star counts generally have a very flat slope (γ ≃ 0.04 dex mag⁻¹), while the galaxy counts have relatively steep slopes γ ≃ 0.21–0.25 dex mag⁻¹, continuing the trend seen at 0.2–1.6 μm wavelengths by Windhorst et al. (2011). As a consequence, galaxies generally dominate the object counts for AB ≳ 18 mag and far outnumber Galactic stars at fainter magnitudes. Hence, reliable identification of faint objects as stars becomes difficult for AB ≳ 26–27 mag, and we have treated all objects fainter than this as galaxies. (The specific limiting values for each filter are shown in the right panels of Figures 6–8.)

Future work may be able to expand the identification of somewhat fainter stars through color–color diagrams and comparison with theoretical stellar loci. In a prior HST study (Windhorst et al. 2011), comparison of the ten-band WFC3 ERS star counts to Galactic-structure model predictions verified the star–galaxy classification procedure a posteriori. A similar check for the TDF data is described by R. Ryan et al. (2022, in preparation). Windhorst et al. (2011) also checked their star counts against spectroscopic ones from the HST ACS R ∼ 100 grism survey "PEARS" (Probing Evolution And Reionization Spectroscopically; Pirzkal et al. 2009). Such a posteriori verification of our star–galaxy classification procedure will be possible for the PEARLS TDF data when we receive the NIRISS grism spectra in all four NEP TDF epochs later in Cycle 1.

Star–galaxy classification in JWST NIRcam images is—ironically—hardest at the bright magnitude levels of AB ≃ 18–20 mag, as can be seen from Figures 6–8. These bright objects are generally unsaturated in the NIRCam images, but their significant diffraction spikes make FWHM an unreliable indicator. The spikes also make it difficult to derive accurate magnitudes and colors for identifying Galactic stars, as was done for the WFC3/IR images of Windhorst et al. (2011). We therefore used the Gaia DR3 (e.g., Gaia Collaboration et al. 2022) catalog to identify stars with AB ≲ 19 mag through their nonzero proper motions and the SDSS DR7 spectroscopic catalog to identify stars with available spectra at AB ≲ 17.5 mag.

As indicated above, the complex shape of the JWST PSF makes star–galaxy separation more difficult for objects with 18 ≲ AB ≲ 21 mag. Another complication is that an AGN can appear as a point source strong enough to hide the host galaxy and make the object appear stellar even in the JWST images. To verify whether the automated star–galaxy separation of Figures 6–8 was done correctly, we therefore proceeded as follows.

(a) Two independent visual observers inspected all objects (both those classified as stars and as galaxies) with 18 ≲ AB ≲ 21 mag to arrive at a consensus on which bright objects are stars and those that are compact galaxies with or without weak AGN. This was done in all 4–8 filters in the JWIDF and El Gordo noncluster fields. In the JWIDF, we found no objects classified among the 20 brightest stars (18 ≲ AB ≲ 24.4 mag) at the shorter wavelengths that were classified as galaxies at the long wavelengths. In El Gordo, we found eight possible galaxies among the 20 brightest objects classified as stellar. So in total the fraction of brighter stellar objects that are misclassified galaxies is about 20%.

(b) We required that stars needed to be classified as such in at least 2 out of 4–8 NIRCam filters using using the method in the middle panels of Figures 6–8. Windhorst et al. (2011) required a stellar object to be classified as such in at least 3 out 10 their HST ACS or WFC3 filters. This method found 69 objects classified as stellar in 2–4 JWIDF filters and 23 objects classified as stellar in two to eight filters in the El Gordo noncluster field. The surface density of stellar objects in Figures 6–8 is thus about 1000–2000 deg⁻²/0.5 mag in the El-Gordo noncluster module and JWIDF fields, respectively. The JWIDF is at lower Galactic latitude, and so has a higher surface density of stars. The results of (a) and (b) were largely consistent, and together reduced the number of bright stars that were misclassified as galaxies. This is reflected in the galaxy counts of Figures 9–10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Despite the above procedures, our objects classified as "stellar" could still be contaminated by faint, compact galaxies or weak AGN. To estimate the contamination level by (weak) AGN, we must consider their expected surface density at near-IR wavelengths. In the optical, the QSO surface density is known to be ≲15 deg⁻²/0.5 mag to B ≲ 21 mag (Boyle et al. 2000) and ≲125 deg⁻²/0.5 mag to B ≲ 23 mag (Koo & Kron 1982). To predict surface densities at the JWST NIRCam wavelengths, we used the UV–far-IR data sets from the GAMA (Bellstedt et al. 2020a) and DEVILS (Davies et al. 2021) surveys with multiwavelength SED fits by Bellstedt et al. (2020b), Thorne et al. (2021), and Thorne et al. (2022). These codes used ProSpect (Robotham et al. 2020), which fits stellar and AGN SEDs with an AGN fraction f_AGN at 1–2 μm wavelengths as a free parameter. We used the GAMA and DEVILS databases to estimate the surface densities of objects with 20 ≲ AB ≲ 25 mag and with SEDs in the VISTA ZYJH filters yielding f_AGN ≳ 0.5–0.9. Similar to the behavior at the optical wavelengths above, the surface density of such weak AGN converges to ≲ 100–50 deg⁻²/0.5 mag for AB(1–2 μm) ≲ 25 mag and f_AGN≳0.5–0.9, respectively. This amounts to ∼10% of our observed surface density of stellar objects of 1000–2000 deg⁻²/0.5 mag in the El-Gordo noncluster and JWIDF fields, comparable to our estimated ∼20% contamination rate of stellar samples by galaxies with weak AGN above. In conclusion, we expect that ∼10%–20% of our stellar samples may be contaminated by extragalactic objects. Future work that includes matched-aperture SED fits and PSF subtraction will be able to make a more accurate assessment of the fraction of compact galaxies or weak AGN remaining in our announced stellar samples.

Figures 6–8 show that the stellar locus moves steadily toward smaller FWHM values as wavelength decreases from 4.5 to 0.9 μm. This is a consequence of the JWST diffraction limit being much better than its requirement at 2.0 μm. That is, the stellar locus keeps moving to smaller object sizes from 2.0 to 1.5 μm and continues to do so down to 1.15 μm, although the stellar region does become somewhat wider in the F090W filter and in some fields in the F115W filter. Achieving a diffraction limit well below 2.0 μm wavelength is a major accomplishment for the JWST Project team, as this had to be planned between 2003 and 2005 without further driving up the Project cost. This was achieved by keeping the diffraction limit requirement at 2.0 μm but polishing the JWST mirrors well enough that the high-frequency wave front error would have a low enough rms to make a diffraction limit below 2.0 μm possible as long as the actuators below the mirrors could remove the main low- and mid-frequency errors well. The telescope in L2 is indeed able to do this (Rigby et al. 2022).

The consequences of this achievement by the JWST Project are far-reaching, as seen in the current paper: the median size of the faintest galaxies in the PEARLS images is about FWHM ≃ 01 as shown in all the SW panels of Figures 6–8. Hence, at JWST's diffraction limited SW resolution at ∼1.1–2.0 μm, essentially all faint galaxies in Figures 6–8 are resolved by NIRCam. In the LW 2.7–4.5 μm panels of Figures 6–8, a significant fraction of faint galaxies remain unresolved and therefore bunch up against the diffraction limits. Therefore, faint galaxies with flat SEDs are more easily detectable with the wider LW PSF of ∼017 FWHM compared to the SW PSF, which has 006–008 FWHM. This is especially visible in the filters F277W and longwards.

Table 1 quantifies the detection limits, using the 80% galaxy-count completeness limits as a fiducial value. These limits were determined from power-law fits to the counts and their extrapolations (Section 4.5). In all NIRCam SW filters, the 80% galaxy count completeness limits typically appear to be about −0.3 to −0.9 mag brighter than the 5σ point-source detection limits predicted by the prelaunch ETC. This is indicated by the ${\rm{\Delta }}{\mathrm{AB}}_{\mathrm{lim}}(80 \% -\mathrm{ETC})$ values on the sixth line for each target in Table 1. A small part of this difference is due to the fact that some SW filters can be ∼10% less sensitive than the prelaunch predictions (Rigby et al. 2022), but for the most part this "apparent loss" in SW point-source sensitivity occurs because the large majority of faint galaxies observed in NIRCam SW are no longer point sources. The simple reverse is true for the faintest galaxies in most NIRCam LW filters: our PEARLS LW galaxy counts appear to be between +0.0 and +0.7 mag more sensitive than the prelaunch ETC prediction for point sources. Part of this difference occurs because most LW filters are 20%–30% more sensitive than the prelaunch predictions of Rigby et al. (2022). But in addition, the detection of the faintest galaxies in the LW images is surely aided by the fact that a significant fraction of the faintest galaxy sizes no longer exceeds the size of the PSF FWHM in the LW filters. The tendency of faint galaxies to bunch up against the HST diffraction limit at brighter flux levels was first suggested based on the Hubble Deep Field images by Odewahn et al. (1996) and Windhorst et al. (1998), and later by Welch et al. (2022c), Windhorst et al. (2021), and references therein based on more recent HST images.

In conclusion, at the FWHM of the NIRCam PSFs delivered by the JWST Project, faint galaxies with flat SEDs are noticeably easier to detect in the LW filters compared to the SW filters. Late-type stars are bluer than galaxies and are point sources in all NIRCam filters, so this PSF-advantage at longer wavelengths for galaxies does not apply to stars. Current and future JWST surveys that search for high-redshift dropouts and other red objects will need to keep this rather strongly wavelength-dependent sensitivity to the typical faint-galaxy sizes into account.

1. The PEARLS 0.9–4.5 μm Galaxy Counts: our first PEARLS galaxy counts are based on the data in two of our four fields observed so far with the best available data in Table 1. We used all 10 NIRCam detectors in the JWIDF plus the four SW and one LW detector in the El Gordo noncluster module. A comparison of the two sets of galaxy counts in the four filters where they overlap then enables us to quantify whether the El Gordo noncluster module was biased in a measurable way by the presence of the z = 0.870 cluster in the adjacent NIRCam module. We defer reporting galaxy counts in the VV 191 field to a later paper, owing to the presence of the two large, targeted galaxies and bright objects in the same nearby galaxy group that cover some of the other detectors. We also do not present object counts in the five medium-band filters of the TNJ-1338−1942 protocluster at z = 4.1 because these medium-band images are shallower than the broadband images in Table 1 and because there are no reference filters available from the ground for comparison to brighter object counts. Nonetheless, we did carry out their star–galaxy classification and object counts, as was done for the JWIDF and El Gordo in Figures 6–8, to check on the reliability of our procedures and to report their 5σ point-source sensitivities in Table 1 as compared to the ETC predictions.

We also inspected all 18 ≲ AB ≲ 21 mag objects in the noncluster field of El Gordo and found 10 galaxies close to the cluster outskirts that were in the noncluster module and have colors very similar to the El Gordo cluster galaxies. These are likely part of the outskirts of El Gordo. We removed these 10 objects from the galaxy counts. The galaxy counts in the JWIDF and El Gordo noncluster module at the JWST bright end (18 ≲ AB ≲ 20 mag) are consistent to within their (large) error bars and are close to the average counts from the previous ground-based, HST, and Spitzer surveys, which have much smaller error bars over this magnitude range. In any case, the error bars on the bright end of the JWST galaxy counts (18 ≲ AB ≲ 20 mag) are large enough that they do not weigh significantly into the spline fits of the galaxy counts (middle and right panels of Figures 9–10).

Our PEARLS galaxy counts in the JWIDF and El Gordo noncluster module are shown in Figures 9–10. These are based on our object catalogs of Section 4.1 with objects identified as stars removed (Sections 4.2–4.3). Error bars reflect the statistical uncertainties in the remaining galaxy counts. At the bright end (AB ≃ 18–21 mag), larger discrepancies are seen in the counts between the two fields due to the uncertainties in the star–galaxy classification procedure (Section 4.2–4.3) and due to cosmic variance. Because we used two NIRCam fields far apart in the sky, their cosmic variance is expected to be ≲9% (e.g., Driver & Robotham 2010, and Section 2.1 above). This estimate uses the JWIDF+El Gordo survey area of Table 1 and the redshift distribution expected for NIRCam objects with 18 ≲ AB ≲ 28 mag, assuming most objects are in the redshift range 0.3 ≲ z ≲ 8. Further details are given in Appendix B.2. For the brighter fluxes of 18 ≲ AB ≲ 21 mag, CV can be as high as 20% for redshifts 0.1 ≲ z ≲ 1, explaining some of the remaining statistical variations seen at the bright end of the counts in Figures 9–10.

The bright end of the 0.9–4.5 μm JWIDF and El Gordo counts in Figures 9–10 are consistent with each other to within their error bars. Hence, we see no evidence that the galaxy counts in the El Gordo noncluster module are significantly higher than those in the JWIDF, which is a random survey field. We will thus proceed with the JWIDF and El Gordo noncluster module counts as representative for the 0.9–4.5 μm galaxy counts, and will use a CV error of ∼9% in our discrete IGL error budget over the entire magnitude range of 18 ≲ AB ≲ 28.5 mag in Figures 9–10.

A few more objects may reside in the large-scale structure associated with the El Gordo cluster at z = 0.870. These may be removed from the galaxy counts in the El Gordo noncluster module when more spectra of the cluster and its surroundings become available. Future work will improve the accuracy of galaxy counts when done over more JWST fields. Given the high surface density of background galaxies detected by NIRCam and the negative magnification bias predicted by the shallow NIR count slopes in Figures 9–10, more accurate background-galaxy counts combined with a weak shear analysis of the same images may improve mass profile measurements of the galaxy clusters (e.g., Umetsu et al. 2011).

2. Comparison to Previous Galaxy Counts at 0.9–4.5 μm: the PEARLS 0.9–4.5 μm galaxy counts are compared to those at brighter levels from the combined GAMA (Driver & Robotham 2010; Driver et al. 2011, 2022) and DEVILS (Davies et al. 2021) surveys at similar wavelengths in Figures 9–10 and at the shorter wavelength (0.9–1.6 μm) also with the deepest available HST counts. A possible check of catalog completeness and reliability is to compare our JWST/NIRCam object counts to the deepest available HST counts in the same or similar filters and see whether the agreement is good to within the known ZP, rms counting, and cosmic-variance errors. This will also help verify the flux level at which catalog incompleteness sets in.

The GAMA and COSMOS/DEVILS survey data were compiled over a wide range of wavelengths by (Driver et al. 2016a, and references therein), building on the GAMA panchromatic data release of Driver et al. (2016b) and the COSMOS data as reanalyzed by Andrews et al. (2017). This compendium was later extended by Davies et al. (2018) and Davies et al. (2021) as part of the DEVILS survey, while Bellstedt et al. (2020a) extended the GAMA data to also include the ESO VST KiDS data. These updates were reported by Koushan et al. (2021), whose work forms the basis of the galaxy-counting data used here. These compendia also include the ESO K bands counts of Fontana et al. (2014), deep Spitzer 3.6 and 4.5 μm counts (e.g., Ashby et al. 2009, 2015; Mauduit et al. 2012), and the WISE counts of Jarrett et al. (2017). The typical combined ZP uncertainties in the combined GAMA+DEVILS surveys in the z band to K band are 2%–3% after bringing the flux scale in every filter onto the flux scale of the VISTA survey filters, which incorporated the GAMA and DEVILS surveys (Table 4 of Koushan et al. 2021).

The deepest panchromatic HST ACS+WFC3 galaxy counts come from the combined HUDF images, whose database has grown considerably over time since the launch of WFC3 in May 2009 (see e.g., Windhorst et al. 2011; Koekemoer et al. 2013; Rafelski et al. 2015, and references therein). Following the rich HUDF database summarized in these papers, the panchromatic galaxy counts were once more repeated on the deepest available HUDF images in 2015 with the same procedures as in Sections 4.1–4.2 and were included by Driver et al. (2016a) and Koushan et al. (2021). This deepest 2015 realization of the HUDF counts is listed in the legend of Figures 9–10 as "(Windhorst et al. 2011; ⁺)" and follows the same methods. The ZP errors in the HST ACS, WFC/UVIS, and WFC3/IR images over the decades resulted in flux scales accurate to 1%–3% as summarized in Section 4.1.5 and Table 6 of Windhorst et al. (2022), i.e., comparable to or slightly better than the 2%–3% ZP accuracy of the VISTA filters of Koushan et al. (2021).

To within these ZP errors, our 0.9–1.5 μm PEARLS galaxy counts that come from shallow NIRCam exposures are consistent with the deeper HUDF counts in the HST ACS and WFC3/IR filters F850LP, F105W/F125W, and F160W (green open circles in Figures 9–10). It is important to realize that our PEARLS galaxy counts come from 1890–3157 s NIRCam exposures in the F090W, F115W, and F150W filters and reach ∼28.5 mag (Table 1), while the HUDF galaxy counts reach ∼29.5–28.5 mag in ∼156–87 HST orbits (≃117–65 hr of net exposure time; Beckwith et al. 2006; Koekemoer et al. 2013) in the F850LP, F105W/F125W, and F160W filters, respectively. Compared to the total HUDF ACS exposure time of 421.6 ks in F850LP, JWST NIRCam thus reaches approximately ∼5× deeper per unit time in its F090W filter, while compared to the total HUDF WFC3/IR exposure time of 236.1 ks in F160W, NIRCam reaches about ∼9–11× deeper per unit time in its F150W filter, respectively. In the light of the discussion in Section 4.4, this is an impressive performance improvement, especially because NIRCam was optimized for performance longwards of 2.0 μm.

Because our PEARLS 0.9–4.5 μm galaxy counts are done in NIRCam filters whose effective wavelengths and bandpasses can be somewhat different from the VISTA, GAMA, and DEVILS surveys, Appendix B.2 addresses whether additional corrections to the NIRCam flux scale are needed to compare the galaxy counts in similar filters. For this, we used the fiducial flux scale of the VISTA filters into which the GAMA and DEVILS surveys were anchored (Koushan et al. 2021). Appendix B.2 shows that the corrections needed to transform the NIRCam AB-mag scale onto the fiducial VISTA/IRAC filters are ≲3%–4% with combined uncertainties of ≲3%–6%, and therefore no corrections to the NIRCam AB-mag scale needed to be applied when plotting the results in Figures 9–10. However, we folded this ≲3%–6% uncertainty into our error budget of the 0.9–4.5 μm IGL of Section 4.6. We now have all the ingredients in place to compare our 0.9–4.5 μm NIRCam galaxy counts to previous work at brighter levels and can do so without further wavelength dependent ZP corrections.

3. The Faint-end Slope of the Combined 0.9–4.5 μm Galaxy Counts: Figures 9–10 compare our PEARLS 0.9–4.5 μm galaxy counts to previous work summarized above and extend it to AB ≲ 29 mag over this wavelength range. Over the AB ≃ 16–29 mag range for which they are available at 0.9–4.5 μm, the Yung et al. (2022) models are consistent with the combined GAMA/DEVILS ground-based, Spitzer/WISE, and our PEARLS NIRCam galaxy counts. The faint-end slopes of our observed galaxy counts have an average value γ ≃ 0.23 ± 0.04 dex mag⁻¹ for 22 ≲ AB ≲ 29 mag, where the counts are a nearly straight power law (Figures 9–10, Table 3). The middle panels of Figures 9–10 showing the 0.9–4.5 μm IGL energy counts were derived from the left panels by dividing by the 0.4 slope. The spline extrapolations (gray lines and error fans) in the middle panels show that energy counts are clearly converging at all these wavelengths. The resulting IGL integral is shown in the right panels of Figures 9–10. Under the assumption that the faint galaxy counts at AB ≳ 29 mag continue as a power law with the same slope as observed for AB ≃ 22–29 mag, these spline extrapolations will form the basis of our IGL values used in Sections 4.6 and 5. Further details on the faint-end slope of the galaxy counts and its wavelength dependence are given by S. Tompkins et al. (2022, in preparation).

A magnitude slope γ ≃ 0.23 mag dex⁻¹ corresponds to a faint end slope α ≃ − 1.58 ± 0.1 in flux units, where α =− 1 − 2.5γ. Ground-based spectroscopic surveys with VLT/MUSE and the spectrophotometric survey 3D-HST with Hubble suggest that faint galaxies with AB ≃ 23–29 mag have a median redshift in the range z_med ≃ 1–2 (e.g., Skelton et al. 2014; Inami et al. 2017) with the caveat that the completeness of these surveys becomes more difficult to quantify at fainter magnitudes. Around this median redshift, the faint end of the galaxy counts samples the power law part of the Schechter luminosity function (LF), which also has a faint-end flux slope α ≃ − 1.4 to −1.5 at z ≃ 1.5 (e.g., Hathi et al. 2010; Finkelstein 2016). In conclusion, over the AB ≃ 22–29 magnitude range sampled by our 0.9–4.5 μm NIRCam images, the PEARLS galaxy counts have a slope consistent with the faint-end slope of the Schechter LF at the median redshift sampled by these objects. Fainter JWST galaxy counts would then be expected to continue with the same slope, if the LF over this redshift range were to continue with the same slope toward fainter luminosities. Upcoming ultradeep JWST NIRCam GTO surveys (M. Rieke PI) are designed to cast light on this issue.

The parameters that best characterize the 0.9–4.5 μm galaxy counts and integrated galaxy light are summarized Figure 11. The errors on these parameters are summarized in Table 3 for each filter as the quadratic sum of the NIRCam ZP uncertainties from Appendix B.1, and the uncertainty in bringing the flux scale of the NIRCam filters onto the effective wavelengths of the VISTA/IRAC filters used as fiducial for the galaxy counts (Section 4.5). Some interesting trends can be seen in the galaxy counts over the wavelength range 0.9–4.5 μm from Figures 9–10. (These trends are best seen if all PNG files of Figures 9–10 are shown on a computer screen at high magnification in rapid succession—all panels are plotted at exactly the same scale for this purpose.) The smooth behavior of the data in Figure 11 suggests that these trends in the IGL are real and meaningful.

• 1.  
  The left panels of Figures 9–10 all show a clear change in slope from a steep nonconverging slope (≥0.4 dex mag⁻¹) to shallow and converging slope (<0.4 dex mag⁻¹). This change sets in around AB ∼ 19.3–20.3 mag at 0.9–4.5 μm (top panel of Figure 11). This peak AB-mag is the flux level where most of the IGL is generated, and is a clear function of wavelength.
• 2.  
  In more detail, the normalized differential counts reaches the highest SB-level around 2–3 μm wavelength, and declines to both longer and shorter wavelengths (first and second panel of Figure 11, where the second panel shows the peak SB-value of the IGL at the AB-magnitude peak of the top panel).
• 3.  
  The magnitude range over which most of the IGL is generated (here called the "IGL FWHM" and measured as the interquartile 25%–75% range of the middle panels in Figures 9–10) decreases from ∼4.5 mag at wavelengths ≲1.25 μm to ∼3 mag at 4.5 μm (third panel of Figure 11). This reflects the luminosity function and redshift distribution of the older, earlier-type galaxies that dominate the 3.5–4.5 μm galaxy counts. These filters sample the rest-frame ∼1.5 μm peak in the stellar emission of early-type galaxies as caused by their stellar mass distribution at z ≃ 1–2. The JWST images at these wavelengths are visibly dominated by earlier-type galaxies. At shorter wavelengths, we sample a larger fraction of galaxies of later-types, which have a wider range of ages and extinction, and extend to lower luminosities and redshifts, causing their larger contribution to the IGL at bluer wavelengths (e.g., Driver et al. 1995, 2016a; Andrews et al. 2018; Koushan et al. 2021). This can be seen in the larger IGL-width at bluer wavelengths in the third panel of Figure 11.
• 4.  
  The discrete IGL with the highest energy (in units of nW m⁻² sr⁻¹) comes from wavelengths between 1 and 2 μm, as shown in the bottom panel of Figure 11. These are derived from the converging integrals from the right panels of Figures 9–10, and are the values we plot in the last Figure of Section 5, which provides further discussion of the total IGL and diffuse light.

Zoom In Zoom Out Reset image size

The new JWST results are most noticeable at the faint end of the AB-magnitude scale plotted in Figures 9–10. The bright end of the galaxy counts can be—and has been—done from the ground (i.e., the z, Y-, J-, H-, and K-band filters) or from space with WISE and Spitzer at L (3.5 μm) and M band (4.5 μm). Nevertheless, JWST NIRCam and also HST WFC3/IR below 1.6 μm have unique filters that are valuable for object counts. These include WFC3 F140W, and NIRCam F277W and F410M, as well as the other medium-band filters in Table 1. These filters have no ground-based counterparts because telluric water vapor blocks these wavelengths. Therefore bright-end galaxy counts to make a full energy integral as in Figures 9–10 are absent for the F277W and F410M as well as the medium-band filters.

Because of this, Carleton et al. (2022) had to interpolate the IGL integral in the WFC3/IR F140W filter from the adjacent WFC3/IR F125W and F160W filters, which have extensive ground-based coverage of the bright end counts. Fortunately, this is straightforward because the IGL SB is flat between 1.25 and 1.65 μm wavelengths (bottom panel of Figure 11). We therefore give low-order spline functions fit to the IGL parameters versus wavelength in Table 3 and plot these in Figure 11. We use these splines to interpolate the IGL SB-values for the JWST NIRCam F277W and F410M filters, which are tabulated as the PEARLS-IGL values in Tables 4–6. This includes the values beyond the PEARLS NIRCam detection limits of AB>28.5 mag, which were derived from the integrals in the right-hand panels in Figures 9–10. This allows us to estimate the IGL as a function of wavelength for the JWST filters at 0.9–4.5 μm wavelength, including the ≲2.5% of the IGL that is not included in our faint object counts to AB ≲ 28.5 mag (see Section 5.2 here and Section 3.4.3 of Carleton et al. 2022 for its procedure). Wider-area JWST surveys such as e.g., COSMOS-Webb (PI J. Kartaltepe) will become available during JWST's lifetime to improve the bright-end of the galaxy counts, as they have for HST ACS and WFC3 during the last two decades (e.g., Windhorst et al.2022).

JWST's ability to work continuously in a dark-sky environment makes it especially suitable for measurements of sky-SB. This is in contrast with HST, which at best gets complete dark time for at most ∼30 minutes of its 96 minutes orbit (e.g., Caddy & Spitler 2021; Caddy et al. 2022; Windhorst et al. 2022). Figures 12–13 summarize the astrophysical foreground and background energy relevant to PEARLS compared to recent data (as summarized by, e.g., Driver et al. 2016a; Koushan et al. 2021; Carleton et al. 2022) and IGL models (e.g., Andrews et al. 2018). Figure 13 shows the PEARLS discrete IGL measurements of Section 4.6 compared to those of D16 and Koushan et al. (2021). The IGL is the sum of the integrated (observed) galaxy counts (iEBL) and extrapolated galaxy counts (eEBL) derived in Section 4.5. The black line shown in Figure 13 is a modification of the Andrews et al. (2018) IGL model for accumulated star formation in spheroids (red dashed), disks (green dashed), and unobscured AGN (purple dashed lines). Here we have adjusted these contributing elements from the published Andrews et al. (2018) model, as described in the caption, to better fit the IGL data including the PEARLS points.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

JWST was meticulously designed and built to have the darkest possible sky as seen from L2. Here we explore its capability to constrain potential levels of diffuse light. Windhorst et al. (2022) stated that over 95% of the 0.6–1.25 μm photons in the HST archive come from the Zodiacal light in the interplanetary dust (IPD) cloud, i.e., from distances <5 au. This can also be seen by comparing the typical Zodiacal light levels (green line) to the IGL counts in Figure 13. The Zodiacal/IGL ratio decreases significantly toward longer wavelengths in the 1.5–3.5 μm wavelength range. This is because the Sun is a zero redshift 5770 K G-star, and the IGL is the summation over multiple stellar populations, including hotter and cool stars, spanning a wide redshift range (e.g., Madau & Dickinson 2014). Longwards of 3.5 μm, thermal radiation from the Zodiacal belt and also from the JWST telescope and instruments make increasing contributions to the SB levels. The wavelength dependence of the diffuse light is precisely what JWST can explore from its first images, bearing in mind the significant JWST and NIRCam calibration uncertainties that we expect (Section 3.3 and Appendices B.1–B.3).

Obtaining diffuse light (DL) estimates requires accurate modeling of the Zodiacal light (ZL) and diffuse Galactic light (DGL). ZL can be ∼10–70× higher than the discrete iEBL+eEBL, dependent on the direction and time of observation (Figures 12–13). Constraints from previous work shown in Figure 13 suggest that there may exist some level of diffuse light at 0.6–1.6 μm, at a level of ∼8–30 nW m⁻² sr⁻¹. (Note that all diffuse light estimates plotted in color in Figure 13 have the full IGL already subtracted, and so truly represent the diffuse light levels or limits reported by various groups.) At this stage, it is not clear whether this diffuse light is due to residual instrumental systematics that have not been accounted for, a dim Zodiacal component (perhaps spherical or spheroidal) seen from 1 au that is not accounted for in the Zodiacal models, a truly diffuse EBL component, or some combination of these possibilities. For details on this topic, please see the discussions by, e.g., Conselice et al. (2016), Matsuura et al. (2017), Sano et al. (2020), Carleton et al. (2022), Korngut et al. (2022), Kramer et al. (2022), Lauer et al. (2022), O'Brien et al. (2022), and Windhorst et al. (2022), and references therein.

Following Equation (2) of Windhorst et al. (2022), the sky-SB level between the detected objects is a sum of ZL, DGL, and residual instrumental systematics including thermal and straylight contributions. For JWST, that equation is:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{SB}(\lambda ,{l}^{\mathrm{Ecl}},{b}^{\mathrm{Ecl}},{l}^{\mathrm{II}},{b}^{\mathrm{II}},t,\mathrm{SA},T)\\ \ =\,\mathrm{Th}(\lambda ,T)+\mathrm{SL}(\lambda ,{l}^{\mathrm{Ecl}},{b}^{\mathrm{Ecl}},t)+\mathrm{ZL}(\lambda ,{l}^{\mathrm{Ecl}},{b}^{\mathrm{Ecl}},t,\mathrm{SA})\\ \ +\,\mathrm{DGL}(\lambda ,{l}^{\mathrm{II}},{b}^{\mathrm{II}})+\mathrm{dEBL}(\lambda ).\end{array}\end{eqnarray} \tag{ 3 }$$

The left term in Equation (3) is the total sky-SB that JWST observes as a function of wavelength λ, Ecliptic coordinates (l^Ecl, b^Ecl), Galactic coordinates (l^II, b^II), time of the year (t or Modified Julian Date), solar elongation angle (SA), and telescope and instrument temperatures, symbolized by T. The terms on the right side include: thermal (Th) signal from blackbody photons in the instruments and telescope that depends on wavelength and temperature; straylight (SL) that depends on wavelength, pointing direction, and observing date (Section 3.3); ZL as seen from L2 that depends on wavelength, Ecliptic coordinates, and observing date via its effect on the SA and on the path through the Zodiacal dust cloud (especially JWST's position above or below the Ecliptic plane; see Appendix C); DGL that depends on wavelength and Galactic coordinates; and any diffuse EBL (dEBL) that is not already included in the discrete object catalogs to AB ≲ 28.5 mag in Section 4, and therefore not yet masked out from our NIRCam images. This last term includes the part of the discrete IGL extrapolated for AB ≳ 28.5 mag (i.e., the eEBL), which is generally small and subtracted below.

At the start of JWST Cycle 1, we only have a limited number of JWST images and so can only explore limits for the sky-SB values observed from L2 thus far by PEARLS. When more JWST images become available, its full database can be studied following the SKYSURF methods of Carleton et al. (2022), who presented the HST sky-SB between discrete objects in 34,000 WFC3/IR images, and set constraints on diffuse light from the subset of HST 1.25–1.6 μm images with the darkest sky-SB values. At this stage, we will use our darkest available JWST images to explore what constraints can be made currently and how these may be improved in the future during JWST's lifetime.

We measured sky-SB in the NIRCam images following the SKYSURF procedures of Windhorst et al. (2022) and Carleton et al. (2022). In short, SKYSURF removes the light from all detected objects from the WFC3/IR images even in short HST exposures (${t}_{\exp }$ ≃ 500 s) to a total-object flux limit of AB ≲ 26.5 mag at 1.25–1.6 μm wavelengths. For JWST NIRCam, we used the same codes to remove the light from all detected objects (Section 3.3 and Appendix B.3) to AB ≲ 27.5–28.5 mag at 0.9–4.5 μm given the detection limits in Table 1, at which flux levels ≳97.5% of the discrete IGL is already detected in the JWST images, as the right panels of Figures 9–10 show. Details on the uncertainty in the estimated sky-SB are given in Appendix B.3.

Tables 4–6 summarize the JWST NIRCam instrumental and astronomical background levels as predicted from, or actually observed from L2 for each of the four PEARLS targets observed as of 2022 July 31. Background components are assumed to be uniform across the field of view but depend on wavelength. Tables 4–6 list the ETC predictions for the L2 Zodi, thermal and straylight, as well as the straylight level of Rigby et al. (2022; using their Figure 5). Details on these predicted sky-SB component values and their uncertainties are given in Appendix C, and we summarize aspects relevant to current discussion here. The sky-SB predictions for all these components and their uncertainties (where relevant) are given in Tables 4–6, and include:

• 1.  
  The ETC-predicted JWST thermal radiation, which is more than 100× lower than the predicted total sky-SB even at 4.5 μm (see Figure 12) and is the dimmest component in Equation (3) for λ ≲ 4μm.
• 2.  
  The L2 model prediction for the Zodiacal sky-SB for each target. This was based on the position and orientation of JWST at the actual time of the observation. These are based on the Kelsall et al. (1998) model, but for the Zodiacal cloud geometry for L2 as seen by JWST at the time of the observation. These predictions are uncertain by at least ∼9%–4% of the dimmest Zodiacal sky-SB observed over the range 1.25–4.5 μm, respectively.
• 3.  
  The IPAC IRSA prediction for the DGL value at the Galactic coordinates of the target. The DGL is generally a factor of 20–100× lower than the total predicted JWST sky-SB, and uncertain by up ∼0.3 dex.
• 4.  
  The actual straylight, which Rigby et al. (2022) noted is likely lower than the preflight predictions. Indeed, using the full Rigby et al. (2022) SL values would make the total sky-SB predictions from Equation (3) exceed the observed values in two of our PEARLS fields in Tables 4–6. As discussed above, the 0.9–3.5 μm Zodiacal component is caused by sunlight scattered off the Zodiacal dust cloud components, and may have been underestimated in some of the models. The 3.5–4.5μm sky-SB is dominated by the thermal contribution from the Zodiacal dust cloud components, while the telescope+instrument thermal components are still negligible in these filters. The thermal Zodiacal Light at λ ≳ 3.5 μm was the key component to be modeled by Kelsall et al. (1998) for their COBE/DIRBE analysis. The minimum sky-SB is predicted to occur around 3.5 μm in wavelength (Figures 12–13), so we will assume that: (a) the thermal Zodiacal components at λ ≳ 3.5 μm are more accurately predicted than the scattered sunlight components at λ ≲ 3.5 μm; and (b) the predicted thermal Zodiacal components should match the values observed at 4.5 μm without exceeding those observed in the minimum at 3.5 μm. Any truly diffuse astrophysical source is expected to be much dimmer than this, so we do not expect it to significantly affect our fitting. In order to not exceed the observed PEARLS sky-SB values, we then find that the implied SL values are generally f ∼ 50%–100% of the Rigby et al. (2022) SL values, with an uncertainty in f of at least ∼20% in Equation (4) below (see also Appendix C). With these adopted SL values in Tables 4–6, our total predicted JWST sky-SB matches the observations in all 13 PEARLS filters in Figure 12 to within the uncertainties summarized above, using Equations (4)–(5) below.
• 5.  
  We use the IGL integral for the seven fiducial filter wavelengths in Section 4.5 from Figures 9–10. For the NIRCam wavelengths for which a full IGL integral is not yet available (F277W, F410M, and the medium-band filters), we used the spline predictions at those wavelengths from Figure 11. The IGL is assumed to be constant across the sky, and therefore to be the same for each PEARLS target. The extrapolated discrete eEBL (eEBL) of objects currently undetected in the JWST images for AB ≳ 28.5 mag is derived from Figures 9–10 (Section 4.6). The extrapolation for AB ≳ 28.5 mag typically amounts to ∼2.5% of the total IGL. This is the only part of the IGL that still need to be subtracted from the sky-SB, as our method already automatically removes ∼97.5% of the light from all brighter objects detected by NIRCam to AB ≲ 28.5 mag (see Carleton et al. 2022, for details). (Section 4.7 and Equation (3) of Windhorst et al. (2022) also corrected this fraction for SB incompleteness, which can be ∼40% at AB ∼ 28 mag, but the IGL correction for known objects at AB ≳ 28.5 mag remains very small.)

Our best prediction of the observed JWST Zodiacal sky-SB is then:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{J}}{\rm{W}}{\rm{S}}{\rm{T}}({\rm{P}}{\rm{r}}{\rm{e}}{\rm{d}}) & = & \mathrm{ETC}(\mathrm{Thermal})+f\times SL(\mathrm{Rigby}\,2022)\\ & & +\,{\rm{Z}}{\rm{o}}{\rm{d}}{\rm{i}}({\rm{L}}2)+{\rm{D}}{\rm{G}}{\rm{L}}+{\rm{e}}{\rm{E}}{\rm{B}}{\rm{L}}.\end{array}\end{eqnarray} \tag{ 4 }$$

Tables 4–6 give this sum as predicted from the above models in all filters at the actual time of PEARLS observations. Finally, these tables list our upper limits to any DL as the difference between the observed JWST sky-SB "JWST(Obs)" and the total prediction of Equation (4) for L2:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Diffuse}\ \mathrm{light}\ \mathrm{limit}\lesssim \mathrm{JWST}(\mathrm{Obs})\\ \quad -\,\mathrm{JWST}(\mathrm{Pred})\pm \mathrm{error}\,\mathrm{budget}.\end{array}\end{eqnarray} \tag{ 5 }$$

The results from Equations (4)–(5) are listed on the bottom lines in each tier of Tables 4–6. This includes the full error budget of ∼6%–8% for JWST(Obs) in the LW–SW modules, respectively, from Appendix B.3, and the combined uncertainty of ∼10% in JWST(Pred) from Appendix C. The total uncertainty in the difference of Equation (5) is thus ∼12%–13% of the total sky-SB in the LW–SW modules, respectively, assuming that the observed and predicted sky-SB values are independent.

With our assumption that the total sky-SB model should fully predict the observed sky-SB values in the four PEARLS filters at 3.5–4.5 μm, the model predictions generally also match the observed sky-SB in the seven PEARLS filters at 0.9–3.5 μm, including the medium-band filters in TNJ1338, within the combined uncertainties. Therefore, to within the error budget of the current assessment, we have no firm detection of remaining DL by JWST NIRCam.

Accordingly, all our diffuse light constraints are plotted as upper limits using the combined uncertainties of Appendices B–C in Figure 13. Brown downward open triangles indicate upper limits from our deepest filter exposures in the JWIDF and the El Gordo noncluster field. Brown downward asterisk and tripod shape indicate the upper limits from the shallower TNJ and VV191 exposures. We excluded in this process the detectors that contained the overlapping nearby galaxy pair of Figure 5 and other large objects. Our PEARLS constraints in Figure 13 indicate upper limits to diffuse light captured by the gray-hashed area, and generally amount to 12%–13% of the total sky-SB observed by NIRCam. Within the current uncertainties in the JWST NIRCam calibration and in the total JWST sky-SB model, we cannot make firmer statements about the diffuse light as seen by JWST. In particular, if the JWST SL were even lower than we adopted here, firmer constraints on DL may be made. For this purpose, future work will require a more accurate assessment of the NIRCam calibration uncertainties, and more accurate models for the JWST straylight, the Zodiacal Light as seen from L2, and for the DGL.

The two JWIDF diffuse light points at 1.5 and 2.0 μm are marginally above the total model predictions in Figure 12(d). Our F150W and F200W NIRCam diffuse light limits in Figure 13 are in line with the 1.1–1.6 μm CIBER detections of Matsuura et al. (2017) and Sano et al. (2020; purple triangles in Figure 13), and with the SKYSURF upper limits in the HST/WFC3 F140W and F160M filters of Carleton et al. (2022). These papers, as well as Tsumura (2018) and Korngut et al. (2022), suggested that some very dim spherical—or nearly spherical—Zodiacal component could be missing from the Kelsall et al. (1998) model. Kelsall et al. (1998) also noted that a dim spherical Zodiacal component could have been missed in their model of the COBE/DIRBE data.

At the longer NIRCam wavelengths of 2.7–4.5 μm, the PEARLS DL limits in Figure 13 are lower in value, reaching as low as 8–12 nW m⁻² sr⁻¹, and consistently so between our four PEARLS fields within the current error budget. This is because the total sky-SB is significantly darker (Figure 12), and the total error budget of the LW modules (Appendix B) and the uncertainties in the sky-SB models are correspondingly smaller at 2.7–4.5 μm (Appendix C).

JWST photometric calibration ("zero-point" or ZP) has to be established in flight and will evolve during the mission. Because standard-star observations were not available in time, the earliest JWST observations had to be calibrated with preflight ZPs. The MAST pipeline for NIRCam began using in-flight ZPs on 2022 July 27. At NIRCam wavelengths ≤2.0 μm, the on-orbit throughput was near prelaunch expectations (Rigby et al. 2022), but it was up to 10% smaller for some filters. At wavelengths >2.0 μm, the on-orbit throughput was about 15%–30% higher than the prelaunch expectations (Rigby et al. 2022; their Figure 8). Rigby et al. (2022) also wrote that the throughput stability is no worse than 4% and is likely much better.

With the new ZPs derived from on-orbit standard star observations in all NIRCam detectors and its main filters (Boyer et al. 2022), we reprocessed all our PEARLS images with jwst_0995.pmap_filters (Section 3). For the record, the results from our original processing with jwst_0916.pmap_filters through jwst_0952.pmap_filters are still available in the first submission of this paper on https://arxiv.org/abs/2209.04119, which also had a table of additional ZP corrections that this earlier processing required. The application of jwst_0995.pmap_filters has made these additional ZP corrections obsolete. It is nevertheless useful to give a brief comparison of the main differences between jwst_0995.pmap_filters (v1) and our earlier PEARLS images (v0.5):

• 1.  
  The latest v1 calibration more accurately corrects for ZP variations between each of the 10 NIRCam detectors (typically by ≲10%–20% per detector). The improvement propagates into our images and science results although in a rather subtle way, because our previous processing already had averaged over 2–8 detectors and 2–4 fields. In general the change reduces the dispersions in the sky-SB values but does not reduce their medians very much, as described below.
• 2.  
  The JWIDF has IRAC observations that allow direct comparison in two NIRCam filters. These IRAC observations were accumulated every two weeks over 15 yr and so are very deep (e.g., Yan et al. 2018). For unsaturated, isolated, and matched objects in the magnitude range 18 ≲ AB ≲ 20 mag in both sets of images, Source Extractor MAG_AUTO gives average total flux differences of F356W–IRAC1 ≃ –0.026 mag and F444W–IRAC2 ≃ –0.029 mag. Small differences between the JWST and Spitzer fluxes can be caused by a combination of filter differences, aperture corrections, and source confusion due to the vastly different PSFs of the two telescopes. The NIRCam PSF has a 2.36² ≃ 5.56× larger area in F444W compared to F090W (Table 1 and Figures 6–8). F444W also has a ∼140× smaller PSF area than the Spitzer IRAC2 filter. Despite these PSF differences, the 2022 October NIRCam F356W and F444W zero-points are consistent (within 2.6%–2.9%) with the deepest Spitzer images available.
• 3.  
  The new calibration also tightened the dispersion between the resulting galaxy counts when compared to the brighter galaxy counts in the ground-based+HST+Spitzer filters and improved our estimates of the integrated galaxy light (Section 4), although these changes are well within the uncertainties quoted in Table 3 and hardly visible in Figures 9–11.
• 4.  
  From the object-free sky-SB measurements in Section 5, we confirm that the ZPs of the NIRCam SW detectors have become ∼10%–20% more sensitive. That is, the ratio of 2022 August to October object-free sky-SB values is typically 1.09–1.22 for the NIRCam SW filters, while this ratio is typically 0.88–0.97 for the LW filters, which have thus become somewhat less sensitive compared to the earlier values. Such ZP changes can affect the colors of discrete objects, which we defer to future papers.
• 5.  
  The new calibration has improved the rms variation between the sky-SB measurements compared to our earlier calibrations. Specifically, the relative errors on the object-free sky-SB values have significantly improved with the new calibrations, with the 2022 October to August ratios of these relative errors typically being a factor 0.63–0.83 for the SW filters and sometimes less than 0.5 for the LW filters. Stated differently, the new ZP calibrations and flat-field corrections have reduced the dispersion in the sky-SB estimates significantly, especially for the LW filters. As a consequence, our limits on diffuse light in Section 5 and Figure 13 have also improved, especially in the LW filters.
• 6.  
  The more accurate detector-to-detector ZPs also improved the overall quality of our sky-SB model fits in Figure 12, as described in Section 5.2 and Appendix C. In particular, the range in scale factors f by which we needed to multiply our adopted Rigby et al. (2022) straylight in Equation (4) has increased from f = 0.3, 0.6, 0.8, 0.9 in our v0.5 reduction to f ≃ 0.5, 0.8, 1.0, 1.0 in v1, i.e., providing a tighter range in the predicted SL values for our four PEARLS fields. The SL spectrum, which we assumed to be constant other than this factor f, in fact somewhat depends on JWST's pointing direction (J. Rigby 2022, private communication; J. Rigby et al. 2022, in preparation). We refer to this work for an in-depth discussion of the JWST straylight. For the very red thermal Zodiacal SED (see Figure 5 of Rigby et al. 2022), the effective wavelength (Equation A15 of Bessell & Murphy 2012) of the F444W filter shifts from ∼4.4 μm for objects with relatively flat NIR-spectra like stars and galaxies to ∼4.57 μm for the total sky-SB, which we accounted for in Figures 12–13.

In conclusion, we independently confirm the Rigby et al. (2022) and Boyer et al. (2022) flux scale and adopt their suggested 4% uncertainty in the current JWST NIRCam flux scale in our error analysis below. While model dependent (Appendix C), our sky-SB analysis does give a nearly PSF-independent check on the zero-points, to the extent that the sky-SB is estimated in areas largely devoid of bright-object PSF-wings. The results from the new calibrations have also been propagated into the error budgets for our sky-SB values in Appendix B.3. The total uncertainties then are ≲8% for the SW filters and ≲6% for the LW filters (Appendix B.3).

In order to compare our PEARLS NIRCam 0.9–4.5 μm object counts in Figures 9–10 to previous work from ground-based telescopes and Spitzer/IRAC, the NIRCam flux densities may need to be transformed to the same flux scale as used for the filters in those previous surveys. The NIRCam filter wavelengths are similar, but not identical, to those used previously. The most extensive survey and number counts for λ < 3 μm is that of Koushan et al. (2021, and references therein), who combined many data sets. For λ > 3 μm, we use the compilation of Driver et al. (2016a, 2016b), who included number counts for AB ≲ 18 mag from WISE (Jarrett et al. 2017) and for 18 ≲ AB ≲ 26 mag from IRAC (Ashby et al. 2009, 2015).

Table 7 lists the effective wavelengths λ_c of the relevant VISTA/IRAC and WISE filters in which these previous counts were done, as displayed in Figures 9–10. We therefore use the effective wavelengths of the VISTA/IRAC filters in Table 7 as fiducial to compare our NIRCam object counts to. The flux scale of the WISE W1 and W2 filters was already transformed to the flux scale of these fiducial VISTA/IRAC filters. We used the ICRAR filter transform tool (Robotham et al. 2020) ⁸¹ to calculate the corrections needed to bring the NIRCam AB-magnitude scale onto that of the VISTA/IRAC filters that are closest in wavelength. The tool uses the filter and telescope transmission curves folded with the detector QE curves for a large number of facilities to perform numerical integration over a range of SEDs and redshifts, and produces AB-flux scale offsets, ΔAB, and their uncertainties, σ_ΔAB, which we represent as:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{NIRCam}\ \ \mathrm{ABmag} & = & \mathrm{VISTA}/\mathrm{IRAC}\ \mathrm{ABmag}\\ & & +{\rm{\Delta }}\mathrm{AB}\pm {\sigma }_{{\rm{\Delta }}\mathrm{AB}}.\end{array}\end{eqnarray} \tag{ 6 }$$

The transformation requires an assumption for the redshift distribution of the galaxy population at AB ≃ 20–28 mag, for which we used a median redshift of z_med ≃ 1–2 (e.g., Skelton et al. 2014; Inami et al. 2017). The resulting values for ΔAB and σ_ΔAB are given in Table 7. ⁸² For most filters, the uncertainty in the transformation is 3%–6%, i.e., similar to or larger than the actual flux scale correction needed, which is –0.04 to +0.03 mag. This is mainly because of the wide redshift range sampled. Therefore, no AB-mag scale corrections were applied to our PEARLS NIRCam number counts to compare them to the VISTA/IRAC counts. But we do add this uncertainty to our error budget.

Table 7 also lists the cosmic variance uncertainty for our two PEARLS fields used in the deep NIRCam galaxy counts thus far (Sections 2 and 4.5). Assuming both ZP uncertainties and the CV uncertainty are independent, we show the combined total error on the bottom line of Table 7. These are our IGL errors used in Figures 12–13. These errors are likely conservative, since other deep fields from HST, VLT and Spitzer fold into the galaxy counts of Figures 9–10 thereby reducing CV, except at the faint end of the counts at wavelengths ≳2.0 μm, where we only have NIRCam.

In summary, the uncertainty in the JWST NIRCam zero-points is at least 4%, while the uncertainty of transforming the NIRCam AB-mag scale onto the fiducial VISTA/IRAC filters that have been used for galaxy counts at brighter levels is ≲3%–6%. The combined uncertainty to compare counts that were done with slightly different filters systems on different telescopes is thus ∼3%–7%. Magnitude offsets of that size are hardly noticeable over the very wide magnitude range plotted in Figures 9–10. Future improvements in the NIRCam ZPs through further standard star monitoring and more detailed comparison to the fluxes in the fiducial VISTA/IRAC filters can provide a more accurate comparison, and observing more JWST fields will decrease the uncertainty in the counts from cosmic variance.

For the uncertainties in our observed JWST NIRCam sky-SB values we need to consider other error sources than those that apply to the flux-scale errors in our galaxy counts in Appendix B.2. We wish to make an estimate of the absolute sky-SB in our 13 NIRCam filters, and so the main sources of error are different. For details of an assessment of this kind, we refer to Section 4 and Table 5 of Windhorst et al. (2022), where the sources of error in the absolute sky-SB as measured by WFC3/IR were summarized for the F125W–F160W filters. In short, their total errors in the estimated WFC3/IR sky-SB were 3%–4% in these filters, and dominated by the flat-field (≲2%) and ZP errors (≲1.5%). This was through careful tracking of the WFC3/IR performance and its calibration over 12 yr in orbit. At this stage, such errors for JWST are surely less well known, so we estimate our error budget by giving conservative limits to each main component that affects our estimated sky-SB values.

• 1.  
  Algorithm to get LES: with the LES algorithm of Windhorst et al. (2022) and O'Brien et al. (2022), we divided the 2048 × 2048 pixel image from each individual NIRCam detector into 32 × 32 boxes of 64 × 64 pixels and used these to determine the lowest estimated sky (LES) values following the percentile clip method of (O'Brien et al. 2022). From Monte Carlo simulations with realistic object densities and CR distributions, they showed that the LES method gives reliable estimates of the object-free sky-SB, to within 0.4% of the simulated sky-SB, even in the presence of 10% gradients across the field. While the object density in the NIRCam images of Section 4.5 is ∼3× higher to AB ≲ 28.5 mag compared to the HST WFC3/IR image density at its detection limit of AB ≲ 26.5 mag, the NIRCam SW and LW pixel size is also ∼16–4× smaller in area compared to WFC3/IR, respectively, so that at least similar amounts of empty sky are available to the current depth in the NIRCam images to measure object-free LES sky-SB values. Hence we adopt an uncertainty of ≲0.4% of the algorithm itself estimates the sky-SB in the object free areas. When applying this algorithm, we find that it does ignore areas with residual wisps and snowballs well (as it flags those as potential objects with positive flux to be avoided in the sky-SB estimate).
• 2.  
  ZP Uncertainties: The 4% NIRCam ZP uncertainties of Table 7 also apply to the observed sky-SB values of Figure 12. Since we plot all data points at their actual effective NIRCam wavelengths the Transform uncertainties of Table 7 do not apply. Most of the JWST sky-SB comes from the Zodiacal belt at distances ≲3–5 au (e.g., Windhorst et al. 2022), and the IGL is ∼10–70× dimmer than than the ZL (Figure 13). Hence, a ≲9% CV error in the IGL (Section 4.6) is very small compared to these other errors in the total sky-SB estimates.
• 3.  
  Flat-field and residual 1/f and pedestal uncertainties: we verified that the LES algorithm of Section 5 and Windhorst et al. (2022) finds the cleanest regions to estimate the sky-SB in each detector after the 1/f corrections of Appendix A. With the most recent reduction of context file jwst_0995.pmap_filters, the flat field uncertainty has improved to ∼2% compared to the 7%–8% uncertainty in our earlier reductions with context file jwst_0942.pmap_filters (B. Sunnquist 2022, private communication). Since flat-field uncertainties can be a dominant component in our error budget for absolute sky-SB estimates, we check for this as following. We find that the LES sky-SB estimates have a ≲2%–4% variation between the LW detectors in our 2.7–4.5 μm filters, including our TNJ1338 medium-band LW filters. However, these variations increase to ≲4%–7% for between the SW detectors in SW 0.9–2.0 μm filters. This is likely due to the larger number of detectors, some of which still have residual offsets after the flat-fielding and pedestal removal procedure of Section 3.1. Hence, we adopt a 7% uncertainty in the flat-field induced sky-SB estimate for all SW filters, and a 4% uncertainty for all LW filters.
• 4.  
  Bias and dark-current frame subtraction uncertainties: in the 12 yr on-orbit data analyzed for the WFC3/IR detectors, these errors were ≲1% following Section 4 and Table 5 of Windhorst et al. (2022). The NIRCam bias and dark-current levels and their uncertainties listed in Section 3.1 and its websites are also very low, typically ≲1.4%–2.1% of the sky-SB at 3.5–2.0 μm in Table 4, respectively. ⁸³ Hence, uncertainties in the NIRCam dark-current removal are much smaller than the ZP and the flat field plus residual pedestal uncertainties above.

In summary, following the discussion of Windhorst et al. (2022), we will assume that the above errors in estimating the sky-SB are independent. This is justified because the standard stars from which the NIRCam ZP are derived are measured over an area much smaller than the above dominant flat-field/residual pedestal errors. The resulting uncertainty in our combined error on the absolute NIRCam sky-SB is thus ∼6% of the observed sky-SB for the LW filters and ∼8% for the SW filters. We propagate these errors for each filter and PEARLS field into the NIRcam sky-SB estimates of Tables 4–6 and Figure 13.