AGN Selection and Demographics: A New Age with JWST/MIRI

Understanding the co-evolution of supermassive black holes (SMBHs) and their host systems requires a comprehensive census of active galactic nuclei (AGN) behavior across a wide range of redshift, luminosity, obscuration level and galaxy properties. We report significant progress with JWST towards this goal from the Systematic Mid-infrared Instrument Legacy Extragalactic Survey (SMILES). Based on comprehensive SED analysis of 3273 MIRI-detected sources, we identify 217 AGN candidates over a survey area of $\sim$34 arcmin$^2$, including a primary sample of 111 AGNs in normal massive galaxies ($M_{*}>10^{9.5}~M_\odot$) at $z\sim$0--4, an extended sample of 86 AGN {\it candidates} in low-mass galaxies ($M_{*}<10^{9.5}~M_\odot$) and a high-$z$ sample of 20 AGN {\it candidates} at $z\sim$4--8.4. Notably, about 80\% of our MIRI-selected AGN candidates are new discoveries despite the extensive pre-JWST AGN searches. Even among the massive galaxies where the previous AGN search is believed to be thorough, 34\% of the MIRI AGN identifications are new, highlighting the impact of obscuration on previous selections. By combining our results with the efforts at other wavelengths, we build the most complete AGN sample to date and examine the relative performance of different selection techniques. We find the obscured AGN fraction increases from $L_{\rm AGN, bol}\sim10^{10}~L_\odot$ to $10^{11}~L_\odot$ and then drops towards higher luminosity. Additionally, the obscured AGN fraction gradually increases from $z\sim0$ to $z\sim4$ with most high-$z$ AGNs obscured. We discuss how AGN obscuration, intrinsic SED variations, galaxy contamination, survey depth and selection techniques complicate the construction of a complete AGN sample.


INTRODUCTION
Supermassive black holes (SMBHs) are found at the centers of most, if not all, galaxies and are thought to be intimately linked to the evolution of galaxies across cosmic time (see reviews by e.g., Alexander & Hickox 2012;Kormendy & Ho 2013).The formation and growth of SMBHs are regulated by the supply of accreting material with the efficiency influenced by various host galaxy properties at a wide range of physical scales (e.g., Hopkins et al. 2008;Hopkins & Quataert 2010;Inayoshi et al. 2020).As SMBHs accrete matter, they release enormous amounts of energy in the form of radiation and powerful outflows that shape the host systems through feedback mechanisms (e.g., Fabian 2012;Heckman & Best 2014).How to unravel the nature of this SMBH-galaxy connection has been a major quest in extragalactic astronomy.The growing phase of SMBHs, commonly known as Active Galactic Nuclei (AGNs), is partic-ularly interesting, as they represent the most active phase of the SMBH-galaxy interaction and have distinctive observational characteristics to separate them out from other objects.
Since the original identifications through optical spectroscopy by Seyfert (1943), extensive ground-based optical surveys have successfully yielded statistical samples of AGN over wide ranges of galaxy properties and cosmic time (e.g., Adams 1977;Filippenko & Sargent 1985;Ho et al. 1997;Richards et al. 2002;Kauffmann et al. 2003;Hao et al. 2005) , particularly with the development of various emission line diagnostics (e.g., Baldwin et al. 1981;Kewley et al. 2001;Kauffmann et al. 2003).Meanwhile, a large number of AGNs missed in the optical have been revealed by new techniques developed at other wavelengths such as radio (e.g., Edge et al. 1959;Heckman & Best 2014;Padovani 2016;Tadhunter 2016), X-rays (e.g., Gursky et al. 1971;Hickox & Alexander 2018) and ultraviolet excess (e.g., Markarian 1977;Schmidt & Green 1983;Hainline et al. 2011).Nevertheless, all these methods have limitations -only a minority of AGNs are prominent in the radio and the other techniques are significantly hampered by absorption by the dust and gas surrounding the AGN.Luckily, most of the luminosity of obscured AGNs emerges in the mid-IR, producing SEDs dominated by emission by warm dust.To mitigate the deficiencies of AGN selections at other wavelengths, various AGN infrared selection techniques, particularly mid-IR broad-band color diagnostics with Spitzer and WISE data have been developed (e.g., Lacy et al. 2004;Stern et al. 2005;Alonso-Herrero et al. 2006;Stern et al. 2012;Assef et al. 2013Assef et al. , 2018;;Mateos et al. 2012;Padovani et al. 2017;Hickox & Alexander 2018;Hviding et al. 2018), however, these methods are mostly effective on power law SEDs (Donley et al. 2012), i.e., on lightly obscured sources, and they leave the heavily obscured population only partially explored (see reviews by Hickox & Alexander 2018;Lyu & Rieke 2022a).In addition, the high-z obscured and faint AGN population (e.g., z ≳ 2) is even less explored given the very limited wavelength coverage and sensitivity beyond ∼8 µm.Consequently, there is wide disagreement in estimates of the fraction of the AGN population that has been overlooked due to strong obscuration, from 10% (Mendez et al. 2013) to 30% (Del Moro et al. 2016) to much more (Ananna et al. 2019), and how it evolves with AGN and/or galaxy properties.As pointed by the Astro2020 Decadal Survey (National Academies of Sciences 2021, Appendix D-Q3c.):"the cosmic census of AGNs is currently patchy, which limits our understanding of the co-evolution of galaxies and their SMBHs.The highest redshifts remain largely unprobed, our knowledge of the most heavily obscured AGNs is incomplete, even at the lowest redshifts, and nuclear activity in the lowest-mass galaxies is poorly constrained." With the successful launch and operation of JWST, we have entered a new era.Notably, from NIRSpec and NIR-Cam spectral observations at λ < 5µm, people have reported a number of high-z AGNs or AGN candidates through restframe UV-optical emission line properties (e.g., Harikane et al. 2023;Kocevski et al. 2023;Larson et al. 2023;Übler et al. 2023) with the highest records at z ∼ 10 ( Maiolino et al. 2023;Goulding et al. 2023).Since most of these searches target preferentially the broad-line systems that do not represent the whole AGN population, it is unclear how many AGNs have been missed.Previous systematic searches of obscured AGNs in the IR have been handicapped because of the very limited number of spectral bands with relevant data past ∼ 8 µm, leaving the SED models under-constrained.This obstacle has been lifted by the Mid-Infrared Instrument (MIRI; Wright et al. 2023) on JWST, which provides nine photometric bands from 5 to 26 µm that continuously sample the AGN hot dust emission features (λ peak ≥ 3µm) up to redshift z ∼ 8. Compared to the Spitzer mission, the MIRI sensivities are about 10-100 times higher with about eight times better spatial resolution at similar wavelengths, allowing the search and characterization of the AGN population in unprecedented precision, depth and redshift range (e.g., Yang et al. 2023) In this paper, we present the first results of AGN selection and demographics using the Systematic Mid-infrared Instrument Legacy Extragalactic Survey (SMILES; Rieke et al. 2017) led by the US MIRI Science Team.This survey is composed of 3×5 MIRI pointings in 8 imaging bands at λ ∼5-27 µm over the central region of the Great Observatories Origins Deep Survey-South (GOODS-S; Giavalisco et al. 2004), including the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006) (see Figure 1 for the survey footprint and sensitivities).In the pre-JWST era, this field was extensively observed across the electromagnetic spectrum, often with the deepest observations anywhere in the whole sky, offering the best ancillary dataset for AGN study.In the UV to the near-IR, it has been visited multiple times by the Hubble Space Telescope (HST) as well as many ground-based facilities that include both imaging and spectral observations.AGNs have been identified through emission line features (e.g., Santini et al. 2009;Silverman et al. 2010) or time variability (e.g., Pouliasis et al. 2019).In the X-ray, this field is covered as part of the Chandra Deep Field South (CDF-S; Giacconi et al. 2002) with a total exposure time about 7 Ms at 0.5-7 keV, enabling the identification of many X-ray detected AGNs (Luo et al. 2017).In the radio-band, super deep data at 3 GHz and 6 GHz with the Jansky Very Large Array (JVLA) have also been obtained and used in AGN searches (Alberts et al. 2020).In the infrared, very deep IRAC and confusionlimited MIPS 24 µm images from Spitzer are available and have been utilized to identify AGNs by mid-IR colors (e.g., Donley et al. 2008) and through SED fitting (Alberts et al. 2020;Lyu et al. 2022a).These multi-wavelength datasets and the corresponding results provide the best benchmark to show the advances brought by JWST and to characterize the relative performance of the MIRI AGN selection.
Besides the pre-JWST multi-wavelength dataset, our MIRI survey also has almost full overlap with the deep NIRCam imaging observations from the JWST Advanced Deep Extragalactic Survey (JADES; Eisenstein et al. 2023) and a large overlap with spectroscopic data from the First Reionization Epoch Spectrosopic COmplete Survey (FRESCO; Oesch et al. 2023), which can help to reveal the nature of these MIRI sources and characterize their properties.Compared to other MIRI imaging surveys in JWST cycles 1 and 2 (e.g., CEERS, COSMOS-Web), SMILES has larger sky coverage and/or more complete wavelength coverage.Combined with very deep JWST/NIRCam and HST data, this MIRI dataset provides a unique opportunity for AGN study over a wide range of redshifts and properties.
Taking advantage of the deep JWST (and HST) photometric data that continuously covers from 0.4 to 26 µm, we will conduct comprehensive SED analyses to identify AGNs with a wide range of luminosity and obscuration levels from z ∼ 0 to z ∼ 8.We will demonstrate the new discovery space offered by MIRI and present a major leap in AGN selection compared to all the previous IR missions.Combined with the super deep X-ray and radio data in the same field, we are able to build the most complete AGN sample in terms of bolometric luminosity and discuss the limits of different wavelengths and selection techniques.These results can establish the foundation of a number of follow-up investigations on AGN-related science such as AGN obscuration, luminosity functions, AGN duty cycles, host galaxy properties and AGN-galaxy relations.
This paper is structured as follows.Section 2 describes the observations, data reduction and photometry measurements.In Section 3, we introduce the AGN selection methods and present the results, which include SED identifications and other techniques with and without using the JWST data.In Section 4, we report the AGN number densities from MIRI, compare them to the AGN selected at other wavelengths and characterize the relative performance of different selection techniques across the electromagnetic spectrum.Section 5 characterizes the obscured AGN fraction from SED analysis and studies its possible evolution.The nature of AGNs missed at different wavelengths is also explored.In Section 6, we discuss various issues that can complicate the construction of a complete AGN sample.A summary is provided in Section 7.

