# 5-25 $\mu$m Galaxy Number Counts from Deep JWST Data

Meredith A. Stone Steward Observatory, University of Arizona, 933 North Cherry
Avenue, Tucson, AZ 85719, USA  Stacey Alberts Steward Observatory, University
of Arizona, 933 North Cherry Avenue, Tucson, AZ 85719, USA  George H. Rieke
Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson,
AZ 85719, USA  Andrew J. Bunker Department of Physics, University of Oxford,
Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK  Jianwei Lyu (吕建伟)
Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson,
AZ 85719, USA  Pablo G. Pérez-González Centro de Astrobiología (CAB), CSIC-
INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain  Irene
Shivaei Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4,
Torrejón de Ardoz, E-28850, Madrid, Spain  Yongda Zhu Steward Observatory,
University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85719, USA

###### Abstract

Galaxy number counts probe the evolution of galaxies over cosmic time, and
serve as a valuable comparison point to theoretical models of galaxy
formation. We present new galaxy number counts in eight photometric bands
between 5 and 25 $\mu$m from the Systematic Mid-infrared Instrument Legacy
Extragalactic Survey (SMILES) and the JWST Advanced Deep Extragalactic Survey
(JADES) deep MIRI parallel, extending to unprecedented depth. By combining our
new MIRI counts with existing data from Spitzer and AKARI, we achieve counts
across 3-5 orders of magnitude in flux in all MIRI bands. Our counts diverge
from predictions from recent semi-analytical models of galaxy formation,
likely owing to their treatment of mid-infrared aromatic features. Finally, we
integrate our combined JWST-Spitzer counts at 8 and 24 $\mu$m to measure the
cosmic infrared background (CIB) light at these wavelengths; our measured CIB
fluxes are consistent with those from previous mid-infrared surveys, but
larger than predicted by models based on TeV blazar data.

††journal: ApJ††facilities: JWST(MIRI)††software: Astropy (Astropy
Collaboration et al., 2013, 2018), Matplotlib (Hunter, 2007), NumPy (van der
Walt et al., 2011), photutils (Bradley et al., 2022), WebbPSF (Perrin et al.,
2014)

##  1 Introduction

The distribution of the number of sources detected in a given filter as a
function of flux, known as source counts or number counts, probes galaxy
evolution on the largest scales and provides a direct observational comparison
to the results of cosmological simulations and semi-analytical models of
galaxy formation and evolution. Unlike luminosity or mass functions, which can
only be constructed with robust redshift constraints (difficult to obtain for
huge samples of galaxies, and with their own uncertainties), deriving number
counts requires only sufficient depth, area, and robust measurements of source
fluxes and is therefore subject to fewer sources of calibration bias when
comparing to simulations.

Number counts at different wavelengths in the mid and far-infrared (IR) in
particular have revealed a great deal about the evolution of galaxies across
cosmic time. The IRAS and ISO missions revealed that the distribution of mid-
and far-IR number counts differs significantly from the optical and near-IR.
While optical and near-IR counts display a “Euclidean” $dN/dS_{\nu}\propto
S_{\nu}^{-2.5}$ slope down to faint fluxes, consistent with a nonevolving
population of spiral and elliptical galaxies, counts at longer mid-infrared
wavelengths display a characteristic “bump” between fluxes of approximately
$0.1$ and $0.5$ mJy before decaying at fainter fluxes, attributed to the
rapidly evolving population of galaxies between the local Universe and $z\sim
2$ (see e.g. Fazio et al., 2004). In the 21st century, mid- and far-IR number
counts have been probed to fainter and fainter limits, and the contributions
of different galaxy populations inferred, with observations from later
generations of IR space telescopes including Spitzer’s IRAC (Fazio et al.,
2004; Treister et al., 2006; Barmby et al., 2008; Ashby et al., 2009, 2015;
Papovich et al., 2016; Lacy et al., 2021), IRS (Teplitz et al., 2011), and
MIPS (Marleau et al., 2004; Papovich et al., 2004; Chary et al., 2004;
Treister et al., 2006; Rocca-Volmerange et al., 2007; Béthermin et al., 2010;
Clements et al., 2011) instruments and the AKARI mission (Pearson et al.,
2010, 2014; Murata et al., 2014; Davidge et al., 2017). Prior to 2022, these
missions had pushed the $80\%$ completeness limit to $<100$ $\mu$Jy at $8$ and
$24$ $\mu$m with Spitzer, and somewhat higher in intermediate bands with
AKARI.

With the launch of the James Webb Space Telescope (JWST), a wealth of new data
will become available to probe number counts in the mid-infrared with MIRI. As
number counts from MIRI data are brand-new and able to reach unprecedented
depth and access sources fainter than ever observed, they provide critical
constraints on the faint end of the distribution, a region that currently
displays the most variation between different cosmological simulations and
semi-analytical models (see e.g. Lacy et al., 2021). These theoretical models
are not yet calibrated to match JWST observations, so JWST data provides an
independent test of model validity.

Already, mid-infrared number counts from MIRI have been published in the
Stephan’s Quintet ERS field (Ling et al., 2022) at 7.7, 10, and 15 $\mu$m and
in a portion of the Cosmic Evolution Early Release Science (CEERS) field in
the Extended Groth Strip (Wu et al., 2023; Kirkpatrick et al., 2023). These
data extend down to $80$% completeness limits ranging from 0.25 $\mu$Jy at 7.7
$\mu$m to 13 $\mu$Jy at 21 $\mu$m, in some filters extending more than two
orders of magnitude deeper than the deepest pre-JWST data available. However,
the small area of these surveys (of order 5 arcmin2 in each field) limited
their ability to connect to previous observations at the bright end and leave
them vulnerable to cosmic variance.

In this work, we present number counts derived from the Systematic Mid-
infrared Instrument Extragalactic Survey (SMILES) in eight MIRI bands,
supplemented by ultra-deep MIRI parallel imaging at 7.7$\,\mu$m observed as
part of the JWST Advanced Deep Extragalactic Survey (JADES). The wider area
($\sim 34$ arcmin2) and complete mid-infrared wavelength coverage of SMILES
complement the incredibly deep data of the MIRI JADES parallel and previously
existing Spitzer data, allowing us to probe infrared emission across more than
five orders of magnitude in flux and suppress cosmic variance by linking MIRI
counts to those from wider Spitzer surveys. We introduce our observations in
Section 2, outlining our data reduction and completeness correction processes
and photometric accuracy tests performed. We describe the construction of the
number counts and present our results in Section 3. In Section 4, we compare
to a number of state-of-the-art theoretical models of galaxy formation and
evolution. We also integrate our 7.7 and 21 $\mu$m counts to constrain the
cosmic infrared background light and compare the result with values deduced
from TeV gamma ray observations. We summarize in Section 5.

##  2 Data

The primary datasets for this work are MIRI imaging in the GOODS-S field from
the Systematic Mid-infrared Instrument Legacy Extragalactic Survey (SMILES,
PID 1207) at 5.6 - 25.5 $\mu$m (G. Rieke et al., in prep, Alberts et al.,
2024) and the ultra-deep 7.7 $\mu$m parallel (S. Alberts, et al., in prep)
from the JWST Advanced Deep Extragalactic Survey (JADES, PID 1180; Eisenstein
et al., 2023). SMILES is comprised of 15 pointings ($\sim 34$ sq arcmin) of
moderate depth exposures (10.7 - 36.4 min) in 8 MIRI filters (excluding only
F1130W), reaching $5\sigma$ point source sensitivities of 0.20, 0.19, 0.38,
0.59, 0.68, 1.7, 2.8 and 16 $\mu$Jy for F560W, F770W, F1000W, F1280W, F1500W,
F1800W, F2100W, and F2550W respectively, measured in apertures containing 65%
of the encircled energy for each PSF and then aperture corrected (see Section
2.2).

The JADES deep MIRI parallel is four pointings ($\sim 10$ sq arcmin) totalling
43 hours of science time per pointing, reaching a $5\sigma$ point source
sensitivity of 21 nJy in F770W. The $5\sigma$ point source sensitivities are
shown in Table 1 and are generally 2-2.5x lower than predicted by ETC v3.0 up
to F1800W.

The majority of the area of these two surveys fall on deep and medium JADES
NIRCam imaging (Eisenstein et al., 2023) as well as the rich ancillary data in
GOODS-S from facilities such as the Hubble Space Telescope (HST, e.g.
Giavalisco et al., 2004; Grogin et al., 2011).

###  2.1 Data Reduction

The data reduction for SMILES and the JADES MIRI parallel are broadly similar;
a full description can be found in Alberts et al. (2024) Both datasets were
reduced using v.1.12.5 of the JWST calibration pipeline (Bushouse et al.,
2023) using JWST Calibration Reference System (CRDS) context jwst_1188.pmap.

We supplement the standard pipeline with custom, external routines to remove
warm pixels missed by the pipeline, “super” background subtraction (Pérez-
González et al., 2024), and apply astrometry corrections, matched to GAIA
through the JADES NIRCam catalog (Rieke et al., 2023). The final astrometry is
accurate to $0.1-0.2\arcsec$ ($1\sigma$, per axis) for F560W-F2100W and
$0.4\arcsec$ for F2550W, as imaging in this filter is relatively shallow
compared to the rest of the survey, with fewer detections suitable for
deriving an astrometric correction (Alberts et al., 2024).

###  2.2 Source Extraction and Photometry

Blind MIRI photometric catalogs are created with a modified version of the
JADES photometric pipeline (Rieke et al., 2023, Robertson et al., in prep).
Object detection is done using a stacked F560W and F770W SNR image for SMILES
and F770W only for the JADES parallel; these detection images are then used to
create a robust segmentation map through an iterative process that starts at a
low threshold for detection and then performs deblending, cleaning to remove
spurious noise spikes, and faint source detection.

Circular aperture photometry (in apertures of diameter $0.5,0.6,0.7$, and
$1.0\arcsec$) is then measured from the final, clean segmentation map using
photutils with the windowed centroiding algorithm from Source Extractor. Kron
photometry (2.5x scaled) is also measured using a stacked signal map to define
source centroids. Kron apertures smaller than our smallest circular aperture
($d=0.5\arcsec$) are replaced with that circular aperture photometry. SMILES
photometry in filters longer than F770W is measured using the apertures
defined by the stacked F560W+F770W detection image. We assess whether this
approach misses long wavelength-only detections and/or changes the derived
long wavelength flux densities; we find that these effects are small and that
our photometric catalog is robust (Alberts et al., 2024).

Aperture corrections are applied to both the circular aperture and Kron
aperture fluxes using custom point spread functions (PSFs). We derive the
aperture corrections for the latter by measuring the fraction of the PSF flux
which falls outside each Kron aperture, to account for the size and shape of
the PSF. For F1000W-F2550W, PSFs are created using WebbPSF (Perrin et al.,
2014). For F560W and F770W, the PSF modeling is complicated by the cross-like
imaging artifact known as the “cruciform” (Gáspár et al., 2021) and so
empirical PSFs are constructed using high dynamic range imaging of stars
obtained during JWST commissioning (A. Gáspár, private communication).
Finally, photometric uncertainties are obtained using randomly placed
apertures across the mosaics, accounting for local exposure time and
correlated pixel noise (e.g. Whitaker et al., 2011; Rieke et al., 2023).

The final SMILES photometric catalog has 3,591 sources with SNR$\,>4$ in
either the F560W or F770W filter. The JADES MIRI parallel catalog has 2,860
sources with SNR$\,>4$ in F770W.

###  2.3 Completeness and Photometric Accuracy

We measure the completeness $-$ the probability of measuring a source at a
given flux density given the mosaic noise properties $-$ in our images by
inserting and measuring the recovery rate of artificial point sources
constructed from our PSFs (Section 2.2). We derive the completeness using
point sources because faint galaxies are predominantly intrinsically small and
un- or marginally-resolved by MIRI, while MIRI-resolved sources tend to be
bright. Bright extended sources are in a flux regime where the completeness
corrections are negligible (see Figure 1), and we do not expect faint,
extended sources (for which our completeness will be overestimated) to be
numerous enough to significantly affect our counts.

For SMILES, artificial sources are inserted in 43 flux bins with a width of
0.1 dex ranging from 0.01-200 $\mu$Jy. Within each bin, the flux of each
artificial source is chosen at random. For each filter and flux bin, we
generate 10 mosaics with 150 artificial sources ($<5\%$ of the real source
catalog) inserted into random positions (with masking of the map edges),
totaling 1,500 artificial sources per filter per flux bin: a number that
ensures robust Gaussian statistics while limiting computational time
requirements. Photometric catalogs are then created for each altered mosaic
following the procedure outlined in Section 2.2. Artificial sources are
considered recovered if they are found within $0.2\arcsec$ of the inserted
position and the measured flux density is detected at $>2\sigma$ and within
$50\%$ of the input flux. The latter accounts for flux boosting of faint
sources from brighter neighbors under the approximation that galaxies are
randomly distributed. As galaxies are in reality clustered, this is an
underestimation, but we do not expect this to have a significant impact on our
measured number counts given the small MIRI beam (Viero et al., 2013).

This procedure is repeated for the JADES MIRI parallel with 53 flux bins
ranging from 0.001 - 200$\mu$Jy with 0.1 dex binwidth. Given the smaller area
of the parallel, we inject 50 artificial sources at a time into 30 mosaics,
again producing 1,500 artificial sources per flux bin. The resulting
completeness curves for both surveys are shown in Figure 1 and the $80\%$
completeness level is reported in Table 1. A completeness of $80\%$ in each
band is reached at approximately the $5\sigma$ detection limit. This high
completeness at faint fluxes is attributable to our careful background
subtraction (Section 2.1). For the rest of this work, only sources above the
$80\%$ completeness limit will be considered.

Table 1: 80% completeness limits from SMILES and JADES in the MIRI bands.  Band |  $5\sigma$ Sensitivity | 80% Completeness | Area
---|---|---|---
| ($\mu$Jy) | Limit ($\mu$Jy) | (arcmin2)
F560W | 0.20 | 0.20 | 34.8
F770W (JADES) | 0.040 | 0.036 | 10.1
F770W | 0.19 | 0.18 | 34.5
F1000W | 0.38 | 0.36 | 34.3
F1280W | 0.59 | 0.54 | 34.5
F1500W | 0.68 | 0.63 | 34.5
F1800W | 1.7 | 1.6 | 34.5
F2100W | 2.8 | 2.2 | 34.5
F2550W | 16 | 15 | 34.5
Figure 1: Completeness fraction as a function of flux for our MIRI data,
including the JADES parallels F770W data (purple, leftmost points) and the
eight bands observed in the SMILES survey. The 80% completeness threshold,
which is the lower limit to which we report our number counts, is shown as a
gray dashed line.

###  2.4 Contamination by Spurious Sources

We expect some fraction of our detected sources to be spurious: not true
astrophysical sources, but peaks in the noise. The contribution of these
“false positives” to our number counts is certainly insignificant at the
bright end of the distribution, but may become important to the total counts
at the faint end.

To quantify the fraction of spurious sources our counts contain, we perform a
source extraction on negative images (see Papovich et al., 2004). By inverting
the images, the flux of all real sources becomes negative, and the only
positive sources remaining in the image are random Gaussian peaks from noise.
As the noise is random, the number of sources in the negative image should
also reflect the number of spurious sources that will be detected in the
original images and removed from the counts.