JWST/MIRI Observations and Image Reduction
Our JWST MIRI imaging observations (JWST ID 1207; PI: George Rieke) comprise 15 separate pointings that produce a 3×5 mosaic at the central 34 square arcminutes of the GOODS-S/HUDF sky area.To improve the PSF sampling, and mitigate cosmic rays and detector artifacts, we used the 4-Point-Sets Dither optimized for point sources with the SMALL pattern.The full MIRI field was used with the FASTR1 readout pattern.For each pointing, the total science exposure was about 2.17 hours with the time spread over 8 MIRI filters: F560W, F770W, F1000W, F1280W, F1500W, F1800W, F2100W and F2550W with the shortest integrations lasting 10.7 min (F1000W) and the longest 36.4 min (F2100W).The exposure time per band was aimed at optimizing detection of a typical AGN or star-forming galaxy at z ∼ 2, except for F2550W which has a shorter exposure time aimed at z ∼ 1 galaxies.Most of these observations were successfully carried out during the scheduled Dec 7-14, 2022 time window with the same position angle (PA), but 2 pointings failed due to telescope guiding failures and the relevant data was re-taken on Jan 1 and 28, 2023.The different observing visits introduce position angle rotations of 26.6 and 48.2 degrees for these two pointings relative to the others and cause some irregularity in the final MIRI mosaic layout.
The data reduction used the JWST Calibration Pipeline v1.10.0 (Bushouse et al. 2022) and JWST Calibration Reference System (CRDS) context jwst 1084.pmap.Starting from the raw ("uncal") data files obtained from MAST 1 , the pipeline module calwebb detector1 was used with default parameters to apply detector-level corrections, including linearity, dark current, first and last frame effects, and jump detection.For the latter, we implemented newly introduced algorithms with default parameters to correct for cosmic ray showers, i.e., large cosmic ray events that can have an impact over a large number of pixels.
Next, the default step calwebb image2 was run to apply instrumental corrections and calibrations to produce "cal" files.Given that MAST data products can still suffer from strong background noise structures, at this point we apply a custom, external background subtraction from the Rainbow Database (Pérez-González et al. 2005, 2008;Barro et al. 2011a,b) JWST pipeline.Briefly, as explained in more detail in Álvarez-Márquez et al. (2023), we use an iterative process that progressively masks sources and median filters out large gradients and striping along detector columns and rows to produce "clean" cal files.For each (unfiltered) cal file, a "super background" is then constructed by median combining and scaling all other clean cal files.We found that better re- .Survey layout from X-ray to radio (left) and the 5σ flux limits of photometric bands at 0.3 -26 µm (right) used for SED fittings in GOODS-S/HUDF.On the left, the footprint of our MIRI survey is shown as the red-shaded region with JADES NIRCam deep (solid green line) and medium (dashed green line) observations, FRESCO NIRCam/grism coverage (solid blue line), HST ACS GOODS-S coverage (thick gray line), Chandra X-ray coverage (light blue shaded region) and JVLA radio observations at 3 GHz (dark green lines) and 6 GHz (dark orange lines).On the right, we denote the pre-JWST filters in gray, JWST/NIRCam filters in green and JWST/MIRI filters in red.Besides the flux limits of these filters as a function of wavelength, we also show some representative SEDs of obscured AGNs (orange), unobscured AGNs (blue) and starburst galaxies (magenta) at z=2 and z=4, where JWST observations are expected to bring major advances to identify and characterize these objects.
sults are achieved by using all (relatively contemporaneous) cal files to construct a super background as opposed to constructing a background per visit.The latter can result in oversubtraction around extended sources due to inadequate availability of background pixels at the source location in detector space.The final super backgrounds were then subtracted from each cal file and an additional 265x265 pix2 box median subtraction is applied to remove any remaining varying background, mainly caused by cosmic ray showers.
The background-subtracted cal files were then astrometrically-corrected on a per visit and per filter basis using tweakreg 2 with custom routines external to the JWST pipeline.The F560W image is first registered to the JADES image mosaic (which is registered to GAIA DR3) by creating matched catalogs of high S/N F444W and F560W sources, which will have similar morphologies and therefore centroids.The corrected F560W positions were then used to correct the next longer wavelength image and so on, so as to minimize the effects of morphological changes3 .The astrometric accuracy achieved by this process is 0.01-0.02′′ (1σ) for all filters except F2550W, which has a low number of high S/N sources and an astrometric accuracy of ∼ 0.04 ′′ .Finally, the astrometry-corrected and background-subtracted cal files were processed through calwebb image3, which produced the final mosaics with a pixel scale of 0.06 ′′ .The final 5σ point source sensitivities are 0.21, 0.21, 0.44, 0.60, 0.68, 1.7, 2.8, and 15 µJy in F560W, F770W, F1000W, F1280W, F1500W, F1800W, F2100W, and F2550W respectively, for apertures encompassing 65% of the encircled energy of the PSFs (see below).Compared to the JWST ETC predictions, our measured detection limits have been improved by a factor of 1.6-2 for F770W-F1800W.
Figure 2 shows a representative cutout of our JWST/MIRI images together with JWST/NIRCam, Chandra X-ray, and JVLA radio data in the same region.Almost all of the Xray and radio sources have been detected in the well-resolved MIRI images (see Section 2.5).The MIRI image resolutions are at 0.2-0.8′′ (PSF FWHM) while Chandra is at ∼0.5 ′′ (field center), JVLA at ∼0.3-1 ′′ and JWST/NIRCam at 0.02-0.16′′ , allowing reliable associations and accurate photometry across the full wavelength range.

MIRI Source Identification and Photometry Measurements
Source detection and photometric catalogs were produced using a modified version of the JADES photometric pipeline (M.Rieke et al. 2023).The object detection uses stacked MIRI F560W and F770W high S/N images and a series of iterative steps, starting with the creation of a low S/N threshold blended segmentation map which is then processed to remove spurious noise detections and deblend sources (see  M. Rieke et al. 2023 for details).The final segmentation map is used to define source centroids and 2.5x scaled Kron apertures, which are then used to measure photometry in all filters.Aperture corrections are applied by determining the fraction of flux outside a given Kron aperture using a model PSF for each filter.Model PSFs are created using WebbPSF (Perrin et al. 2014) 4 for the F1000W-F2550W filters.The F560W and F770W PSFs include an extended cross-like imaging artifact known as a "cruciform" (Gáspár et al. 2021) which is not yet adequately modeled by WebbPSF.As such, we constructed empirical PSFs using high dynamic range imaging of stars taken during commissioning (Gaspar, private communication).Photometric uncertainties are determined by placing random apertures across the mosaics to account for correlated pixel noise (e.g., Whitaker et al. 2011;M. Rieke et al. 2023).
The nature of our source detection does not impose a hard S/N threshold and our resulting catalog contains both real and spurious detections down to faint limits.We conduct visual inspections with the aid of deep NIRCam and HST images to throw out spurious sources.Finally, we have 3273 MIRI sources with 3-σ detection in at least one band.
By comparing with previous Spitzer photometry of the same sources, we found that our flux measurements in the F560W and F770W bands are systematically higher and concluded the origin is likely due to an underestimation of the cruciform during the flux calibration step.To mitigate the problem, we derived the correction factors as follows.We fitted the ASTRODEEP photometry up to IRAC band-25 of 25 isolated stars with single stellar population models from the Flexible Stellar Population Synthesis (FSPS) code6 and computed the synthetic MIRI fluxes in F560W and F770W.The median offsets between the observed (aperture-corrected) Kron flux and the model flux were then computed; the values are 1.26 for F560W and 1.04 for F770W.Finally, we divide the F560W and F770W fluxes by these factors during the SED fitting in Section 3. The photometry is archived in the v0.4.2 MIRI photometry catalog.

JWST/NIRCam, HST/ACS, HST/WFC3 Photometry and SED Constructions
Most of the SMILES survey area overlaps with the NIR-Cam image footprint from JADES (Eisenstein et al. 2023).This includes the JADES GOODS-S NIRCam/Deep (ID: 1180; PI: Daniel Eisenstein), partly taken from September 29 to October 10, 2022 and already publicly released (M.Rieke et al. 2023), and the JADES GOODS-S NIRCam/Medium (ID 1210; PI: Nora Luetzgendorf), taken from October 20 to 24, 2023.All these datasets include deep F090W, F115W, F150W, F200W, F277W, F335M, F356W, F410M, F444W images, offering critical insights into the nature and properties of the MIRI sources.In addition, two additional NIR-CAM pointings in the HUDF were conducted with the JWST Extragalactic Medium-band Survey (JEMS; Williams et al. 2023), offering F182M, F210M, F430M, F460M, F480M coverage over 15.6 arcmin 2 .In contrast to the previous multiband photometric analysis with data collected by different missions/instruments across the years, these NIRCam observations were carried out no more than 3 months earlier than the MIRI data, which together with the time dilation (typical redshifts are z=1-3), greatly mitigates the possible influence from AGN variability on the source SEDs. 7 Although the previous HST data is not as deep as JADES NIRCam data at the same wavelength, it is still indispensable as HST provides measurements at λ < 0.8 µm that are missed by NIRCam.The addition of HST bands also improves the SED analysis such as constraining the photoz.The HST/ACS F435W, F606W, F775W, F814W, F850LP and HST/WFC3 F105W, F125W, F140W, F160W images of this field have been reprocessed by the JADES team, including registering the astrometry against JWST/NIRCam images and measuring the HST and NIRCam fluxes in a consistent way.
For the NIRCam and HST photometry in the JADES field, we have adopted the PSF-convolved KRON flux from the JADES v0.8.1 catalog.We searched for the NIRCam counterparts for all the MIRI sources within a radius of 0.3 arcsec.For MIRI sources outside the current JADES footprint, we searched for the nearest CANDELS source within a radius of 0.6 arcsec and adopted the corresponding ASTRODEEP (Merlin et al. 2021) photometry up to IRAC band 2 (4.5 µm).
We have verified the quality of the JADES+MIRI photometry measurements with synthetic JWST photometry based on the SED fittings of the same sources with the AS-TRODEEP+Spitzer IRS 16 µm and MIPS 24 µm data.On a statistical level, they agree with each other within 1-4%.