We invert the mosaic in each filter and mask the noisy edges before performing
source detection (using Photutils’ find$\\_$peaks algorithm), extracting
photometry with a $1.2\arcsec$-diameter aperture. We find that the vast
majority of spurious sources have integrated fluxes well below our
completeness limits, with only a few noise peaks bright enough to contribute
to the faint end of the counts. Normalizing by the areas of the SMILES and
JADES surveys, we find that we expect $\sim$0.05 to 0.25 counts/arcmin2 from
spurious sources at the 80$\%$ completeness limits of SMILES, and even fewer
in the JADES parallel. Contributions from spurious sources at brighter fluxes
are, as expected, completely negligible, and we therefore make no corrections
accounting for their presence.

##  3 Results

###  3.1 Constructing Number Counts

Differential and cumulative number counts are derived from the photometric
catalogs described in Section 2.2. Precise individual areas are measured for
each mosaic by counting unmasked pixels (Table 1). The exposure times are
largely uniform across the mosaics: $70-80\%$ of pixels have exposure times
within $5\%$ of the median.

We consider counts in each filter down to the $80\%$ completeness limits,
which are comparable to the $5\sigma$ point source sensitivities (Table 1).
The number counts discussed in the following sections are completeness-
corrected, and the error bars reflect pureply statistical error.

###  3.2 Number Counts at 7.7$\mu$m

Figure 2: Cumulative (left) and Euclidean-normalized differential (right)
number counts at 7.7 (MIRI) and 8 (IRAC) $\mu$m. The deep MIRI F770W counts
from the JADES and SMILES surveys are shown as purple and blue circles,
respectively, and their 80% completeness limits are marked by colored, dotted
vertical lines. For visual clarity, points at the bright end with high scatter
and large error bars are plotted with greater transparency. We also plot the
differential F770W counts in two MIRI pointings in the CEERS field measured by
Wu et al. (2023) (black crosses in right panel), which agree very well with
our data. We link our MIRI data to the bright end by including the IRAC 8
$\mu$m counts from the Spitzer Deep Wide-Field Survey (SDWFS; Ashby et al.,
2009) as gray squares: together these produce number counts spanning more than
five orders of magnitude in flux. The Ashby et al. (2009) IRAC counts have had
the contribution from Galactic stars subtracted according to predictions from
Fazio et al. (2004). We also include for comparison predicted number counts at
8 $\mu$m from the GALFORM semi-analytical model (Cowley et al., 2018, gray
solid line), the SHARK semi-analytical model (Lagos et al., 2019, gray dashed
line), the SPRITZ simulations (Bisigello et al., 2021, gray dot-dashed line),
and the models of Rowan-Robinson (2009, 2023) (gray dotted line). For ease of
comparison, we normalize each model to the observed counts at 102 $\mu$Jy.

Our deepest data is at 7.7 $\mu$m. At low redshift, F770W traces the mid-
infrared PAH features, and probes the rest-frame peak of stellar emission at
higher redshift, $z\sim 2-5$. At all redshifts, 7.7 $\mu$m—and other mid-
infrared—counts will additionally contain a population of dust-obscured AGN
and composite (star-forming + AGN) galaxies. When corrected for the
contribution of AGN, 7.7 $\mu$m counts can therefore test the distribution of
stellar mass, particularly at high redshift, of galaxy formation models.

The similarity between the MIRI F770W and Spitzer/IRAC Channel 4 (8 $\mu$m)
filters also allows us to connect our very deep counts with those from wide-
field Spitzer surveys, probing both the bright and faint ends of the
distribution. However, the two filters are not identical. While the MIRI F770W
filter has slightly higher throughput, the IRAC-4 bandpass is wider on the red
side, extending to $\sim 9.5$ $\mu$m compared to F770W’s $\sim 8.75$ $\mu$m.
The largest difference between an individual source’s IRAC-4 and F770W fluxes
will likely occur for sources at redshifts where a bright PAH feature is
shifted out of the F770W bandpass but remains in the IRAC-4 bandpass, but this
occurs in only two small redshift ranges ($z\sim 0.13-0.27$ and $0.41-0.52$
for the 7.7 and 6.2 $\mu$m PAH features respectively). Early JWST results
comparing MIRI to IRAC counts indicate that the effect of this difference is
very small (see e.g. Yang et al., 2023; Wu et al., 2023) in the flux regime
where MIRI and IRAC counts overlap.

We present our differential and cumulative 7.7 $\mu$m counts from SMILES and
JADES in Figure 2, alongside bright-end IRAC-4 data from the Spitzer Deep
Wide-Field Survey (SDWFS) (Ashby et al., 2009). Together, these deep MIRI data
and wide Spitzer data provide the deepest 8 $\mu$m counts ever observed,
spanning more than five orders of magnitude in flux from $<0.1$ to more than
$10^{4}$ $\mu$Jy.

The MIRI F770W counts in the SMILES and JADES fields are slightly offset from
each other. We have verified that this offset is not due to any differences in
the data reduction process. Instead, this is likely due to cosmic variance:
the Hubble Ultra Deep Field (HUDF, which forms the center of the SMILES
footprint) is known to be slightly underdense, while a number of overdensities
have already been catalogued in the small JADES parallel field (Alberts et
al., 2023). The difference between the SMILES and JADES counts likely encodes
the typical variation that will be observed between surveys (see also the
variation in counts for fluxes probed by Spitzer, Béthermin et al., 2010).

The differential counts (right panel of Figure 2), which have been normalized
to a Euclidean universe (multiplied by $S_{\nu}^{2.5}$), display the
characteristic power-law decline toward faint fluxes observed by other
telescopes, and flatten at the bright end when we remove the contribution from
Galactic stars from the model of Fazio et al. (2004).

For comparison, we also plot in Figure 2 the 8 $\mu$m number count predictions
from a number of simulations and semi-analytical models of galaxy formation,
namely the GALFORM semi-analytical model and GRASIL radiative transfer code
(F770W, Cowley et al., 2018), the SHARK semi-analytical model (IRAC-4, Lagos
et al., 2019), the SPRITZ simulations (F770W, Bisigello et al., 2021), and the
IRAC-4 model of Rowan-Robinson (2009, 2023). We discuss these models in
greater depth in Section 4.1.

###  3.3 Number Counts at 21$\mu$m and $25\mu$m

Figure 3: Number counts at 21 and 25 $\mu$m. We plot our cumulative (left) and
Euclidean-normalized differential (right) MIRI SMILES F2100W and F2550W counts
as colored circles (lighter red and darker brown, respectively). As in Figure
2, points at the bright end with high scatter and large error bars are plotted
with greater transparency for visual clarity. Colored, dotted vertical lines
mark the 80% completeness limits of the SMILES 21 and 25 $\mu$m counts. We
link to the bright end by including Spitzer/MIPS 24 $\mu$m counts from
Papovich et al. (2004) (gray squares). These counts are not significantly
contaminated by stars. For comparison, we also plot the differential F2100W
counts in a portion of the CEERS field from Wu et al. (2023) at right as black
crosses. We additionally include modeled number counts from the literature,
both from the GALFORM semi-analytical model (Cowley et al., 2018, gray solid
line) and the models of Rowan-Robinson (2009, 2023) (gray dotted line). These
models have been normalized to the counts at 103 $\mu$Jy.

We plot our SMILES cumulative and differential counts at both 21 and 25 $\mu$m
in Figure 3, as the counts in both bands have very similar distributions. Our
new JWST counts also align with counts from deep 24 $\mu$m Spitzer surveys
(Papovich et al., 2004; Béthermin et al., 2010, plotted as gray squares and
black plus-signs respectively), which we use to extend our distribution to the
bright end. Together, these 21/24 $\mu$m counts extend over approximately four
orders of magnitude, from less than $100$ to more than $10^{5}$ $\mu$Jy, and
probe deeper than any $\sim$24 $\mu$m counts previously published.
Additionally, due to the unprecedented spatial resolution of JWST, our MIRI
$\sim$24 $\mu$m counts are not confusion-limited, unlike the Spitzer counts
also shown in Figure 3.

The flux at $\sim 25$ $\mu$m probes the mid-infrared warm dust continuum until
$z\sim 1.5$, when the bright 7.7 $\mu$m PAH feature begins to shift into the
band. We observe the characteristic bump at $\sim 100$ $\mu$Jy in the
Euclidean-normalized differential counts. Despite the minimal overlap of the
MIRI F2100W and F2550W bands, the counts display very similar distributions
down to the limit of the F2550W data.

The MIPS 24 $\mu$m filter is much broader than the MIRI filters, encompassing
most of the wavelength range covered by the MIRI F2100W and F2550W filters.
This discrepancy, however, does not appear to severely affect the magnitude or
shape of the source counts’ distribution, as the MIPS 24 $\mu$m data from
Papovich et al. (2004) and Béthermin et al. (2010) display fair agreement,
within the uncertainties and scatter, with the MIRI data. No correction has
been applied to either the MIRI or MIPS counts to account for the very
different bandpasses and bring the counts to alignment.

In addition to the Spitzer counts, we can compare our results to the counts
obtained in two MIRI pointings of the CEERS survey by Wu et al. (2023). While
the Wu et al. (2023) counts agree well with ours down to $\sim 100$ $\mu$Jy,
their differential counts flatten at the faintest end of their data, between
$\sim 10$ and 100 $\mu$Jy, rather than continuing to decrease as a power-law
like our SMILES counts or model predictions. This may be an artifact of those
data’s lack of background subtraction, as we observe a similar effect to a
lesser extent in other bands longward of 10 $\mu$m (discussed below in Section
3.4).

In general, when normalized near 103 $\mu$Jy, the models of Rowan-Robinson
(2009, 2023) reasonably closely follow our observed counts down to the faint
end, and reproduce the shape of the differential number counts distribution
well (right panel of Figure 3). The GALFORM model, on the other hand,
underpredicts the number of faint sources and does not exhibit the peak in the
differential counts near 100 $\mu$Jy seen in both JWST and Spitzer data.
Further discussion of these comparisons can be found in Section 4.1.

###  3.4 Number Counts in other Bands

Figure 4: Cumulative number counts for the remaining MIRI bands. MIRI data
from the SMILES survey is plotted as colored circles, with colored dotted
vertical lines marking the 80% completeness limit in each band. We also
include model number counts in each band from the GALFORM semi-analytical
model (Cowley et al., 2018) as gray solid lines. All models have been
normalized to the counts near the bright end (102 $\mu$Jy at 5.6, 10, and 12.8
$\mu$m; 103 $\mu$Jy at 15 and 18 $\mu$m). At 5.6 and 10 $\mu$m, we include the
predictions of the SPRITZ simulations (Bisigello et al., 2021, gray dot-dashed
line). We extend the counts to slightly brighter fluxes by including data from
the Spitzer Deep Wide-Field Survey at 5.6 $\mu$m (Ashby et al., 2009, gray
squares) and the AKARI NEP-Wide (Pearson et al., 2010) and NEP-Deep surveys
(Pearson et al., 2014) at 15 and 18 $\mu$m respectively (gray stars). We have
corrected the IRAC 5.6 $\mu$m counts for the contribution from Galactic stars
(Fazio et al., 2004). At 10, 12.8, 15, and 18 $\mu$m, we include MIRI number
counts in two MIRI pointings of the CEERS field as black crosses (Wu et al.,
2023). At 10 $\mu$m, counts from more extensive data in the CEERS survey are
shown as gray triangles (Kirkpatrick et al., 2023).

We plot the cumulative number counts for all remaining MIRI bands (F560W,
F1000W, F1280W, F1500W, and F1800W) in Figure 4. We include, when available,
other data and models in the bands, including AKARI counts (Pearson et al.,
2010, 2014), MIRI counts in the CEERS field (Wu et al., 2023; Kirkpatrick et
al., 2023), and the GALFORM (Cowley et al., 2018), SPRITZ (Bisigello et al.,
2021), and Rowan-Robinson (2009, 2023) models.

Counts from the AKARI North Ecliptic Pole survey at 15 and 18 $\mu$m (Pearson
et al., 2010, 2014) extend our MIRI counts to the bright end, and are in very
good agreement with our counts where they overlap. As a result, the counts at
15 and 18 $\mu$m extend over nearly four orders of magnitude in flux. These
AKARI counts have been corrected for the contributions of stars, but not for
the difference between the MIRI F1500W/F1800W and AKARI IRC $L15$/$L18$
bandpasses. The AKARI bandpasses are significantly broader than the
corresponding MIRI bandpasses: for example, while the two filters have roughly
the same central wavelength, $L15$’s half-power width (6.8 $\mu$m) is more
than twice that of F1500W ($\sim 3.1$ $\mu$m). Nonetheless, the excellent
agreement of the AKARI counts with our MIRI data demonstrates that the
different filter shapes do not strongly affect the counts.

At the fainter end, we can compare our counts to the slightly shallower counts
measured in two of the MIRI pointings of the CEERS field by Wu et al. (2023).
These counts, in F1000W, F1280W, F1500W, and F1800W, drop off quickly at $\sim
100$ mJy due to the smaller area of the survey, and do not go quite as deep.
While they agree fairly well with our counts, some of the bands (particularly
F1000W and F1800W) display steeper slopes between 1 and 100 $\mu$Jy. This may
be a background subtraction artifact, similar to that observed near the faint
end at 21 $\mu$m (right panel of Figure 3). This hypothesis is reinforced by
the 10 $\mu$m counts from Kirkpatrick et al. (2023): these data, also from the
CEERS survey, agree well with our 10 $\mu$m counts.

###  3.5 Cosmic Variance

The error bars in Figures 2, 3, and 4 reflect only statistical errors.
Potentially larger errors occur because of cosmic variance: the differences in
source characteristics among different fields on the sky (e.g., Trenti &
Stiavelli, 2008). Cosmic variance is reduced significantly as the number of
distinct fields surveyed increases, and also with their area. In addition, the
effects of cosmic variance are more pronounced for the brighter sources than
the fainter ones (e.g., Robertson, 2010; Driver & Robotham, 2010; Moster et
al., 2011).

By combining counts from the large-area, moderate-depth surveys of Spitzer and
AKARI with the (to date) small areas surveyed to great depth by JWST, the
cosmic variance effects at play can be understood and mitigated. We discuss
cosmic variance for each wavelength in turn.

5.6 $\mu$m: The top left panel of Figure 4 compares our 5.6 $\mu$m counts to
those in IRAC-3 reported in $\sim$ eight square degrees in the Bootes field by
Ashby et al. (2009). The Ashby et al. (2009) counts demonstrate excellent
agreement with those in the 0.38 square degree survey in the Extended Groth
Strip surveyed to $\sim$ five times deeper reported by Barmby et al. (2008),
and with counts from the IRAC Shallow Survey (Fazio et al., 2004). The large
areas surveyed and the agreement of counts from two widely separated fields
indicate that cosmic variance does not have a strong effect on these counts.
Their close agreement with our JWST-based counts in the region of overlap
indicates, despite the small area surveyed, that the latter counts are also
representative without a large bias due to cosmic variance.