Redshift Constraints for MIRI Sources
7 For a typical AGN with a BH mass ∼ 10 7 M ⊙ , the light crossing time for the UV/optical disk is 0.06-6 days; the observed relative continuum time lags in the rest-frame UV/optical band are a few days (see review by e.g., Cackett et al. 2021).Meanwhile, the rest-frame near-IR emission of the AGN torus is typically time lagged to the optical emission by τ ∼ 100 × (L ⊙ /10 11 L ⊙ ) 0.5 days with the variability signals quickly going away at longer wavelengths (see review by Lyu & Rieke 2022b).Given the fact that most AGNs are obscured and the time lag is further diluted by redshift with a factor of (1 + z), AGN variability should not impact the JADES+SMILES photometry for SED analysis in general.
To aid the SED analysis, redshifts of the MIRI sources are collected from various sources that can be grouped into three categories: • spectroscopic redshifts: These include high confidence spectroscopic redshift measurements published in the literature over last twenty years (Kodra et al. 2023;H. Hathi, private communication), the redshift catalog of the MUSE HUDF surveys published this year (Bacon et al. 2023), and also the recently released redshifts from JADES NIRSpec data release 1 (Bunker et al. 2023).
• grism redshifts: Most of the MIRI field has also been covered by FRESCO grism spectral observations with NIRCam/F444W (3.8-5.0 µm) at R∼1600 (Oesch et al. 2023).Since many objects only show a single emission line, we need to guide the redshift measurements with photo-z from JADES NIRCam as priors.Typically, the confidence level of the prism redshift is the same as the spectroscopic ones above given the spectral resolution and depth of the NIRCam grism observations.When a FRESCO grism redshift is not available, we search for the 3D-HST grism redshift from WFSC3/G141 (1.1-1.65 µm) at R∼130.Given the shorter wavelength coverage and poorer resolution, the latter is mostly only useful for objects at low redshifts.
• photometric redshifts: Thanks to the extensive multiband HST and JWST/NIRCam image coverage of the JADES, reasonably accurate photo-z measurements are possible for objects without spectroscopic observations.The JADES team has updated EAZY (Brammer et al. 2008) configurations and provided the photo-z measurements based on fitting HST and JWST/NIRCam photometry (see details in Hainline et al. 2023).We adopt this product by default.For objects outside the JADES/NIRCam coverage, we adopt the optimized photometric redshifts for CANDELS published in Kodra et al. (2023).
Figure 3 gives a summary of the redshift distribution of the MIRI sources.Despite the rather short exposure times of these MIRI bands, we detected over 30 objects at z > 6.0.About 45% of the sample have spectroscopic or grism redshifts.The EAZY-based photometric redshifts based on NIR-Cam data alone have 5% outliers and a scatter (excluding outliers) of ⟨z spec − z phot ⟩ of 0.024, 1σ (M.Rieke et al. 2023).
With the combination of MIRI photometry at longer wavelength, we expect the updated photo-z computed from our SED fitting with Prospector (see Section 3.1) to be improved further with even less outliers.Following the strategies in Lyu et al. (2022a), we crossmatched our MIRI sources with the CDF-S and JVLA catalogs to look for AGN signatures within a search radius of 3.0 ′′ .As the astrometry of the Chandra data was registered in the pre-Gaia era, we investigated the coordinate differences between the closest-matched MIRI and X-ray sources and applied a systematic shift of ∆RA=0.128′′ and ∆Dec=−0.288′′ to all the CDF-S sources to remove the relative astrometric offsets.There are 202 X-ray sources within the MIRI footprint but four of them do not have MIRI counterparts.Upon visual inspection, their X-ray detections are below 3-σ.For the radio sources, there are 183 objects but one does not have a MIRI counterpart.The nature of this radio source will be explored in the future.
In Section 3.4, we will discuss AGN selections for sources detected in X-ray or radio.

AGN IDENTIFICATIONS
With the data described above, we conducted SED analysis to search for AGNs.Section 3.1 introduces the semiempirical SED models used for AGN identifications.In Section 3.2, we conduct the initial round of SED analysis to get an overview of the sample properties of all the MIRI sources and check the performance of our models.Section 3.3 describes how the AGNs are identified from SED analysis with the removal of possible false positives.In Section 3.4, we summarize the AGN samples selected by other methods in the same field, including those utilizing the JWST data and those in the pre-JWST era.

SED Models and Fitting Setup
As discussed in detail in Lyu & Rieke (2022a), reliable AGN identification through SED fitting requires careful min-imization of the free parameters of the models while achieving accurate fits to the data.To meet this requirement, we introduced a fitting package for AGN SED selection and analysis in Lyu et al. (2022a) that features a modified version of the Prospector code (Johnson et al. 2021) with the Flexible Stellar Population Synthesis (FSPS; Conroy et al. 2009;Conroy & Gunn 2010) model for the stellar component; semiempirical SED models for the AGN component (Lyu et al. 2017;Lyu & Rieke 2017, 2018;Lyu & Rieke 2022a); and a well-proven template for dust emission from the star-forming galaxies (Rieke et al. 2009;Lyu et al. 2022a).The same tool is also used for this work with a few updates described below.
In Figure 4, we show the SEDs of different components and illustrate some of their variations.
The configuration of the stellar component is very similar to Lyu et al. (2022a) -we assumed the same Kroupa initial mass function and delayed-tau star-formation history but adopted the Calzetti attenuation curve with a flexible slope as in Kriek & Conroy (2013).Since there are objects at very high redshifts, IGM absorption is also included.In total, there are seven free parameters to describe the stellar component: (1) stellar mass formed (mass); (2) stellar metallicity (logzsol); (3) stellar age (tage); (4) e-folding time of the star formation history (tau); ( 5) the attenuation level at 5500 Å (dust2); ( 6) the power law index of the attenuation (dust index); (7) the IGM absorption factor (igm f actor).For stellar populations older than ∼ 10 Myr, regardless of the metallicity and star formation history, the near-to-mid-IR infrared photospheric emission is dominated by red giants and supergiants with very similar spectra, so any degeneracy of the parameters (2)-(4) will not influence the AGN identification.
For the star forming galaxy (SFG) dust emission component, by default we also adopt the template for log(L IR /L ⊙ ) = 11.25, which is the best description of the average properties of normal SFGs at a wide range of redshifts that constitute our primary sample (see Appendix A in Lyu et al. 2022a for details).Furthermore, as shown later, our MIRI data is deep enough to probe the dwarf galaxy regime, where the systems can be low-metallicity and their dust emission SEDs likely behave differently from those of more luminous SFGs.These galaxies constitute our extended sample.We therefore adopt the Haro 11 SED template developed in Lyu et al. (2016) as the choice to match the dust emission SED for dwarf galaxies.This object has been selected from the Dwarf Galaxy Survey (Rémy-Ruyer et al. 2013) to represent the extreme SED behavior for low-metallicity starforming galaxies and it features a much higher dust temperature (T =46.5 K) than normal dust SFGs (T ∼20-30K) with a warmer mid-IR SED and weak PAH features.De Rossi et al. (2018) have also shown that the Haro 11 SED template is the best choice among other options for young, presumably sub-solar metallicity, galaxies undergoing extremely vigorous episodes of star formation.In Section 3.2, we will test the performance of these two sets of galaxy dust emission templates with MIRI data and explore the situation where the Haro 11 template is preferred for our SED fits.Once the template is decided, there is only one free parameter L IR obs to scale it to match the data together with other components.
Compared to Lyu et al. (2022a), several improvements have been made for the AGN component: • Based on the SDSS DR7 optical spectral archive (Abazajian et al. 2009) and near-IR AGN templates in the literature (Glikman et al. 2006;Hernán-Caballero et al. 2016), we increased the resolution of the AGN template from the UV to the near-IR and added AGN emission lines to the model.These emission lines are further decomposed into narrow and broad components with the separation at FWHM=1200 km/s and these two groups of lines can be modified in different ways (e.g., obscuration at different physical scales).
• We introduced a hybrid extinction configuration for the AGN component that considers the possibly different dust grain properties at different scales or angular directions.This setup features a Small Magellanic Cloud (SMC) curve for the commonly seen UV-optical extinction in type-1 AGNs and an empirical attenuation law constructed for the AGN IR obscuration.The former extinction law is associated with small dust grains at large physical scales far away from the central engine and the latter is for the circumnuclear obscuration, such as the dusty torus, which is likely to be dominated by larger dust grains because of the harsh environment (Baskin & Laor 2018).
Figure 4 presents this updated AGN model with various emission lines and shows how the SED can be modified by the different levels of extinction characterized by the two attenuation laws parameterized by tau opt and tau ir.In our fitting, the AGN continuum and broad emission line component are bound together and can be obscured by both attenuation laws (i.e., tau opt and tau ir) while the relative strength of the narrow line component can be only changed by the SMC extinction (tau ir).We adopt a uniform prior for tau opt = [0, 3.0] and tau ir = [0, 20.0] as these ranges cover the majority of AGN obscuration without too much redundancy.For the purpose of AGN selections, we also adopt the normal AGN template following  HST).For objects that only have JADES or 3D-HST photoz, we will re-compute the photo-z with the new SED fittings that include the MIRI data points.Accordingly, we have up to 12 (11 if the redshift is fixed) free parameters for the SED fitting.We use Dynamic Nested Sampling to find the bestfitting results and assess the model degeneracies.Readers are highly encouraged to read the relevant sections in Lyu et al. (2022a) for details.A summary of the fitting parameters is given in Table 1.
Our approach has a number of advantages.First, our models for the AGN component and the SFG dust component are derived from empirical observations and have been extensively tested, as elaborated in detail in Lyu et al. (2022a); Lyu & Rieke (2022a).This allows us to narrow the number and range of free parameters while retaining realistic fits, in contrast to approaches that are more strongly based on theoretical models.In addition, the SED model of our AGN com-ponent not only provides a good fit to the AGN observations, but can change continuously from the unobscured phase to the most heavily obscured phase, offering an unbiased view on the distribution of AGN obscuration properties.Moreover, compared to most other tools, our AGN model includes broad and narrow emission lines, a typically over-looked but now very needed feature for the much improved photometric data in the JWST era.Finally, although preliminary, our SED fittings have the option to use an empirical low-metallicity galaxy template to fit the IR SEDs of the dwarf galaxy population.
More details of this SED fitting package, including its public release, will be presented in a separate paper in the future.

SED Fitting of MIRI Sources
We fit the SEDs for all the MIRI sources with the redshift priors described in Section 2.4.Within the current JADES footprint, the fits can utilize a maximum of 27 photometric bands that include five HST/ACS bands at 0.44 -0.9 µm, fourteen JWST/NIRCam bands at 0.9 -4.4 µm and eight JWST/MIRI bands at 5.6 -25.5 µm.For the 89 MIRI sources outside the current JADES NIRCam catalog footprint, besides the MIRI photometry, we have 23 bands from the AS-TRODEEP catalog at 0.3 -4.5 µm, which includes HST ACS and WFC3, ground-based deep optical and near-IR observations and Spitzer/IRAC bands 1 and 2. In the first round, we adopted the normal luminous SFG dust template for the fitting to get an overall picture of the sample properties and to test the performance of our SED models.