7.7 $\mu$m: Figure 2 compares our counts at 7.7 $\mu$m to those from Ashby et
al. (2009) at 7.8 $\mu$m. As at 5.6 $\mu$m, the Ashby et al. (2009) 7.8 $\mu$m
counts agree well with counts measured in other fields with Spitzer (Barmby et
al., 2008; Fazio et al., 2004), and agree with our MIRI F770W counts where
they overlap. We therefore do not expect cosmic variance to seriously bias our
7.7 $\mu$m counts.

15 & 18 $\mu$m: The counts in Figure 4 at 15 and 18 $\mu$m use results from
the AKARI North Ecliptic Pole (NEP) surveys111ISO counts at 15 $\mu$m were
reported by Oliver et al. (1997) but have been supplanted by those from
AKARI.. These observations were designed to reduce cosmic variance to the
$5-10$% level (Pearson et al., 2010)222The coverage of even larger (and more)
fields in IRAC Band 4 and MIPS 24 $\mu$m reinforces our expectation that these
surveys are not significantly affected by cosmic variance.. They cover areas
of $\sim$ 0.4 - 0.6 deg2 and 5.8 deg2 to depths of 117 (150) and 250 (300)
$\mu$Jy respectively at 15 (18) $\mu$m (Pearson et al., 2010, 2014). The
excellent agreement with the counts in the SMILES field in the region of
overlap again implies that the entire range shown in Figure 4 is unlikely to
be compromised by cosmic variance. On the other hand, MIRI counts from Wu et
al. (2023) fall low at the bright end, suggestive of an influence of variance
on those results due possibly to the significantly smaller field available
($\sim$ 6 arcmin2).

21 - 25 $\mu$m: Figure 3 compares JWST/MIRI counts at 21 and 25.5 $\mu$m with
Spitzer/MIPS counts at 24 $\mu$m. The counts from Papovich et al. (2004) are
from five widely separated fields with $\sim$ 0.61 deg2 surveyed to a depth of
$\sim$ 80 $\mu$Jy and $>$ 10 deg2 to a depth of at least 270 $\mu$Jy.
Béthermin et al. (2010) report surveys of $\sim$ 0.64 deg2 to $\sim$ 80
$\mu$Jy and 53 deg2 to at least $\sim$ 400 $\mu$Jy in seven additional fields.
These surveys agree very well with each other and with our counts down to
their limits of 80 $\mu$Jy (see Figure 3). Again, this indicates that the full
range of counts presented here are not strongly affected by cosmic variance.

Table 2: Cumulative MIRI number counts from SMILES and JADES

Flux |  | F770W | F770W |  |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
($\mu$Jy) | F560W | (SMILES) | (JADES) | F1000W | F1280W | F1500W | F1800W | F2100W | F2550W
3.98e-2 | - | - | 218 $\pm$ 5 | - | - | - | - | - | -
6.31e-2 | - | - | 167 $\pm$ 4 | - | - | - | - | - | -
1.00e-1 | - | - | 130 $\pm$ 4 | - | - | - | - | - | -
1.58e-1 | - | - | 104 $\pm$ 3 | - | - | - | - | - | -
2.51e-1 | 71.3 $\pm$ 1.4 | 55.0 $\pm$ 1.3 | 81.8 $\pm$ 2.9 | - | - | - | - | - | -
3.98e-1 | 52.2 $\pm$ 1.2 | 41.7 $\pm$ 1.1 | 64.5 $\pm$ 2.6 | 41.3 $\pm$ 1.1 | - | - | - | - | -
6.31e-1 | 38.7 $\pm$ 1.1 | 32.6 $\pm$ 1.0 | 50.8 $\pm$ 2.3 | 31.1 $\pm$ 1.0 | 35.6 $\pm$ 1.0 | 39.3 $\pm$ 1.0 | - | - | -
1.00e+0 | 30.0 $\pm$ 0.9 | 25.3 $\pm$ 0.9 | 41.2 $\pm$ 2.0 | 25.0 $\pm$ 0.9 | 27.1 $\pm$ 0.9 | 30.0 $\pm$ 0.9 | - | - | -
1.58e+0 | 23.0 $\pm$ 0.8 | 20.6 $\pm$ 0.8 | 31.6 $\pm$ 1.8 | 20.0 $\pm$ 0.8 | 22.0 $\pm$ 0.8 | 24.1 $\pm$ 0.8 | - | - | -
2.51e+0 | 18.0 $\pm$ 0.7 | 16.1 $\pm$ 0.7 | 23.6 $\pm$ 1.6 | 15.8 $\pm$ 0.7 | 18.0 $\pm$ 0.7 | 19.6 $\pm$ 0.8 | 22.7 $\pm$ 0.8 | 27.7 $\pm$ 0.9 | -
3.98e+0 | 13.7 $\pm$ 0.6 | 11.7 $\pm$ 0.6 | 17.6 $\pm$ 1.3 | 12.1 $\pm$ 0.6 | 14.5 $\pm$ 0.6 | 15.6 $\pm$ 0.7 | 18.5 $\pm$ 0.7 | 21.0 $\pm$ 0.8 | -
6.31e+0 | 9.73 $\pm$ 0.53 | 8.48 $\pm$ 0.50 | 12.5 $\pm$ 1.1 | 8.81 $\pm$ 0.51 | 11.2 $\pm$ 0.6 | 12.7 $\pm$ 0.6 | 15.0 $\pm$ 0.7 | 16.4 $\pm$ 0.7 | -
1.00e+1 | 6.60 $\pm$ 0.44 | 5.41 $\pm$ 0.40 | 8.01 $\pm$ 0.90 | 5.96 $\pm$ 0.42 | 8.67 $\pm$ 0.50 | 9.96 $\pm$ 0.54 | 12.6 $\pm$ 0.6 | 13.5 $\pm$ 0.6 | -
1.58e+1 | 4.37 $\pm$ 0.35 | 3.17 $\pm$ 0.30 | 4.62 $\pm$ 0.69 | 4.16 $\pm$ 0.35 | 6.17 $\pm$ 0.42 | 7.81 $\pm$ 0.48 | 9.58 $\pm$ 0.53 | 10.7 $\pm$ 0.6 | 13.5 $\pm$ 0.6
2.51e+1 | 2.46 $\pm$ 0.27 | 1.79 $\pm$ 0.23 | 2.06 $\pm$ 0.46 | 2.92 $\pm$ 0.29 | 4.79 $\pm$ 0.37 | 5.99 $\pm$ 0.42 | 7.26 $\pm$ 0.46 | 8.29 $\pm$ 0.49 | 9.13 $\pm$ 0.51
3.98e+1 | 1.27 $\pm$ 0.19 | 0.85 $\pm$ 0.16 | 1.13 $\pm$ 0.34 | 1.86 $\pm$ 0.23 | 3.58 $\pm$ 0.32 | 4.41 $\pm$ 0.36 | 5.58 $\pm$ 0.40 | 6.23 $\pm$ 0.43 | 6.50 $\pm$ 0.43
6.31e+1 | 0.46 $\pm$ 0.12 | 0.35 $\pm$ 0.10 | 0.31 $\pm$ 0.18 | 1.27 $\pm$ 0.19 | 2.53 $\pm$ 0.27 | 3.03 $\pm$ 0.30 | 4.03 $\pm$ 0.34 | 4.12 $\pm$ 0.35 | 4.41 $\pm$ 0.36
1.00e+2 | 0.12 $\pm$ 0.06 | 0.21 $\pm$ 0.08 | 0.10 $\pm$ 0.10 | 0.83 $\pm$ 0.16 | 1.73 $\pm$ 0.22 | 1.91 $\pm$ 0.24 | 2.58 $\pm$ 0.27 | 2.79 $\pm$ 0.28 | 2.70 $\pm$ 0.28
1.58e+2 | 0.06 $\pm$ 0.04 | 0.18 $\pm$ 0.07 | - | 0.53 $\pm$ 0.12 | 0.91 $\pm$ 0.16 | 0.94 $\pm$ 0.17 | 1.56 $\pm$ 0.21 | 1.47 $\pm$ 0.21 | 1.53 $\pm$ 0.21
2.51e+2 | - | 0.09 $\pm$ 0.05 | - | 0.24 $\pm$ 0.08 | 0.56 $\pm$ 0.13 | 0.53 $\pm$ 0.12 | 0.68 $\pm$ 0.14 | 0.79 $\pm$ 0.15 | 0.76 $\pm$ 0.15
3.98e+2 | - | 0.06 $\pm$ 0.04 | - | 0.21 $\pm$ 0.08 | 0.26 $\pm$ 0.09 | 0.21 $\pm$ 0.08 | 0.29 $\pm$ 0.09 | 0.35 $\pm$ 0.10 | 0.32 $\pm$ 0.10
6.31e+2 | - | 0.03 $\pm$ 0.03 | - | 0.12 $\pm$ 0.06 | 0.09 $\pm$ 0.05 | 0.09 $\pm$ 0.05 | 0.12 $\pm$ 0.06 | 0.09 $\pm$ 0.05 | 0.18 $\pm$ 0.07
1.00e+3 | - | - | - | 0.06 $\pm$ 0.04 | 0.06 $\pm$ 0.04 | 0.03 $\pm$ 0.03 | 0.03 $\pm$ 0.03 | 0.03 $\pm$ 0.03 | 0.06 $\pm$ 0.04

Note. — All reported counts (2-10) are cumulative (N $>$ $S_{\nu}$), with
units arcmin-2.

##  4 Discussion

###  4.1 Models

Figures 2 - 4 show the predictions of number count models between 5 and 25
$\mu$m, many of which diverge from observed counts. Many models, even those
published recently, use SED templates created before Spitzer IRS data were
available and/or use relatively abstract templates that do not match the
results with full Spitzer data (e.g., as in Rieke et al., 2008; Magdis et al.,
2012; Schreiber et al., 2018; Bernhard et al., 2021). This could be important
because the aromatic bands dominate the infrared output of star-forming
galaxies from 6 through 12 $\mu$m and fall within the wavelength range of JWST
and Spitzer number counts for $z\sim 0-3$. For example, Cowley et al. (2018)
utilize a combination of models from GRASIL (Silva et al., 1998) and GALFORM
(Granato et al., 2000). However, in these models the PAH bands are represented
very schematically, and in the latter one the region between 15 and 40 $\mu$m
is not represented accurately and the silicate absorption features at 10 and
18 $\mu$m are not included. The modeled number counts by Rowan-Robinson (2023)
are based on mid-infrared SEDs from Rowan-Robinson et al. (2004, 2005). These
are constrained by Spitzer photometry, but the PAH bands are from ISO
spectroscopy of a very small number of galaxies (Rowan-Robinson, 2001) and
again are quite schematic. This relatively limited sample of PAH spectra can
cause problems: for an example, compare the proposed spectra in Figure 6 in
Rowan-Robinson et al. (2005) for “Arp 220 starbursts” with the actual spectrum
of Arp 220 from Spoon et al. (2004).

Two sets of models, SHARK (Lagos et al., 2019) and SPRITZ (Bisigello et al.,
2021), do include accurate spectra of the PAH bands making use of Spitzer IRS
spectra. The SHARK models are based on a semi-analytic model of galaxy
formation and evolution. These models are compared with counts in the IRAC
bands, and their good matching of the 7.7 $\mu$m counts presented here to very
faint levels is very encouraging in terms of such models. Obviously, using
this approach to model counts at the longer wavelengths is very desirable,
continuing with accurate template SEDs. The SPRITZ models are empirically
based on a global analysis of galaxy properties. They provide good fits out to
8 $\mu$m, including to our very deep 7.7 $\mu$m counts. At the longer
wavelengths, their results are subject to substantial errors that undermine a
quantitative assessment of their accuracy. For both sets of models, the
agreement with our 7.7 $\mu$m counts demonstrates that representing the PAH
bands accurately is critical to interpreting the number counts presented in
this paper.

###  4.2 The Cosmic Infrared Background

The Cosmic Infrared Background (CIB) is the fraction of the Extragalactic
Background Light (EBL) emitted in the infrared, and arises mainly from dust-
obscured star formation (with contributions from dust-obscured AGN) in
galaxies across cosmic time (Hauser et al., 1998; Gispert et al., 2000;
Lagache et al., 2005). With more sensitive telescopes, more and more of the
CIB can be resolved into individual sources, e.g. the sources making up our
galaxy number counts. Our discussion below follows on from the treatment in
Driver et al. (2016).

To estimate the total 8 $\mu$m flux from sources detectable in our MIRI data,
we combine and integrate 8 $\mu$m number counts from the Spitzer Deep Wide-
Field Survey (SDWFS, Ashby et al., 2009) at the bright end, with stars
removed, and our JADES MIRI parallel at the faint end. To combine the data
sets, we exclude all JADES MIRI parallel counts at fluxes brighter than
$S_{\nu}=50$ $\mu$Jy, because the counts in this flux regime are affected by
cosmic variance due the small area of the survey. The wider-area SDWFS has
data in this flux range, which will be less affected by cosmic variance.
Correspondingly, we exclude all counts at fluxes fainter than $S_{\nu}=30$
$\mu$Jy from the SDWFS IRAC data: this flux corresponds to the 8 $\mu$m SDWFS
80% completeness limit. In the small remaining region of overlap, the MIRI and
IRAC counts agree well (see Figure 2). We also do not include the two
brightest SDWFS IRAC points ($S_{\nu}>2\times 10^{4}$ $\mu$Jy) in our
integration, as they deviate sharply from the flat, Euclidean trend expected
at very bright fluxes in the right panel of Figure 2 due to cosmic variance.
After integrating and propagating through the statistical errors on each bin,
we obtain an 8 $\mu$m CIB flux of $3.02\pm 0.09$ nW m-2 sr-1, slightly larger
than estimates from 8 $\mu$m Spitzer surveys (2.6 nW m-2 sr-1, Fazio et al.,
2004).

Driver et al. (2016) estimated the CIB at 8 $\mu$m using Spitzer/IRAC data
from Barmby et al. (2008), who surveyed 0.38 deg2 in the Extended Groth Strip.
Their uncertainty from cosmic variance ($\pm 0.08$ mag) translates to
approximately $\pm 0.035$ mag for the 10 deg2 area we have used to anchor our
counts at bright levels (Ashby et al., 2009). They overstate the zero point
error; the absolute calibration of the IRAC data is accurate to $\sim$ 3%. The
ultra-deep JWST data basically remove uncertainties due to extrapolation to
low fluxes. The net error in our value should therefore be $\sim$ 10%.

We calculate number counts at 21/24 $\mu$m in the same manner, by combining
Spitzer/MIPS counts (Papovich et al., 2004) and our new MIRI SMILES counts. We
remove the brightest MIRI 21 $\mu$m points ($S_{\nu}>200$ $\mu$Jy) which are
most affected by cosmic variance. The 80% completeness limit for the Papovich
et al. (2004) MIPS data lies at $S_{\nu}\sim 100$ $\mu$Jy, so we do not
include any fainter MIPS counts in our integration. The SMILES/MIRI and MIPS
points agree very well (see Figure 3) near this remaining small region of
overlap. We integrate from the faintest MIRI point (3.25 $\mu$Jy) to the
brightest MIPS point (4.76$\times 10^{4}$ $\mu$Jy), and do not extrapolate;
this integrated value is therefore a lower limit on the 21/24 $\mu$m CIB flux.
Propagating through only the statistical errors on the counts, we obtain a
flux of $2.85\pm 0.04$ nW m-2 sr-1. Since we have taken our bright 24 $\mu$m
counts from the same source as Driver et al. (2016), their estimate of a
cosmic variance error of $\pm 0.07$ also applies to our result. Most other
errors listed in Driver et al. (2016) should be reduced by use of the JWST
data. As before, the zero point error is significantly overstated; the
absolute calibration of MIPS is accurate to $\sim$ 3% (Rieke et al., 2008). We
conclude that the total error in our value is, as at 7.8 $\mu$m, $\sim$10%.