Testing the SED Models for the Massive Galaxy Population
To demonstrate the robustness of our SED models, we present some representative SED fittings of dusty starforming galaxies (DSFGs) and quiescent galaxies (QGs) from z ∼ 0.2-4 in Figure 5.As we require good MIRI coverage to test our SED templates, most of these galaxies are relatively massive for their redshifts, by selection (i.e., M * ≳ 10 9.5−10 M ⊙ ) and they form our primary sample.For DSFGs, various SED features are successfully matched by our templates, showing that the MIRI bands can be used to trace the various PAH features at e.g., 3.3, 6.2, 8.6, 11.3 µm.For the QG population, the MIRI photometry is deep enough to constrain the stellar Rayleigh-Jeans tail at rest-frame 5 µm even for galaxies at z ∼ 4, as shown in the bottom-right panels of Figure 5.
The high-quality MIRI photometry as well as the satisfactory performance8 of our SED models for massive DSFGs and QGs establish a solid foundation for identifying AGN signatures through SED analysis in such systems.

Testing the SED Models for the AGN Population
Figure 6 shows some representative examples of MIRI sources with AGN evidence revealed by SED fittings.Compared to galaxies without AGNs shown in Figure 5, the SEDs of these objects present much larger variations that can be explained by the different level of AGN-galaxy contrast and the SED variations of the AGN component itself.Our examples include AGNs with a wide range of redshift, obscuration level and galaxy contamination as revealed by the very successful SED decompositions.

The Complicated Situation of the Dwarf Galaxy Population
Besides reproducing the photometric SEDs, our fittings provide measurements of the galaxy stellar mass.As shown in Figure 7, this MIRI sample includes a considerable number of dwarf galaxies (M * < 10 9.5 M ⊙ ), where our templates for normal SFGs may not work.To determine where this may become an issue, we visually inspected the fittings of all the MIRI sources and picked out those pure SFGs  where the updated Rieke et al. (2009) template gives a good match (orange crosses).These "normal" SFGs have a stellar mass distribution peaked between 109.5 and 10 10 M ⊙ , which is about 0.5 dex above the median stellar mass of the whole MIRI galaxy population.However, we found that the IR SED features (e.g., PAH strength) of many low-mass galaxies are not captured by these "normal" SFG templates.Compared to the normal galaxy populations, these low-mass galaxies could have different IR SED behavior associated with lowmetallicity ISMs (e.g., Rémy-Ruyer et al. 2013;Shivaei et al. 2020Shivaei et al. , 2022)).In particular, the reduced PAH emission strengths and warmer dust continuum emission associated with the dwarf galaxy component can mimic the mid-IR colors of obscured AGNs to some level (e.g., Hainline et al. 2016), making the AGN identification with classical SED models not convincing.Thus, we need to introduce a method to remove most, if not all, false positives of AGN candidates in this extended sample.
Based on empirical observations, a correlation between PAH strength and metallicity has been revealed from the lowz Universe (e.g., Engelbracht et al. 2005Engelbracht et al. , 2008;;Madden et al. 2006) to cosmic noon at z ∼ 2.0 (Shivaei et al. 2017).While the PAH intensity and relative strength of different bands stay almost the same at higher metallicity, below ∼ 12 + log(O/H) = 8.2 the PAH strength drops suddenly (Engelbracht et al. 2005(Engelbracht et al. , 2008;;Marble et al. 2010;Aniano et al. 2020). 9Given the correlation between the galaxy stellar mass and metallicity observed over a wide range of stellar mass and redshift (e.g., Tremonti et al. 2004;Lee et al. 2006;Zahid et al. 2013), we use the stellar mass from our SED fittings to approximate the gas metallicity 12 + log(O/H) following the empirical fittings in Ma et al. (2016), and adopt the dwarf galaxy dust template (see Section 3.1) when the metallicity (or stellar mass) threshold is not met.According to fittings in Ma et al. (2016), the corresponding mass threshold M * ,lim as a function of redshift z for 12 + log(O/H) ∼8.2 is log(M * ,lim /M ⊙ ) = 10.714+ 2.657 exp(−0.43z)(1) This curve is shown as the red solid line in Figure 7.To account for the uncertainties of e.g., stellar mass measurements and the dispersion of the mass-metallicity relation, we will compare the normal and dwarf galaxy dust template fittings for all objects within ±0.3 dex of M * ,lim (dashed red lines), and use the normal galaxy dust model results for galaxies above +0.3dex of this mass threshold and use the dwarf model below -0.3 dex of it.In addition, this threshold is not a hard cut as we will compare the fittings with the normal SFG and Haro 11 templates for objects within ±0.3 dex of the mass cut, corresponding to a 12+(O/H) range from 8.3 to 8.1.
Despite these steps, it is possible that the AGN selections in our primary sample are mildly contaminated by cases with stellar-heated hot dust.In the case of our extended sample, Haro 11 is one of the most extreme cases with significant warm dust emission and weak PAH strength among the Dwarf Galaxy Survey (Rémy-Ruyer et al. 2013).The usage of this template therefore may have biased our results against AGN identifications in this sample, as discussed in Section 3.3.2.

AGN Identifications through SED Analysis
We now describe our AGN identifications through SED analysis for the three sub-samples of our full object list.Following the same strategies in Lyu et al. (2022a), we picked out AGN candidates by analyzing the quality of the SED fittings, including inspections of any degenerate models from the posterior distributions of the fitted parameters.Typically, the AGN solution will be accepted if the distribution of the fitted AGN luminosity parameter L AGN has one single isolated peak and the AGN-dominated bands are not all noisy (i.e., S/N<3.0).In addition, the high spatial resolution provided by NIRCam and MIRI allows us to further check if the AGN evidence inferred from the SED analysis is supported by the galaxy having a point-like core.Depending on galaxy stellar mass and redshift, the AGN identifications have some additional considerations described below.

AGNs in Massive Galaxies (Primary Sample)
For typical massive galaxies (M * ≳ 10 9.5 M ⊙ ) at z =0-4, the MIRI data has good coverage of the AGN warm dust features and our templates are well-tested.As a result, we treat the selected objects as SED-identified AGNs.In total, we identify 111 AGNs under this category.
In the top six panels of Figure 6, we show some representative SED fittings together with some NIRCam and MIRI images at the source location.In general, once the AGN emission dominates in a spectral band according to our SED decomposition, the corresponding radial profile of the source in that band is consistent with the instrument PSF, supporting the validity of our SED interpretation.

AGNs in Dwarf Galaxies (Extended Sample)
As discussed in Section 3.2.3,low-metallicity dwarf galaxies have IR SED shapes that can mimic the relatively hot dust emission features of AGN (see one example in Figure 8).For such galaxies, we replace the massive "normal" SFG template with a dwarf galaxy template based upon Haro 11 to fit the SEDs and look for AGN candidates.An additional challenge is the reduction of MIRI coverage at longer wavelengths at the lower infrared luminosities of these objects.We highlight all AGNs identified with the normal SFG template regardless of the object redshift and stellar mass in light blue on the left, and the refined AGN sample after adopting the Haro 11 dwarf SFG template to reject possible false positives among low-mass galaxies in dark blue on the right.In both two panels, we show the normal SFG galaxies as orange crosses and other galaxies as gray dots.The red solid line represents the stellar mass threshold for 12+log(O/H)∼8.2 with the ±0.3 dex variation range as red dashed line.If we adopt the same normal galaxy template for all stellar masses (left panel), there would be 500 AGN candidates selected.In contrast, we end up with 217 AGN candidates after considering the dwarf galaxy complication (right panel), as described in Sections 3.2.3 and 3.3.2 These two effects restrict our identification of AGNs in these galaxies to cases with SEDs significantly bluer (hotter) than that of Haro 11.Our selection is probably missing some obscured AGNs with SEDs similar to that of Haro 11.However, using the Haro 11 template greatly reduces the chance of a source being mis-classified as an AGN, especially for SF-AGN composite systems.We have shifted the mass threshold to test whether the value of our adopted mass threshold is reasonable.If we shift it upward by 0.3 dex (corresponding to higher gas metallicity), some previous AGNs selected with the normal SFG template are re-classified as pure star-forming (dwarf) galaxies as expected.However, most such objects are either X-ray sources or show pointlike morphology in the NIRCam or MIRI images, indicating that they might be real AGNs and the additional shift of the threshold is not desirable. 10As a result, we think the current mass threshold in Equation 1 is valid.
We adopted the definition of dwarf galaxies to be M * ≲ 10 9.5 M ⊙ , regardless of the object redshift.In the first round, 187 AGN candidates that belong to this category are picked out.We have visually inspected all of these fits and rejected cases where the results are ambiguous due to low signal to noise, inconsistencies in the photometry from band to band, etc.Finally, 86 convincing AGN detections are left.This value is a lower limit given the limitations discussed above, particularly if AGNs with SEDs similar to that of Haro 11 have been overlooked.Methods to improve the AGN identifications throughout this extended sample need further development.
In Figure 6, the second row from the bottom shows two example dwarf galaxies where the AGN revealed by the MIRI SEDs is consistent with the obscured AGN template.The SED of MIRI 431 even shows evidence of strong silicate absorption, a typical IR tracer for a heavily embedded nucleus.

AGNs at z > 4 (High-z Sample)
For galaxies at z > 4, the MIRI data at the longest wavelengths are typically not deep enough to constrain the AGN hot dust SED shape well and the other MIRI bands are probing the rest-frame near infrared (or even optical) bands where some of the signatures of obscured AGNs are no longer accessible.Therefore, the identification of high-z AGNs is relatively less robust and likely only yields a subset.Given the possibility that the IR SEDs of galaxies at z ≳ 4 resemble that of Haro 11 (De Rossi et al. 2018), we adopt the Haro 11 template for the SED fittings of these high-z systems to make a conservative selection.In total, we identify 20 AGN candidates at z > 4 with the highest redshift at z =8.4.
Figure 9 presents six examples of AGNs identified at z > 4. All of these galaxies have spectroscopic redshifts and the AGN nature of MIRI 1104 is confirmed by the broad Hα emission revealed by the FRESCO spectrum (Matthee et al. 2023).The infrared excesses in all of these cases extend to wavelengths significantly shorter than the excesses of starforming galaxies.These examples demonstrate that MIRI observations can find obscured AGNs at z > 4, although the samples will be incomplete.* * * Combining the results above, we found 217 AGNs (or candidates) from SED analysis in total as listed in Table 2.In Figure 10, we present the AGN luminosity distribution as a function of redshift.Compared to the AGN sample selected in the pre-JWST era (Lyu et al. 2022a), the new sample contains more objects with lower bolometric luminosities and/or higher redshifts.
As mentioned above, there are several caveats in the AGN search among the extended and high-z samples and our selections are conservative to reduce the number of false positives.In the future, we will improve these selections with more comprehensive analysis and additional data.

Other Selections in the SMILES Footprint
Besides SED fitting of MIRI sources, there are many other methods to identify AGNs.Now we describe the identification of radio-loud AGNs, broad-line AGNs and mid-IR variable AGNs from JWST data with the combination of datasets at other wavelengths when necessary, and then summarize AGN samples identified without the usage of JWST data within the same footprint.

Radio-Loud AGNs with MIRI and JVLA
Since we now have deeper F2100W photometry than the previous Spitzer/MIPS 24µm data, it is possible to improve the previous radio-loud AGN selections based on the ratio between the radio and mid-IR fluxes (i.e., [rtype rl] in Lyu et al. 2022a).As the F2100W data is deeper than F2550W, we adopt q 21,obs = log(S 21µm,obs /S 1.4 GHz,obs ) to reveal outliers from the radio-infrared correlations for dusty starforming galaxies, where S 21µm,obs is the observed MIRI F2100W flux and S 1.4 GHz,obs is the observed radio flux at 1.4 GHz.For the latter, we assume a power-law index of α = −0.7 to estimate the 1.4 GHz flux density from 3 GHz, similar to Lyu et al. (2022a).Following the same strategy in Alberts et al. 2020, we determine the locus of SFGs by calculating q 21,obs at z =0-3 from the Rieke et al. ( 2009) SFG templates at log(L IR /L ⊙ ) = [11, 11.5, 12].An outlier from the SFG radio-infrared relation is defined as being 0.5 dex below the midpoint of SFG q 21,obs locus.Using this criterion, we identify 8 radio-loud AGNs among 167 z < 3.0 MIRI sources that have been detected at 3 GHz, as shown in Figure 11.In total, there are 39 AGNs in this radio-detected sample and the radio loud fraction is about 20%.

Spectroscopically-identified broad-line AGNs with JWST
FRESCO has conducted deep NIRCam/grism observations with the F444W filter for the 62 arcmin 2 area of GOODS-S that has a considerable overlap with our SMILES survey (see Figure 1), offering the chance to identify AGNs through broad-line emission.At z > 5, only one AGN has been reported at z=5.481 in GOODS-S, based on the broad Hα emission line (Matthee et al. 2023).Our MIRI SED fitting also identified this AGN (MIRI ID v0.4.2: 1104).We have also looked for AGN candidates at lower redshifts based on broad emission lines from Paα, Paβ and He II (Lyu et al, in prep).All the sources we found have been previously identified in Lyu et al. (2022a).
In addition, a few AGNs in GOODS-S have been identified from JADES/NIRSpec spectral observations (Bunker et al. 2023).In terms of broad-line AGNs, there are only two: 10013704 (z=5.92) and 8083 (z=4.64)(Maiolino et al. 2023).In contrast to the FRESCO AGN, these objects are much fainter by selection.In our MIRI images, they are only marginally detected in F560W and/or F770W without additional measurements in other bands and thus we cannot conduct robust SED fittings for AGN selection.

Mid-IR Variable AGNs with NIRCam and IRAC
Given the deep IRAC and NIRCam images available across the field, we can also look for long-term mid-IR variable AGNs.We have degraded the JADES/NIRCam F356W images to match the PSF of Spitzer/IRAC band-1 mosaics and carried out local difference photometry to look for variable objects (Lyu et al., 2023, in prep).Nevertheless, no new AGNs have been revealed.We color-code different sources on the main plot and show their relative distribution of q 21,obs in the right panel.

AGNs with weak Pa α
Potential identification of AGN candidates from weak Pa α compared with mid-infrared continuum is discussed in Appendix A. We have not yet applied this method; with current data it adds a small number of candidates.

AGN Selected in the Pre-JWST Era
In Lyu et al. (2022a), we have presented an AGN sample selected based on the pre-JWST data from X-ray to radio with eight different methods that included X-ray luminosity, X-ray-to-radio luminosity ratio, UV-to-mid-IR SED fitting, mid-IR IRAC color, radio-to-mid-IR ratio (for radioloud AGN), radio slope (for flat spectral radio sources), and time variability in the X-ray and optical bands.We can now replace the AGN sample selected by IRAC color and SED fittings from the Lyu et al. (2022a) pre-JWST catalog with the new SED-identified AGNs obtained in this work using the improved JWST photometry.In the MIRI survey footprint, there are 203 pre-JWST candidate AGNs in total; only five of them do not have a MIRI counterpart.Among these pre-JWST AGNs without MIRI source association, four are identified by variability and one is identified by the X-ray to radio flux [x2r].The latter object is also an X-ray source with very low S/N detection that does not have a MIRI counterpart, which therefore may be spurious.
We can also compare the SED fittings with the Spitzer and JWST data for the same sources.Two examples are given in Figure 12.The first object (MIRI 32) does not show evidence for an AGN, since the SED constraints for wavelengths longer than IRAC band 2 are very limited in the shallower Spitzer data.In contrast, the new MIRI data points reveal this object is a composite galaxy with an obscured nucleus in the mid-IR.The second object (MIRI 1691) seems to be an obscured AGN based on the limited wavelength coverage by the Spitzer data while the new fittings with MIRI data prefer the star forming galaxy solution.In fact, the AGN solution can be ruled out due the extended MIRI source morphology.
For the other selections based on X-ray or radio properties, the new mid-IR JWST data do not have any influence.

AGN STATISTICS AND COMPARISON OF SELECTION TECHNIQUES
With the various selection methods presented above, we can now compile all known AGNs in the MIRI footprint and compare the relative performance of different selections.As described above, four variable AGN candidates and one Xray-to-radio AGN candidate do not have MIRI counterparts and they are likely to be spurious; thus we will ignore them and focus the discussion on the MIRI-detected sources.

MIRI AGN Number Density
Given the survey area of ∼34 arcmin 2 , the total number density for the robust MIRI AGNs in the primary sample is 3.4 per arcmin 2 , or nearly eight AGNs per MIRI pointing.In addition, there are 2.9 AGNs per arcmin 2 from the extended sample.For each band, the total MIRI-selected AGN (candidates) number densities are 6.4, 6.3, 5.6, 5.0, 4.7, 3.7, 3.2 and 2.1 arcmin −2 from F560W to F2550W (3σ detections, corresponding to flux limits at 0.16, 0.13, 0.26, 0.36, 0.41, 1.0, 1.7 and 9.0 µJy).These estimates will evolve as we develop more sophisticated methods to differentiate the emission by warm dust in extended-sample galaxies from that of embedded AGNs.For now, these values should be treated as lower limits for total numbers of AGNs detected by MIRI since (1) due to the reduced sensitivities at longer wavelengths, we are likely missing obscured AGNs at high-z; and (2) for the AGNs in the extended-sample galaxy population, we are using an extreme SFG template based on Haro 11 and it is possible that some AGNs are missed due to the conservative selection; (3) there are AGNs that have been identified by other selection techniques and that are detected by MIRI but are missed by the SED analysis, as shown in the upcoming section.
In Figure 13, we plot the number density of the various MIRI sources as well as the total AGN fraction as a function of observed flux for all the MIRI bands.In every band, the AGNs in the primary sample of massive galaxies have dominated the distribution at the bright end while those in the extended sample of low-mass galaxies and high-z galaxies contribute mainly at the faint end.The total AGN fraction typically dominates the bright end with a fraction >40-60% and gradually drops toward the faint end to about 13-16%.For F560W, F770W and F1000W, the AGN fraction presents an apparent drop at the faint end, which is expected given our reduced selection efficiency since these sources are not detected at longer wavelengths.For F1280W, F1500W and F1800W, the AGN fraction seems to present a trend that decreases from the bright end to the intermediate flux followed by a gradual upturn towards the faint end.Given that the redshift distribution of our MIRI sources is peaked around z ∼1-1.5, such trends are likely caused by the strong 7.7µm PAH features from DSFGs.

Comparison to Previous AGN Samples
In Figure 14, we present the Venn diagram to show how the MIRI SED-identified AGN sample intersects with the CDF-S AGN catalog (Luo et al. 2017) and the pre-JWST AGN catalog (Lyu et al. 2022a).Over the same area, there are 94 X-ray detected AGNs reported from the 7Ms CDF-S catalog with an AGN number density about 2.8 per arcmin 2 .Among the 217 MIRI AGNs, only 48 of them (22%) have X-ray detections in the 7Ms CDF-S catalog and 38 of them (18%) are reported as AGN by Luo et al. (2017).In other words, our relatively shallow MIRI survey (∼2.17 hours per pointing) yields 2.4 times more AGNs than the deepest X-ray survey (∼ 1800 hr) and 80% are new discoveries.
Within the same footprint, there are 181 sources in the pre-JWST AGN catalog that include not only X-ray detected but also radio-and mid-IR detected AGNs.Still, about 78% of the MIRI AGN sample does not overlap with the pre-JWST AGN sample, demonstrating the huge discovery space offered by JWST.Among the 168 AGNs only identified by MIRI SEDs, there are 64 AGNs in massive galaxies at z <4.0, 85 AGN candidates in dwarf galaxies and 19 highz AGN candidates.The latter two populations are mostly found by MIRI, as there are only 1 dwarf AGN candidate and 1 high-z AGN candidate also identified by other selections.
Meanwhile, there are 56 CDF-S AGNs and 134 pre-JWST AGNs not identified by the MIRI SED analysis, which contribute to ∼60% and ∼74% of the corresponding catalog within the MIRI survey footprint.As pointed out in Lyu et al. (2022a), there is no single method or wavelength can identify all the AGNs.The results here provide very clear support to this argument -despite the much improved AGN hunting capability offered by MIRI, there are always some AGN population(s) left out and the combination of all possible selections across the electromagnetic spectrum is necessary to complete the AGN census.
To illustrate the new discovery space provided by MIRI in terms of AGN bolometric luminosity, Figure 15 shows the detection limits of the Chandra X-ray and JVLA radio bands as a function of redshift.In the X-ray, the Chandra data can only probe the most luminous MIRI AGNs, which is just the tip of the iceberg.This is consistent with the fact that only 33 out of 217 MIRI AGNs (15%) have reported X-ray luminosity that passes the AGN selection criteria.Most of the MIRI AGN sample are much fainter and the MIRI limit is about one order of magnitude deeper than the Chandra one.
Regarding the radio emission, there is a wide range of intrinsic variation of its strength from radio-loud to radio-quiet populations and we only show the AGN bolometric luminosity limits from the JVLA data from some average templates (see Lyu et al. 2022a for details about the template models).It is very clear that the current JVLA data is not deep enough to detect the radio emission from the radio-quiet MIRI AGNs and any detected radio emission is most likely coming from the host galaxies (Alberts et al. 2020).The JVLA data can probe the radio properties of bright radio-loud MIRI AGNs but would miss relatively faint radio-loud AGNs.To compare the yields of different selection methods and the possible overlaps, we used the UpSet visualization tool (Lex et al. 2014) to show the intersection of different AGN samples in Figure 16.Similar to the conclusions reached Lyu al. the top most efficient selection techniques are SED fittings, X-ray to radio luminosity cut ([X2R]) and X-ray luminosity cut, which identify 97.5% of the whole sample.For the remaining AGNs, 1.6% are identified by variability and 0.9% by the JWST-JVLA radio-loud AGN selection.

The Relative Performance of Different Selection Methods
Despite the substantial improvement of the AGN SED identification with JWST data, about 34% of the known AGNs in the field have been missed and most of them have been identified by [X2R] (23% of the whole AGN sample), illustrating the need to combine multi-wavelength selections to reach a more complete AGN census.(We will discuss the nature of these objects in Section 5.2.)Meanwhile, 53% of the AGN are only revealed by the JWST SED analysis.The other selection methods such as [xtype lumcut], [var], [jwst rl] all include unique AGNs that are solely picked out by one method, despite their much smaller contribution compared to the [jwst sed] sample.Given the discussion in Sec-tion 4.2, the yields of all the other selection methods are limited by the depth of X-ray or radio data and there should be additional galaxies detected by current MIRI data that have AGNs missed by the current SED analysis.
Regarding the different AGN groups (primary, extended and high-z), the vast majority of extended-sample AGN candidates (89%) and high-z AGN candidates (100%) have been identified through the MIRI SED analysis.For all the normal AGNs in massive galaxies at z < 4, the MIRI SED analysis has identified 114 out of 211 these objects (∼54%).Notably, after merging all AGN samples together, 18 out of 20 high-z AGN candidates (∼92%), 88 out of 109 dwarf AGN candidates (81%) and 72 out of 211 normal AGNs (34 ± 4%) are only identified by MIRI SED fittings, where the error is from the sample size.In our previous study (Lyu et al. 2022a), which was less effective at identifying AGNs due to less complete data but covered a five times larger area, we found that ∼ 20% of the AGNs were very strongly obscured and 26% were undetected in the deepest Chandra data, qualitatively in agreement with the new estimate of the fraction missed in previous studies.The numbers in the brackets are the total numbers of AGN from each source.For objects identified by MIRI, we also plot a horizontal stacked bar chart to show the number of normal AGNs (green), dwarf AGNs (red) and high-z AGNs (yellow) for each subset.

Comparison to MIRI AGN Studies in other fields
Pre-JWST efforts to use infrared photometry to identify AGNs are reviewed by Lyu & Rieke (2022a).Although they began with IRAS, they came into their own with the wide and deep surveys obtained with Spitzer and WISE (e.g., Lacy et al. 2004;Stern et al. 2005;Alonso-Herrero et al. 2006;Donley et al. 2008;Mateos et al. 2012;Asmus et al. 2020;Hviding et al. 2022).The possibilities are greatly expanded with JWST.
Based on four MIRI pointings from CEERS, Yang et al. (2023) reported 102 SF-AGN mixed systems and 25 AGNs from a total of 560 MIRI detections over a sky coverage of ∼ 9 arcmin 2 by fitting the galaxy SEDs with CIGALE.If we combine their SF-AGN mixture and AGN categories, their reported AGN number density is ∼14 arcmin −2 .Our fits include SF-AGN mixtures, so they can be compared directly.We find significantly lower numbers: 3.4 arcmin −2 for the main sample and 2.9 arcmin −2 for the extended one, for a total of ∼ 6.3 arcmin −2 .The high AGN numbers in Yang et al. (2023) presumably come with a caution because of the possibility of significant contamination by low-mass galaxies.Some other important distinctions between the two studies are that: (1) we cover a larger field, with five times as many overall MIRI detections and with more ultra-deep ancillary data; (2) we have more accurate redshifts, with ∼ 40% spectroscopic ones and somewhat more accurate photometric ones for the rest; and (3) we separate galaxies assumed to have PAH emission from those without on the basis of mass/metallicity, whereas they leave this difference as a free parameter. 11ased on the same CEERS dataset, Kirkpatrick et al. (2023) found ∼ 3 AGNs per arcmin −2 at F1000W based on mid-IR color selections, compared to our ∼5 AGNs per arcmin −2 in the same band.There are two likely reasons for this low yield: (1) their template fitting is based on a single set of templates from Kirkpatrick et al. (2015) in which AGN SEDs are represented by power laws with a narrow range of slopes and thus do not represent the full range of AGN SED variations; (2) color-color plots are less diagnostic than SED fitting, as discussed extensively in Lyu & Rieke 2022a.Consequently, the AGN numbers identified in this way are expected to be less than those from SED analysis.Another limitation of the Kirkpatrick et al. (2023) study is that their templates include strong PAH emission for the normal SFGs and do not capture the IR behavior of lower-metallicity galaxies, which may result in misclassifying some dwarf galaxies with strong warm dust emission as obscured AGNs.

THE OBSCURED AGN POPULATION
Given the AGN sample constructed from this work, we now discuss the obscured AGN population seen at different wavelengths.

The Fraction of Obscured AGNs from SED analysis
We first characterize the obscured AGN fraction of MIRIselected AGNs and explore its dependence on the AGN lu-  minosity and redshift.Our discussion will focus on the 114 MIRI-selected AGNs at L AGN,bol > 10 10 L ⊙ , M * ≳ 10 9.5 M ⊙ and z=0-4.These criteria ensure a reasonably large sample to carry out statistical studies that cover a large range of source properties while minimizing the complications caused by dwarf and high-z AGN candidates where the selections are less robust.
Depending on the observed wavelengths and available spectral features, the definition of what is an obscured AGN varies in the literature.To first order, the obscuration of an AGN, i.e., whether the AGN can be seen or not, at a given band is caused by two major effects: (1) how much of the AGN emission can be diminished by obscuring materials (e.g., a dusty torus); (2) how strongly the host galaxy emission dilutes the AGN light.Our SED fittings provide direct measurements of both effects at rest-frame UV to mid-IR.We define AGN obscuration in three ways: • optically obscured: regardless of the AGN intrinsic luminosity and dust attenuation, if the observed AGN light fraction at rest-frame 0.1-1 µm is less than 20% of the whole galaxy, we count it as an optically obscured AGN; • IR obscured: the near-to mid-IR attenuation level of the AGN component is described by tau ir, which can be transferred to the silicate strength following equation S10 = 0.2 − tau IR /5.5 (negative value for absorption) (Lyu et al. 2022a).We define an AGN as IR obscured if S10 < 0, corresponding to tau ir>1.1; • IR heavily obscured: we define an AGN as IR heavily obscured if tau ir>3.85.This corresponds to the silicate strength of S10 < −0.5, indicating significant dust attenuation (e.g., Goulding et al. 2012).
In a statistical sense, the significance of AGN obscuration increases from optically obscured, IR obscured to IR heavily obscured, though the definition of the first case is complicated by AGN-host galaxy contrast as well as the SMC extinction law.
In Figure 17, we show how the obscured AGN fractions change with different source properties.All three obscuration fractions first increase with AGN bolometric luminosity from L AGN,bol ∼ 10 10 L ⊙ to L AGN,bol ∼ 10 11.5 L ⊙ , with the slope increasing from optical obscured, obscured to IR heavily obscured, and then dropping towards higher luminosity in a consistent way.For redshift evolution, all the obscured AGN fractions increase as a function of object redshift.We discuss the implications of these trends below.
Previous studies have repeatedly shown that the obscured AGN fraction decreases with increasing AGN luminosity (e.g., Lawrence 1991;Simpson 2005;Lusso et al. 2013), so the trend in the left panel of Figure 17 is a surprise at first glance.However, most of these studies have been limited to L AGN,bol ≳ 10 11 L ⊙ .At low-luminosity end, the obscured fraction has been found to increase with AGN luminosity in the X-ray (e.g., Burlon et al. 2011).In particular, Buchner et al. (2015) did show a similar trend with the Xray obscuration level peaking at L X ∼ 10 43.5 erg s −1 and declining both above and below this value around Cosmic Noon (see also, e.g., Brightman & Nandra 2011 for a similar trend reported at low redshifts).Adopting the bolometric luminosity correction of K X ∼ 11 (Duras et al. 2020), this corresponds to L AGN,bol ∼ 10 11 L ⊙ , consistent with the peak luminosity of our SED-inferred obscuration.
From a theoretical perspective, AGN obscuration occurs from the gas inflow that feeds the central SMBH, where the disk becomes unstable due to turbulence that produces a torus-like structure or dusty wind (e.g., Wada 2012Wada , 2015;;Hopkins et al. 2012; see review by Netzer 2015).Based on high-resolution hydrodynamic simulations, the obscuring gas column density presents positive correlations with the AGN luminosity (e.g., Blecha et al. 2018;Trebitsch et al. 2019).When the AGN becomes sufficiently bright, strong feedback is expected to blow out the surrounding material, reducing the chance of obscuration (e.g., Sanders et al. 1988;Hopkins et al. 2008).
Based on the analysis of the rest-frame UV to mid-IR SEDs, our study suggests an increasing fraction of obscured AGNs towards higher redshifts.Previous work in the X-ray based on the gas column density N H has also revealed an increasing fraction of obscured AGNs from the local Universe to Cosmic Noon (e.g., Ueda et al. 2003;La Franca et al. 2005;Treister & Urry 2006;Buchner et al. 2015;Vijarnwannaluk et al. 2022).Such results are possibly associated with the large amount of gas in high redshift galaxies (e.g., Carilli & Walter 2013) as well as the more vigorous growth of SMBHs in the early Universe (Inayoshi et al. 2020).
Finally, AGN selections in our extended and highz samples (dwarf and z > 4 galaxies) are less secure or complete so we only report the integrated numbers of the obscured fractions for reference.For extended sample AGN candidates, the fraction of optically obscured, IR obscured and IR heavily obscured AGN is 86%, 76% and 42%.These numbers are subject, of course, to further confirmation of the AGN nature of individual members of this sample.For the AGN candidates at z > 4, the fraction of optically obscured, IR obscured and IR heavily obscured AGN is very roughly 54%, 75% and 46%.Given that the MIRI data is shallower than that from NIRCam, the obscured AGN fraction among these galaxies may be underestimated.

X-ray Bright but Mid-IR Faint AGNs
As shown in Section 4.2, our MIRI+NIRCam data is much deeper than other multi-wavelength data used for AGN selection in terms of source bolometric luminosity.Naively, one would expect that the superb JWST data can pick out all the AGNs found at other wavelengths.However, this is not the case.Notable examples are 23 X-ray bright AGNs not identified by SED fittings.In other words, these systems are X-ray detected bright AGNs but very faint in the optical to the mid-IR.
In Figure 18, we illustrate the locations of such objects on the X-ray luminosity vs the rest-frame 6 µm AGN luminosity plot.Previous studies have established a strong correlation between the absorption-corrected X-ray luminosity and the AGN luminosity at 6 µm for most AGNs (e.g., Asmus et al. 2015;Stern 2015;Mateos et al. 2015;Chen et al. 2017).For AGNs identified by both MIRI SED fitting and X-rays, such a trend also exists despite a large scatter.However, there are many X-ray bright AGNs with very weak AGN 6µm emission (left side of Figure 18).In fact, even if we ignore the SED decomposition and assume all the observed flux is from the AGN, eleven objects cannot be moved near the typical L X,int − L 6µm,AGN relation.We have carefully inspected these objects.Besides 8 objects with noisy SEDs or degenerate fitting results, 15 do not show any AGN evidence in the UV-to-mid-IR SEDs.In Figure 19, we present some example SEDs with the cutouts from NIRCam, MIRI and Chandra.Besides the unresolved SED analysis, all these systems are extended in both NIRCam and MIRI images, consistent with the lack of AGN signatures.Meanwhile, their emission is very strong in the X-ray.Similar AGNs have been reported previously (LaMassa et al. 2019;Lyu et al. 2022a).
One possible explanation could be dust-deficient AGNs that lack the hot/warm dust emission (Lyu et al. 2017).If we limited the sample to X-ray detected AGNs, the fraction of such objects is about 16-25%, consistent with the typical fraction of dust deficient AGN in an IR-unbiased sample as reported in Lyu et al. (2017).Another possibility is that the X-ray emission of these AGNs is boosted by, e.g., jets.Typically, the X-ray jet emission is associated with radio loud AGNs.Among these 23 X-ray bright AGNs, 15 are detected by the JVLA at 3 GHz -one is a confirmed radio-loud AGN, two are ambiguous cases and 12 are consistent with typical SFGs.It is therefore likely that the X-ray bright but mid-IR faint AGNs arise from both the radio-loud and radio-quiet populations.
Lastly, as pointed out in Section 4.3, many [X2R] AGNs are not identified via JWST SEDs (see also Figure 16).These AGNs have been picked out by higher X-ray to radio luminosity ratios than is expected from stellar processes in a galaxy.It is possible that they are similar to the X-ray bright and IR weak AGNs described above but with even lower AGN power.

MIRI AGNs without X-ray Detections
As described in Section 4.2, a large fraction of our MIRI AGNs do not have X-ray detections (see also Figure 18).The intrinsically faint AGNs will not be detected in X-ray even if they are unobscured due to the shallower Chandra detection limits compared to MIRI.For relatively bright AGNs, the X-ray non-detected AGNs can be interpreted as being Compton-thick or intrinsically X-ray weak.We now focus on the first possibility and put some constraints on the Comptonthick AGN fractions with the aid of the MIRI SED results.
Since our SED fittings provide a measurement of the AGN bolometric luminosity, we can convert them to the equivalent X-ray luminosity for different gas column densities and compare with observations.In Figure 20, we show the ratio between the SED-derived bolometric luminosity and the Xray luminosity at 2-10 keV against the AGN bolometric luminosity for MIRI-identified AGNs at different redshifts.We compute an upper limit to the 2-10 keV luminosity assuming no extinction for objects without X-ray detections.We also plot the bolometric correction for the extinction-corrected AGN X-ray luminosity at 2-10 keV from Duras et al. (2020) as well as the expected correlation if the X-ray emission is Compton-thick with gas column density N H ∼ 10 24 cm −2 .The relative location of an AGN on this diagram is determined by its intrinsic luminosity as well as the gas column density.AGNs above the red curve can be counted as Compton-thick AGNs.
The weak X-ray emission of most AGNs at relatively low redshift is well-explained in terms of absorption (Saade et al. 2022;Wang et al. 2022).Above z ∼1, we found the obscured AGN fraction gradually increases with redshift, possibly due to an increase in the number of Compton thick sources.Our finding of ∼ 34% X-ray-undetected AGNs with SEDs suggestive that they are Compton thick aligns well with the typical explanation of diffuse X-ray background (e.g., Akylas et al. 2012).Although there are some variations, a summary of its properties is: (1) 92.7 ± 13.3 % of the hard X-ray (2-7 keV) background is resolved into individual sources by the Chandra 7 Msec exposure (Luo et al. 2017;Cappelluti et al. 2017); (2) models to fit the background are centered on the AGN population at redshifts of 0.5 -2; and (3) many of these models invoke a population of Compton thick AGNs not detected in X-rays to fit the background around 30 keV (e.g., Gilli et al. 2007;Moretti et al. 2009;Treister et al. 2014).It  (Luo et al. 2017), the intrinsic X-ray luminosity log(LX,int/erg s −1 ) and the corresponding gas column density log(NH /cm −2 ) of each object are: MIRI 419 -41.58, 23.47;MIRI 757 -42.70, 23.65;MIRI 779 -42.72, 21.84;MIRI 1382 -43.12, 22.46.should be possible to advance our understanding of the Xray background significantly through deep JWST surveys for infrared-discovered AGNs, although such studies are beyond the scope of this paper.
Although the increase of the fraction of X-ray undetected AGN population apparently continues for z >3, Compton absorption becomes weak at the relevant rest-frame energies12 .Rather than obscuration, some of these AGNs may be intrinsically weak in the X-ray.Indeed, the existence of such objects has been suggested previously by, e.g., Leighly (2004); Simmonds et al. (2016); Lyu et al. (2022a).Stacking the X-ray observations for all six of the AGNs in Figure 9 yields non-detections of 3.6 ± 5.0 × 10 −18 erg cm −2 s −1 for the 0.5 − 2 keV band and −1.13±1.28×10−17 erg cm −2 s −1 for the 2 − 7 keV band.These stacked values reinforce the result in Figure 20 that the high redshift AGNs are in general very X-ray weak.
We now explore the implications of this result.The Xray output of high redshift quasars in general does not show any deficiency of X-ray emission with an average α OX of ∼ −1.5 for those with measured values (omitting two outliers), using the results from Li et al. (2021).This does not exclude weak X-rays from the undetected members, which are the majority of the sample, although these authors con-clude that the relation of the X-ray outputs to the other parameters of the quasars is similar to that for lower redshift quasars.The low X-ray fluxes for the lower luminosity highz AGNs could arise due to Compton thick absorption with n H ≳ 10 25 cm −2 , but such high absorbing columns are very rare locally.The infrared AGN-driven SED shapes for the six sources in Figure 9, extending down to 1-2 µm and even shorter wavelengths, are similar to those of typical local Type-2 Seyfert galaxies, rather than extreme Compton thick examples like NGC 1068.A hint toward an alternative explanation is provided by the very X-ray weak AGNs in a number of low-metallicity dwarf galaxies (Dong et al. 2012;Simmonds et al. 2016;Cann et al. 2021;Bohn et al. 2021;Latimer et al. 2021), often two orders of magnitude below the values expected by analogy with AGNs in massive galaxies.That is, the behavior might be associated with low metallicity in modest mass high redshift host galaxies.
Hard X-rays from AGNs are produced by Compton scattering in the coronal regions above their accretion disks.Therefore, a process that disrupts the coronal regions could suppress the X-ray outputs.The dwarf galaxies have strong AGN-driven outflows that are believed to be effective at clearing gas from the galaxies (Bohn et al. 2021); these outflows are also candidates to disrupt the coronal regions.A possible issue with this explanation is that a small sample indicates that the outflow and coronal line strengths may be correlated (Bohn et al. 2022).Overall, this is a complex situ-  ation requiring much larger samples than provided in this paper to test, particularly since AGN variability is likely to influence the observational properties of any sample of AGNs.
Besides the possibilities discussed above, the lack of Xray detections for high-z AGNs might be also explained by the possible evolution of the X-ray spectral shape.For example, Zappacosta et al. (2023) found the X-ray spectra of ten quasars at z > 9 have steeper average photon index (Γ ≈ 2.4 ± 0.1) than classical values (Γ ∼1.8-2).For such systems, even their X-ray luminosity can be high, their redshifted spectrum at the commonly observed X-ray energy band could be a factor up to 4-5 times fainter than standard AGNs, which can greatly reduce the chance of detection.

THE BLIND MEN AND THE ELEPHANT: HOW TO BUILD A COMPLETE AGN SAMPLE?
As demonstrated in this work, JWST has made a huge step forward in finding AGNs, including those at very high redshift, highly obscured, or residing in low-mass galaxies.However, this does not mean the current AGN census is complete.As an epilogue, we discuss the lessons learned so far and give some future outlook.
The process of selecting AGNs is like the story of the blind men and the elephant.Every study is limited by the wavelength regime and survey depth, so the results are inevitably biased and incomplete.In fact, the AGN phenomena is more complicated than elephants and astronomers typically work on different samples built at different wavelengths.Due to e.g., the intrinsic variation of AGN properties, the apparent obscuration caused by dust and gas, the efficiency of the selection method for AGN, and the data quality that the selection has to be based on, AGNs seen in one study can be totally missed in another study (e.g., see Section 4 on the existence of X-ray bright AGNs without relevant IR signatures in the JWST SEDs and bright mid-IR AGNs without X-ray detections).Thus, a complete picture can be only obtained by combining the results across the electromagnetic spectrum together.
In terms of selection techniques, the so-called gold standards for AGN identification, such as very broad emission lines, bright X-ray emission or clear hot dust signatures, are definitive but certainly not complete indicators.Given how complicated and diverse the AGN population is, we could miss a substantial fraction of AGNs if the selections rely on these most obvious features, as seen in Section 5. When multi-wavelength data is available, it is reasonable to relax the criteria and use multiple (relaxed) selections to complete the AGN census.From the AGN hints from other wavelengths, selections at one band can be calibrated and im-proved (e.g., Donley et al. 2012;Hviding et al. 2022).However, each selection technique has limitations.For example, we should not adopt the classical AGN color criteria or SED fittings with normal galaxy templates to look for AGN in the dwarf galaxy population (e.g., Hainline et al. 2016).As we enter a new regime, relevant tools need to be updated or developed.Lastly, besides refining the AGN selection methods and expanding their usage, there must also be a balance between robustness and completeness.
Given the various factors that influence the observed AGN behavior, whether an AGN sample is complete is actually determined by how we define the term "AGN" observationally, which can be very controversial (e.g., the extended discussion about how to classify LINERs, now understood to be a mixed category).Our ultimate goal is not to collect all kinds of AGN samples but to understand the physics behind the AGNs themselves and the SMBH-galaxy interaction.In this sense, one intriguing future direction is to combine the observed AGN statistics with the predicted observables from cosmological simulations that involve reasonable treatments of various selection biases.This approach has more physics behind it and can mitigate the issues of obscured AGNs that may not be easily revealed by the data.If one day we could build such a physical framework that matched all the key empirical correlations between SMBHs and galaxies with good predictive power for the observed AGN statistics at different wavelengths, our long journey to complete the AGN census might reach its final end.

SUMMARY
In this work, we presented the first results on AGN selection and demographics from SMILES, a JWST Cycle 1 GTO program that has targeted the central region of GOODS-S with eight MIRI filters from 5.6 to 26 µm.Combining the MIRI data with JADES NIRCam and HST photometry at shorter wavelengths, we conducted comprehensive SED analyses of 3273 MIRI-sources and reported 217 AGNs over a footprint of 34 arcmin 2 .This MIRI-selected AGN sample includes 111 AGNs in the primary sample of normal massive galaxies, 86 AGN candidates in the extended sample of dwarf galaxies and 20 high-z AGN candidates at z > 4. Our major conclusions are: • We reached a total number density of SED-identified AGNs of ∼ 6.6 arcmin −2 , which is over two times larger than that for the confusion-limited Spitzer/MIPS 24 µm observations of the same field (Lyu et al. 2022a).In addition, compared to the deepest X-ray survey by Chandra in the same footprint, our relatively shallow MIRI surveys (∼2.17 hours per pointing) yield 2.4 times more AGNs, of which 80% are new discoveries.
• AGNs in our primary sample of massive galaxies should be sufficiently luminous for detection by previous searches in other wavebands, assuming they are not obscured.However, 34% of the AGNs found in our study do not have previous identifications.This result indicates that 34 ± 6% of AGNs have been missed due to obscuration, where the error arises from the modest sample size (114); • For the first time at z ∼ 2, JWST/MIRI is detecting significant numbers of AGNs in low-mass galaxies [log(M * /M ⊙ ) < 9.5] -from our conservative SED identifications, they are comparable in number with the AGNs we have found in the more massive primary sample; • We confirm the existence of X-bright but mid-IR faint AGNs and suggest that they may be dust-deficient AGNs that do not present strong hot dust SED features, or AGNs with X-ray emission boosted by e.g. a jet, or a mixture of both populations; • ∼80% of the MIRI-selected AGNs do not have X-ray detections.Most of them are too faint to be detected by Chandra, some of which are likely to be Comptonthick AGNs.For MIRI AGNs without X-ray detections at higher redshifts (z ∼4-6), they may be intrinsically weak in the X-ray, possibly due to the lack of a strong corona component due to their lower metallicity environment.
• Based on SED analysis, we find the obscured AGN fraction increases with AGN luminosity towards L AGN ∼ 10 11 L ⊙ and then drops at higher luminosity.The obscured AGN fraction also increases with redshift, regardless of the detailed definitions.For AGNs at z > 4, it is likely half of the population is heavily obscured.
• Despite the comprehensive MIRI photometric data and the substantially improved SED selections, about 28% of the known AGN (candidates) in the field have been missed by the SED analysis, indicating a huge diversity of AGN properties that challenges the completion of the AGN census.In other words, no single method or wavelength can identify all the AGNs and a combination of multi-wavelength dataset and selection techniques is always desired; Our results demonstrate the unique power of JWST MIRI for AGN selection and characterization.This paper summarizes first results from a modest initial survey and we look forward to larger and more ambitious studies.a 1 = from IR SED; 2 = from IR colors; 3 = from X-ray luminosity; 4 = from X-ray to radio flux ratio; 5 = from radio loudness; 6 = from optical spectrum; 7 = from variability Figure1.Survey layout from X-ray to radio (left) and the 5σ flux limits of photometric bands at 0.3 -26 µm (right) used for SED fittings in GOODS-S/HUDF.On the left, the footprint of our MIRI survey is shown as the red-shaded region with JADES NIRCam deep (solid green line) and medium (dashed green line) observations, FRESCO NIRCam/grism coverage (solid blue line), HST ACS GOODS-S coverage (thick gray line), Chandra X-ray coverage (light blue shaded region) and JVLA radio observations at 3 GHz (dark green lines) and 6 GHz (dark orange lines).On the right, we denote the pre-JWST filters in gray, JWST/NIRCam filters in green and JWST/MIRI filters in red.Besides the flux limits of these filters as a function of wavelength, we also show some representative SEDs of obscured AGNs (orange), unobscured AGNs (blue) and starburst galaxies (magenta) at z=2 and z=4, where JWST observations are expected to bring major advances to identify and characterize these objects.

Figure 2 .
Figure 2. Showcase of multi-wavelength data used for AGN searches in SMILES with a field of view 3.5 arcmin×1.0arcmin.The yellow circles indicate the AGNs or AGN candidates identified from this work.The JWST NIRCam and MIRI data are shown as three-color images with the long wavelength in red, intermediate wavelength in green and short wavelength in blue.

Figure 3 .
Figure3.Redshift distribution of MIRI sources.We highlight the sources of the final adopted redshifts with different colors.

FFFigure 4 .
Figure4.Illustration of how the AGN SED shape is changed by different levels of attenuation as well as the templates for the galaxy stellar and dust components in our model.The upper panel shows the AGN SED obscured by the SMC extinction curve, mainly used for UV-optical reddened AGNs.The lower panel presents the AGN SEDs modified by the empirical IR attenuation law.We show a galaxy stellar SED with a stellar age of 500 Myr (light green thick line), and the pure dust emission templates for normal star-forming galaxies (light orange line) and low-metallicity dwarf galaxies (light magenta line) to compare with the AGN SEDs in both panels.See text for details.

Figure 5 .
Figure5.Example JWST NIRCam+MIRI SED fittings of dusty star-forming galaxies (left) and IR-quiescent galaxies (right).The object redshifts increase from top to bottom.Below each SED plot, we also show 3 ′′ ×3 ′′ cutouts of the source NIRCam and MIRI images and comparisons of the normalized source light profile (black line) and the corresponding ePSF (red line).( The x-axis of the light profile panels spans from 0 to 6.0 arcsec.)

Figure 6 .
Figure 6.Similar to Figure 5 but for AGNs.In the top six panels, we show normal AGNs with a range of nuclear obscuration and host galaxy contamination at z < 4. The four panels in the bottom row are example AGNs in low-mass galaxies.

Figure 7 .
Figure7.The stellar mass as a function of redshift for MIRI sources.We highlight all AGNs identified with the normal SFG template regardless of the object redshift and stellar mass in light blue on the left, and the refined AGN sample after adopting the Haro 11 dwarf SFG template to reject possible false positives among low-mass galaxies in dark blue on the right.In both two panels, we show the normal SFG galaxies as orange crosses and other galaxies as gray dots.The red solid line represents the stellar mass threshold for 12+log(O/H)∼8.2 with the ±0.3 dex variation range as red dashed line.If we adopt the same normal galaxy template for all stellar masses (left panel), there would be 500 AGN candidates selected.In contrast, we end up with 217 AGN candidates after considering the dwarf galaxy complication (right panel), as described in Sections 3.2.3 and 3.3.2

Figure 8 .
Figure 8. Example SED fittings of a low-mass galaxy.The top panel shows the object is identified as an AGN if we use the typical templates for normal SFGs and the bottom panels shows the fittings with the Haro 11 template, where the AGN evidence is gone.

Figure 9 .
Figure 9. Similar to Figure 5 but for high-z AGN candidates.The limited rest frame infrared coverage still allows identification of obscured AGNs, although not as completely as for lower redshifts.

Figure 10 .
Figure 10.AGN luminosity as a function of redshift.The MIRI AGNs are shown in blue and the pre-JWST AGNs are shown in gray.

Figure 11 .
Figure 11.The observed IR-to-radio flux q 21,obs as a function of redshift of all the sources.The red dotted line is 0.5 dex below the mid-point of the radio-infrared correlations in the Rieke et al. (2009) SFG templates (the green solid line); we define the source as a radio-loud AGN if it falls below this line at ≥ 2-σ significance.We color-code different sources on the main plot and show their relative distribution of q 21,obs in the right panel.

Figure 12 .
Figure 12.Comparison of the SED fittings with pre-JWST data (left) and JWST data (right).Two examples are shown.

Figure 13 .
Figure 13.The MIRI source number density and AGN fraction as a function of observed flux in each MIRI band.

Figure 14 .
Figure14.Venn diagram of the various AGN samples within the SMILES footprint.We compare the X-ray detected AGN sample inLuo et al. (2017), the pre-JWST AGN sample compiled inLyu et  al. (2022a)  and the MIRI SED-identified AGN sample in this work.The numbers in the brackets are the total numbers of AGN from each source.For objects identified by MIRI, we also plot a horizontal stacked bar chart to show the number of normal AGNs (green), dwarf AGNs (red) and high-z AGNs (yellow) for each subset.

Figure 15 .
Figure 15.L AGN,bol detection limits as a function of redshift for different bands.

Figure 16 .
Figure16.Comparison of AGN sample selected with different methods and their overlaps in the SMILES footprint, visualized with the UpSet technique.The sample intersections of AGNs selected with different methods are plotted as a matrix with bar charts on the left and top to show the corresponding sample sizes.Each row corresponds to one selection method, and the bar charts on the left show the number of AGN identified with a given method.Each column corresponds to a possible intersection: the filled-in cells indicate which selection method is part of the intersection.The lines connecting the filled-in cells show in which direction the plots should be read.The bar charts on the top give the size of the AGN sample identified with the corresponding intersection of various selection techniques.We further break each vertical bar into different colors to show the relative fraction of normal, dwarf and high-z AGNs within each subset.We highlight the subsets contributed by JWST SED selections in red.

Figure 17 .
Figure 17.The obscured AGN fraction as a function of AGN luminosity (left) and object redshift (right).See text for the details.

Figure 18 .
Figure18.The relation between X-ray luminosity and AGN continuum luminosity at rest-frame 6 µm for AGN in the SMILES footprint.

Figure 20 .
Figure 20.Bolometric to X-ray luminosity fraction as a function of AGN luminosity at different redshifts for MIRI-identified AGNs.On the bottom of each panel, we denote the total number of AGNs as N and Compton-thick AGN fraction fCT for z < 3.0 and X-ray Intrinsic Weak AGN fraction fXIW for z > 3.0.

Table 1 .
Parameter Setup of the Modified Prospector Code

Table 2 .
List of MIRI-selected AGN in the SMILES Footprint classification code for different selection methods.'o' is given if the object is classified as AGN by this method and 'x' means not.From NOTE-Only a portion of this table is shown here to demonstrate its form and content.A machine-readable version of the full table is available.

Table 3 .
AGNs identified from weak Pa α