Béthermin et al. (2010) found a CIB flux of 2.29$\pm$0.09 nW m-2 sr-1 from
Spitzer 24 $\mu$m counts probing fluxes $>35\mu$Jy. However, when the counts
were extrapolated at the faint end, they calculated an integrated 24$\mu$m CIB
flux of 2.86${}^{+0.19}_{-0.16}$ nW m-2 sr-1, in excellent agreement with our
value. Dole et al. (2006) also measured the 24 $\mu$m CIB by adding up the
flux from all sources brighter than 60 $\mu$Jy and using the counts of
Papovich et al. (2004) to determine the contribution from fainter sources,
finding an EBL value of 2.7 $\pm$ 0.26 nW m-2 sr-1. Because Dole et al. (2006)
do not quote any additional error from the extrapolation to very faint limits,
this uncertainty is likely somewhat underestimated.

TeV gamma rays are attenuated as they traverse intergalactic space through
pair production with the ambient ultraviolet, optical, and infrared photons.
Semi-analytic models (SAMs) of the history of galaxy formation and galaxy
properties can be combined with TeV gamma ray measurements to predict the
spectrum of the EBL. Our values at both 7.8 and 24 $\mu$m are in mild tension
with the fiducial EBL model based on TeV gamma ray sources by Gilmore et al.
(2012), which predicts values of 2.4 nW m-2 sr-1 at both wavelengths. As shown
by Abeysekara et al. (2019), the highest weight values based on TeV gamma ray
observations tend to fall on or below the Gilmore fiducial, although the
disagreement with our values is not huge. It should be kept in mind that the
EBL estimates from source counts are strictly lower limits, since they exclude
any diffuse emission and discriminate against extended sources. The
discrepancy might therefore be larger than indicated here. These accurate
determinations of the mid-IR EBL suggest that some tuning of the SAMs used to
predict the EBL may be desirable.

##  5 Summary and Conclusions

We measure JWST/MIRI source counts from 5 to 25 $\mu$m using data from the
SMILES and JADES surveys. We combine these unprecedentedly deep counts with
existing, wider data from Spitzer and AKARI to produce counts spanning, in
some bands, five or more orders of magnitude in flux. The uncertainties on our
counts from cosmic variance and spurious sources are small.

When comparing our observed counts to predictions from a number of semi-
analytical models of galaxy formation, we find that the accuracy of the
predicted counts is highly dependent on the templates used to construct model
SEDs. In particular, models incorporating Spitzer/IRS spectroscopy of the PAH
features from 6-12 $\mu$m more closely match our observed counts.

We integrate our deep JWST + wider-area Spitzer source counts at 8 and 24
$\mu$m to constrain the cosmic infrared background at these wavelengths. Our
results are consistent with those from shallower infrared surveys which
extrapolated their data at the faint end, but our deep data significantly
reduces the associated uncertainty on the integrated EBL. Our results are
strictly lower limits since our number counts are not sensitive to diffuse or
extended emission. Nonetheless, we obtain larger EBL fluxes than predicted by
fiducial EBL models based on TeV blazar data, suggesting that some tuning of
the semi-analytic models used for these predictions is needed.

Robust infrared source counts spanning several orders of magnitude in flux can
also be directly compared to predicted counts from models of galaxy formation.
State-of-the-art models have large uncertainties at wavelengths longer than 8
$\mu$m. Our accurate counts should encourage generation of more robust models
for comparison. Such modeling is also timely, since JWST’s unique capabilities
compared to any previous infrared telescope will allow it to continue to
provide critical constraints on models of galaxy formation and extragalactic
light as larger areas of the sky are imaged to unprecedented depth.

MS, SA, GR, JL, and YZ acknowledge support from the JWST Mid-Infrared
Instrument (MIRI) grant 80NSSC18K0555, and the NIRCam science support contract
NAS5-02105, both from NASA Goddard Space Flight Center to the University of
Arizona. The JWST data presented in this paper were obtained from the Mikulski
Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.
The observations can be accessed viahttp://dx.doi.org/10.17909/zfxs-ye40
(catalog DOI.) PGP-G acknowledges support from grant PID2022-139567NB-I00
funded by Spanish Ministerio de Ciencia e Innovación
MCIN/AEI/10.13039/501100011033, FEDER Una manera de hacer Europa. AJB
acknowledges funding from the “FirstGalaxies” Advanced Grant from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (Grant agreement No. 789056).

## References

  * Abeysekara et al. (2019) Abeysekara, A. U., Archer, A., Benbow, W., et al. 2019, ApJ, 885, 150, doi: 10.3847/1538-4357/ab4817
  * Alberts et al. (2023) Alberts, S., Williams, C. C., Helton, J. M., et al. 2023, arXiv e-prints, arXiv:2312.12207, doi: 10.48550/arXiv.2312.12207
  * Alberts et al. (2024) Alberts, S., Lyu, J., Shivaei, I., et al. 2024, arXiv e-prints, arXiv:2405.15972.  https://arxiv.org/abs/2405.15972
  * Ashby et al. (2009) Ashby, M. L. N., Stern, D., Brodwin, M., et al. 2009, ApJ, 701, 428, doi: 10.1088/0004-637X/701/1/428
  * Ashby et al. (2015) Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2015, ApJS, 218, 33, doi: 10.1088/0067-0049/218/2/33
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Barmby et al. (2008) Barmby, P., Huang, J. S., Ashby, M. L. N., et al. 2008, ApJS, 177, 431, doi: 10.1086/588583
  * Bernhard et al. (2021) Bernhard, E., Tadhunter, C., Mullaney, J. R., et al. 2021, MNRAS, 503, 2598, doi: 10.1093/mnras/stab419
  * Béthermin et al. (2010) Béthermin, M., Dole, H., Beelen, A., & Aussel, H. 2010, A&A, 512, A78, doi: 10.1051/0004-6361/200913279
  * Bisigello et al. (2021) Bisigello, L., Gruppioni, C., Feltre, A., et al. 2021, A&A, 651, A52, doi: 10.1051/0004-6361/202039909
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.6.0, 1.6.0, Zenodo, Zenodo, doi: 10.5281/zenodo.7419741
  * Bushouse et al. (2023) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, JWST Calibration Pipeline, 1.12.5, Zenodo, doi: 10.5281/zenodo.10022973
  * Chary et al. (2004) Chary, R., Casertano, S., Dickinson, M. E., et al. 2004, ApJS, 154, 80, doi: 10.1086/423307
  * Clements et al. (2011) Clements, D. L., Bendo, G., Pearson, C., et al. 2011, MNRAS, 411, 373, doi: 10.1111/j.1365-2966.2010.17689.x
  * Cowley et al. (2018) Cowley, W. I., Baugh, C. M., Cole, S., Frenk, C. S., & Lacey, C. G. 2018, MNRAS, 474, 2352, doi: 10.1093/mnras/stx2897
  * Davidge et al. (2017) Davidge, H., Serjeant, S., Pearson, C., et al. 2017, MNRAS, 472, 4259, doi: 10.1093/mnras/stx1935
  * Dole et al. (2006) Dole, H., Lagache, G., Puget, J. L., et al. 2006, A&A, 451, 417, doi: 10.1051/0004-6361:20054446
  * Driver & Robotham (2010) Driver, S. P., & Robotham, A. S. G. 2010, MNRAS, 407, 2131, doi: 10.1111/j.1365-2966.2010.17028.x
  * Driver et al. (2016) Driver, S. P., Andrews, S. K., Davies, L. J., et al. 2016, ApJ, 827, 108, doi: 10.3847/0004-637X/827/2/108
  * Eisenstein et al. (2023) Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv e-prints, arXiv:2306.02465, doi: 10.48550/arXiv.2306.02465
  * Fazio et al. (2004) Fazio, G. G., Ashby, M. L. N., Barmby, P., et al. 2004, ApJS, 154, 39, doi: 10.1086/422585
  * Gáspár et al. (2021) Gáspár, A., Rieke, G. H., Guillard, P., et al. 2021, PASP, 133, 014504, doi: 10.1088/1538-3873/abcd04
  * Giavalisco et al. (2004) Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93, doi: 10.1086/379232
  * Gilmore et al. (2012) Gilmore, R. C., Somerville, R. S., Primack, J. R., & Domínguez, A. 2012, MNRAS, 422, 3189, doi: 10.1111/j.1365-2966.2012.20841.x
  * Gispert et al. (2000) Gispert, R., Lagache, G., & Puget, J. L. 2000, A&A, 360, 1, doi: 10.48550/arXiv.astro-ph/0005554
  * Granato et al. (2000) Granato, G. L., Lacey, C. G., Silva, L., et al. 2000, ApJ, 542, 710, doi: 10.1086/317032
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35, doi: 10.1088/0067-0049/197/2/35
  * Hauser et al. (1998) Hauser, M. G., Arendt, R. G., Kelsall, T., et al. 1998, ApJ, 508, 25, doi: 10.1086/306379
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Kirkpatrick et al. (2023) Kirkpatrick, A., Yang, G., Bail, A. L., et al. 2023, CEERS Key Paper VII: JWST/MIRI Reveals a Faint Population of Galaxies at Cosmic Noon Unseen by Spitzer, arXiv, doi: 10.48550/arXiv.2308.09750
  * Lacy et al. (2021) Lacy, M., Surace, J. A., Farrah, D., et al. 2021, MNRAS, 501, 892, doi: 10.1093/mnras/staa3714
  * Lagache et al. (2005) Lagache, G., Puget, J.-L., & Dole, H. 2005, ARA&A, 43, 727, doi: 10.1146/annurev.astro.43.072103.150606
  * Lagos et al. (2019) Lagos, C. d. P., Robotham, A. S. G., Trayford, J. W., et al. 2019, MNRAS, 489, 4196, doi: 10.1093/mnras/stz2427
  * Ling et al. (2022) Ling, C.-T., Kim, S. J., Wu, C. K. W., et al. 2022, MNRAS, 517, 853, doi: 10.1093/mnras/stac2716
  * Magdis et al. (2012) Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6, doi: 10.1088/0004-637X/760/1/6
  * Marleau et al. (2004) Marleau, F. R., Fadda, D., Storrie-Lombardi, L. J., et al. 2004, ApJS, 154, 66, doi: 10.1086/422584
  * Moster et al. (2011) Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113, doi: 10.1088/0004-637X/731/2/113
  * Murata et al. (2014) Murata, K., Pearson, C. P., Goto, T., et al. 2014, MNRAS, 444, 2346, doi: 10.1093/mnras/stu1617
  * Oliver et al. (1997) Oliver, S. J., Goldschmidt, P., Franceschini, A., et al. 1997, MNRAS, 289, 471, doi: 10.1093/mnras/289.2.471
  * Papovich et al. (2004) Papovich, C., Dole, H., Egami, E., et al. 2004, ApJS, 154, 70, doi: 10.1086/422880
  * Papovich et al. (2016) Papovich, C., Shipley, H. V., Mehrtens, N., et al. 2016, ApJS, 224, 28, doi: 10.3847/0067-0049/224/2/28
  * Pearson et al. (2010) Pearson, C. P., Oyabu, S., Wada, T., et al. 2010, A&A, 514, A8, doi: 10.1051/0004-6361/200913382
  * Pearson et al. (2014) Pearson, C. P., Serjeant, S., Oyabu, S., et al. 2014, MNRAS, 444, 846, doi: 10.1093/mnras/stu1472
  * Pérez-González et al. (2024) Pérez-González, P. G., Barro, G., Rieke, G. H., et al. 2024, arXiv e-prints, arXiv:2401.08782, doi: 10.48550/arXiv.2401.08782
  * Perrin et al. (2014) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 91433X, doi: 10.1117/12.2056689
  * Rieke et al. (2008) Rieke, G. H., Blaylock, M., Decin, L., et al. 2008, AJ, 135, 2245, doi: 10.1088/0004-6256/135/6/2245
  * Rieke et al. (2023) Rieke, M. J., Robertson, B., Tacchella, S., et al. 2023, ApJS, 269, 16, doi: 10.3847/1538-4365/acf44d
  * Robertson (2010) Robertson, B. E. 2010, ApJ, 716, L229, doi: 10.1088/2041-8205/716/2/L229
  * Rocca-Volmerange et al. (2007) Rocca-Volmerange, B., de Lapparent, V., Seymour, N., & Fioc, M. 2007, A&A, 475, 801, doi: 10.1051/0004-6361:20065217
  * Rowan-Robinson (2001) Rowan-Robinson, M. 2001, New A Rev., 45, 631, doi: 10.1016/S1387-6473(01)00142-7
  * Rowan-Robinson (2009) —. 2009, MNRAS, 394, 117, doi: 10.1111/j.1365-2966.2008.14339.x
  * Rowan-Robinson (2023) —. 2023, arXiv e-prints, arXiv:2312.03007, doi: 10.48550/arXiv.2312.03007
  * Rowan-Robinson et al. (2004) Rowan-Robinson, M., Lari, C., Perez-Fournon, I., et al. 2004, MNRAS, 351, 1290, doi: 10.1111/j.1365-2966.2004.07868.x
  * Rowan-Robinson et al. (2005) Rowan-Robinson, M., Babbedge, T., Surace, J., et al. 2005, AJ, 129, 1183, doi: 10.1086/428001
  * Schreiber et al. (2018) Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30, doi: 10.1051/0004-6361/201731506
  * Silva et al. (1998) Silva, L., Granato, G. L., Bressan, A., & Danese, L. 1998, ApJ, 509, 103, doi: 10.1086/306476
  * Spoon et al. (2004) Spoon, H. W. W., Moorwood, A. F. M., Lutz, D., et al. 2004, A&A, 414, 873, doi: 10.1051/0004-6361:20031656
  * Teplitz et al. (2011) Teplitz, H. I., Chary, R., Elbaz, D., et al. 2011, AJ, 141, 1, doi: 10.1088/0004-6256/141/1/1
  * Treister et al. (2006) Treister, E., Urry, C. M., Van Duyne, J., et al. 2006, ApJ, 640, 603, doi: 10.1086/500237
  * Trenti & Stiavelli (2008) Trenti, M., & Stiavelli, M. 2008, ApJ, 676, 767, doi: 10.1086/528674
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Viero et al. (2013) Viero, M. P., Moncelsi, L., Quadri, R. F., et al. 2013, ApJ, 779, 32, doi: 10.1088/0004-637X/779/1/32
  * Whitaker et al. (2011) Whitaker, K. E., Labbé, I., van Dokkum, P. G., et al. 2011, ApJ, 735, 86, doi: 10.1088/0004-637X/735/2/86
  * Wu et al. (2023) Wu, C. K. W., Ling, C.-T., Goto, T., et al. 2023, MNRAS, 523, 5187, doi: 10.1093/mnras/stad1769
  * Yang et al. (2023) Yang, G., Papovich, C., Bagley, M. B., et al. 2023, ApJ, 956, L12, doi: 10.3847/2041-8213/acfaa0