# CEERS Key Paper. VI. JWST/MIRI Uncovers a Large Population of Obscured AGN at High Redshifts

G. Yang, K. I. Caputi, C. Papovich, P. Arrabal Haro, M. B. Bagley, P.

Behroozi, E. F. Bell, L. Bisigello, V. Buat, D. Burgarella, et al.

► To cite this version:

G. Yang, K. I. Caputi, C. Papovich, P. Arrabal Haro, M. B. Bagley, et al. CEERS Key Paper. VI. JWST/MIRI Uncovers a Large Population of Obscured AGN at High Redshifts. The Astrophysical journal letters, 2023, 950, 10.3847/2041-8213/acd639. insu-04479080

HAL Id: insu-04479080

https://insu.hal.science/insu-04479080

Submitted on 5 Mar 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

CEERS Key Paper. VI. JWST/MIRI Uncovers a Large Population of Obscured AGN at
High Redshifts

G. Yang (杨光)<sup>1,2</sup>①, K. I. Caputi<sup>1,3</sup>①, C. Papovich<sup>4,5</sup>①, P. Arrabal Haro<sup>6</sup>①, M. B. Bagley<sup>7</sup>①, P. Behroozi<sup>8,9</sup>①, E. F. Bell<sup>10</sup>①,
L. Bisigello<sup>11,12</sup>①, V. Buat<sup>13</sup>①, D. Burgarella<sup>13</sup>①, Y. Cheng<sup>14</sup>①, N. J. Cleri<sup>4,5</sup>①, R. Davé<sup>15,16</sup>①, M. Dickinson<sup>6</sup>①,
D. Elbaz<sup>17</sup>①, H. C. Ferguson<sup>18</sup>①, S. L. Finkelstein<sup>7</sup>①, N. A. Grogin<sup>19</sup>①, N. P. Hathi<sup>18</sup>①, M. Hirschmann<sup>20</sup>①,

B. W. Holwerda<sup>21</sup>①, M. Huertas-Company<sup>22,23,24</sup>①, T. A. Hutchison<sup>25,34</sup>①, E. Iani<sup>1</sup>①, J. S. Kartaltepe<sup>26</sup>①, A. Kirkpatrick<sup>27</sup>①,
D. D. Kocevski<sup>28</sup>①, A. M. Koekemoer<sup>19</sup>①, V. Kokorev<sup>1</sup>①, R. L. Larson<sup>29,35</sup>①, R. A. Lucas<sup>19</sup>①, P. G. Pérez-González<sup>30</sup>①,
P. Rinaldi<sup>1</sup>①, L. Shen<sup>4,5</sup>①, J. R. Trump<sup>31</sup>①, A. de la Vega<sup>32</sup>①, L. Y. A. Yung<sup>25</sup>①, and J. A. Zavala<sup>33</sup>①

Abstract

Mid-infrared (mid-IR) observations are powerful in identifying heavily obscured active galactic nuclei (AGN) that have weak emission in other wavelengths. Data from the Mid-Infrared Instrument (MIRI) on board the James Webb Space Telescope provides an excellent opportunity to perform such studies. We take advantage of the MIRI imaging data from the Cosmic Evolution Early Release Science Survey to investigate the AGN population in the distant universe. We estimate the source properties of MIRI-selected objects by utilizing spectral energy distribution (SED) modeling, and classify them into star-forming galaxies (SFs), SF-AGN mixed objects, and AGN. The source numbers of these types are 433, 102, and 25, respectively, from four MIRI pointings covering ~9 arcmin<fnref n="2" />. The sample spans a redshift range of ≈0–5. We derive the median SEDs for all three source types, respectively, and publicly release them. The median MIRI SED of AGN is similar to the typical SEDs of hot dust-obscured galaxies and Seyfert 2s, for which the mid-IR SEDs are dominated by emission from AGN-heated hot dust. Based on our SED-fit results, we estimate the black hole accretion density (BHAD; i.e., total BH growth rate per comoving volume) as a function of redshift. At z < 3, the resulting BHAD agrees with the X-ray measurements

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands; gyang206265@gmail.com</footnote>

<footnote n="2">² SRON Netherlands Institute for Space Research, Postbus 800, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="3">³ Cosmic Dawn Center (DAWN), Copenhagen, Denmark</footnote>

<footnote n="4">⁴ Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA</footnote>

<footnote n="5">⁵ George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA</footnote>

<footnote n="6">⁶ NSF's National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="7">⁷ Department of Astronomy, The University of Texas at Austin, Austin, TX, USA</footnote>

<footnote n="8">⁸ Department of Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="9">⁹ Division of Science, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="10">¹⁰ Department of Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI 48109-1107, USA</footnote>

<footnote n="11">¹¹ Dipartimento di Fisica e Astronomia “G.Galilei,” Università di Padova, Via Marzolo 8, I-35131 Padova, Italy</footnote>

<footnote n="12">¹² INAF—Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy</footnote>

<footnote n="13">¹³ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="14">¹⁴ University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA</footnote>

<footnote n="15">¹⁵ Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh, EH9 3HJ, UK</footnote>

<footnote n="16">¹⁶ Department of Physics and Astronomy, University of the Western Cape, Robert Sobukwe Road, Bellville, Cape Town 7535, South Africa</footnote>

<footnote n="17">¹⁷ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="18">¹⁸ Space Telescope Science Institute, Baltimore, MD, USA</footnote>

<footnote n="19">¹⁹ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="20">²⁰ Institute of Physics, Laboratory of Galaxy Evolution, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland</footnote>

<footnote n="21">²¹ Physics & Astronomy Department, University of Louisville, Louisville, KY 40292, USA</footnote>

<footnote n="22">²² Instituto de Astrofísica de Canarias, La Laguna, Tenerife, Spain</footnote>

<footnote n="23">²³ Universidad de la Laguna, La Laguna, Tenerife, Spain</footnote>

<footnote n="24">²⁴ Université Paris-Cité, LERMA—Observatoire de Paris, PSL, Paris, France</footnote>

<footnote n="25">²⁵ Astrophysics Science Division, NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA</footnote>

<footnote n="26">²⁶ Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA</footnote>

<footnote n="27">²⁷ Department of Physics and Astronomy, University of Kansas, Lawrence, KS 66045, USA</footnote>

<footnote n="28">²⁸ Department of Physics and Astronomy, Colby College, Waterville, ME 04901, USA</footnote>

<footnote n="29">²⁹ The University of Texas at Austin, Department of Astronomy, Austin, TX, USA</footnote>

<footnote n="30">³⁰ Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850 Madrid, Spain</footnote>

<footnote n="31">³¹ Department of Physics, 196 Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269, USA</footnote>

<footnote n="32">³² Department of Physics and Astronomy, University of California, 900 University Ave, Riverside, CA 92521, USA</footnote>

<footnote n="33">³³ National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote>Received 2023 March 21; revised 2023 May 13; accepted 2023 May 15; published 2023 June 12</footnote>

<footnote n="34">³⁴ NASA Postdoctoral Fellow.</footnote>

<footnote n="35">³⁵ NSF Graduate Fellow.</footnote>
in general. At $z > 3$, we identify a total of 27 AGN and SF-AGN mixed objects, leading to that our high-$z$ BHAD is substantially higher than the X-ray results (~0.5 dex at $z \approx 3-5$). This difference indicates MIRI can identify a large population of heavily obscured AGN missed by X-ray surveys at high redshifts.

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); High-redshift galaxies (734); Supermassive black holes (1663); Star formation (1569)*

*Supporting material: data behind figure*

# 1. Introduction

Active galactic nuclei (AGN) are powered by growing supermassive black holes (BHs) typically harbored at the centers of galaxies. AGN luminosities can be used to estimate BH growth rates. From studying the AGN demography at different redshifts, the evolution history of the BH population can be inferred (e.g., Merloni & Heinz 2008; Yang et al. 2018; Aird et al. 2019, 2022; Ananna et al. 2019). To depict an unbiased picture of the BH cosmic evolution, it is critical to conduct a complete census of the AGN population minimizing selection biases.

Multiwavelength techniques of AGN selection have been intensively developed in the past ~60 yr since their discovery (e.g., Padovani et al. 2017; Hickox & Alexander 2018; Lyu et al. 2022). Among these techniques, optical AGN selection is often the most straightforward and can be performed with ground-based telescopes. Color-color diagrams based on broadband photometry are often employed to select AGN candidates (e.g., Richards et al. 2002). Follow-up spectroscopic analyses identifying broad lines or using line ratios (e.g., Kauffmann et al. 2003) can verify the AGN nature. However, optical color-color selections are biased toward selecting relatively unobscured type 1 AGN. Compared to optical selections, X-ray AGN selections are more complete in general. Strong X-ray emission is ubiquitous among the AGN population, and X-ray photons, especially hard X-rays, have excellent penetrating power. These advantages allow X-ray methods to select not only unobscured AGN but also obscured ones (e.g., Brandt & Alexander 2015; Brandt & Yang 2022). However, even the hardest X-ray photons with $E > 10$ keV can be Compton scattered and absorbed when the obscuration is heavy with column density $N_H \gtrsim 10^{24} \text{ cm}^{-2}$. Indeed, evidence suggests that X-ray surveys might miss many heavily obscured AGN, especially at high redshifts ($z \gtrsim 3$) where the AGN obscuration is generally strong, because gas tends to be more abundant and compact toward high redshift (e.g., Yang et al. 2021a; Gilli et al. 2022).

Infrared (IR) observations have the distinct advantage of selecting obscured AGN. Obscured AGN are often dust-rich, and these dust particles are heated to ~100–1000 K by the emission from the central engine of the accretion disk. The dust reemits the “waste heat” primarily in the mid-IR wavelengths of ~3–30 µm. By detecting this mid-IR radiation, AGN can be identified. Color-color methods based on Spitzer/Infrared Array Camera (IRAC) have been developed and widely applied to IRAC surveys (e.g., Stern et al. 2005; Lacy et al. 2007; Donley et al. 2012). Studies based on Spitzer data have also worked on the selection of AGN based on the presence of a power law in the IRAC bands (e.g., Alonso-Herrero & Pérez-González 2006; Caputi 2013). These techniques have successfully identified numerous obscured AGN, some of which are missed by even the deepest X-ray surveys (e.g., Lyu et al. 2022). However, due to the limitations of IRAC sensitivity and

wavelength coverage (below 8 µm), the techniques are ineffective in identifying extremely faint AGN, especially at high redshifts, where AGN radiation is redshifted beyond the IRAC wavelength coverage (Mendez et al. 2013; Lyu & Rieke 2022). The Multiband Imaging Photometer for Spitzer (MIPS) and Herschel cover longer wavelengths, but they are not comparably suitable for distant-AGN searches due to their lower sensitivities compared to IRAC.

The Mid-Infrared Instrument (MIRI) on board the recently launched James Webb Space Telescope (JWST) has the potential to become a game changer in the detection and selection of AGN especially at high redshifts. MIRI has continuous photometry coverage of ~5–25 µm, reaching an unprecedented sensitivity level about 10–100 times deeper than Spitzer. Besides, thanks to the large aperture size of JWST, the angular resolution of MIRI reaches an FWHM of point-spread functions being ~0.′′2–0.′′8, about 8 times better than Spitzer/IRAC at similar wavelengths. Indeed, simulations suggest that MIRI is able to capture AGN activity down to a faint level of, e.g., an Eddington ratio around 0.01 at $z \sim 1–2$ (e.g., Kirkpatrick et al. 2017; Yang et al. 2021b).

In this Letter, we take advantage of the MIRI imaging data from the Cosmic Evolution Early Release Science Survey (CEERS; Finkelstein et al. 2017), which targets the Extended Groth Strip (EGS) field. We also compile the multiwavelength data available in EGS and perform AGN versus star-forming (SF) classifications with spectral energy distribution (SED) modeling. We produce median SEDs for different source types, respectively, and publicly release these empirical SEDs to facilitate future science. Finally, we investigate the BH accretion density (BHAD) as a function of redshift, based on the MIRI-identified AGN.

The structure of this Letter is as follows. In Section 2, we describe and analyze the MIRI imaging data. We also perform AGN selections based on the MIRI and other multiwavelength photometric data. In Section 3, we discuss the physical implications based on our analysis results. We summarize our results and discuss future prospects in Section 4.

Throughout this Letter, we assume a cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$. We adopt a Chabrier initial mass function (IMF; Chabrier 2003). Quoted uncertainties are at the $1\sigma$ (68%) confidence level. All magnitudes are in AB units (Oke & Gunn 1983), where $m_{AB} = -48.6 - 2.5 \log(f_\nu)$ for $f_\nu$ in units of erg s$^{-1}$ cm$^{-2}$ Hz$^{-1}$.

# 2. Data and Analysis

## 2.1. MIRI and Other Multiwavelength Photometry

The CEERS survey has eight MIRI pointings: four “blue” and four “red.” Each blue pointing has two bands of F560W and F770W. Two of the red pointings have six bands from F770W to F2100W, and the other two have four bands from F1000W to F1800W. In this Letter, we focus on the four red pointings, as the blue pointings lack the long mid-IR
<table><caption>**Table 1**
CIGALE Model Parameters</caption><thead><tr><th>Module</th><th>Parameter</th><th>Symbol</th><th>Values</th></tr></thead><tbody><tr><td rowspan="2">Star formation history<br>sfhdelayed</td><td>Stellar e-folding time</td><td>τ<sub>star</sub></td><td>0.5, 1, 2, 3, 4, 5 Gyr</td></tr><tr><td>Stellar age</td><td>t<sub>star</sub></td><td>0.1, 0.2, 0.5, 1, 2, 3, 4, 5 Gyr</td></tr><tr><td rowspan="2">Simple stellar population<br>bc03</td><td>Initial mass function</td><td>...</td><td>Chabrier (2003)</td></tr><tr><td>Metallicity</td><td>Z</td><td>0.02</td></tr><tr><td rowspan="2">Nebular emission<br>nebular</td><td>Ionization parameter</td><td>log U</td><td>-2.0</td></tr><tr><td>Gas metallicity</td><td>Z<sub>gas</sub></td><td>0.02</td></tr><tr><td rowspan="2">Dust attenuation<br>dustatt_modified_starburst</td><td>Color excess of nebular lines</td><td>E(B-V)<sub>line</sub></td><td>0, 0.02, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, 0.9, 1.0</td></tr><tr><td>Ratio between line and continuum E(B-V)</td><td>E(B-V)<sub>line</sub>/E(B-V)<sub>cont</sub></td><td>1</td></tr><tr><td rowspan="3">Galactic dust emission<br>dl2014</td><td>PAH mass fraction</td><td>q<sub>PAH</sub></td><td>0.47, 2.5, 7.32</td></tr><tr><td>Minimum radiation field</td><td>U<sub>min</sub></td><td>0.1, 1.0, 10, 50</td></tr><tr><td>Fraction of PDR emission</td><td>γ</td><td>0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.9</td></tr><tr><td rowspan="4">AGN (UV-to-IR) emission<br>skirtor2016</td><td>Average edge-on optical depth at 9.7 µm</td><td>τ<sub>9.7</sub></td><td>3.5, 7.9, 11</td></tr><tr><td>Viewing angle</td><td>θ<sub>AGN</sub></td><td>70°</td></tr><tr><td>AGN contribution to IR luminosity</td><td>frac<sub>AGN</sub></td><td>0, 0.01, 0.03, 0.05, 0.1–0.9 (step 0.1), 0.99</td></tr><tr><td>Wavelength range where frac<sub>AGN</sub> is defined</td><td>λ<sub>AGN</sub></td><td>3–30 µm</td></tr><tr><td>Redshift+IGM<br>redshifting</td><td>Source redshift</td><td>z</td><td>0.01–8.0 (50 steps)</td></tr></tbody></table>

Note. For parameters not listed here, default values are being adopted.

wavelength coverage for AGN identification (Section 1). The MIRI photometry is extracted with T-PHOT (Merlin et al. 2015), using Hubble Space Telescope (HST)/F160W as the high-resolution prior. The details of the MIRI observations, data reduction, photometry extraction, and quality assessment will be presented in a dedicated paper (G. Yang et al. 2023, in preparation), and some methods are discussed in Yang et al. (2021b) and Papovich et al. (2023).

We also compile other multiwavelength photometric data, specifically the EGS catalog by Stefanon et al. (2017) from the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011). These data include 17 broadband photometries from the U band to the IRAC 8 µm.

Because this Letter is focused on the MIRI properties of galaxies, we restrict our sample to objects with at least two MIRI bands with detections of signal-to-noise ratio (S/N) above 3. This leads to a parent sample of 560 sources. This criterion guarantees that we have at least one robust MIRI color for each source. In general, shorter-wavelength MIRI bands have higher detection rates thanks to their relatively deeper sensitivity. For example, the S/N > 3 detection rates of F770W and F2100W are 97% and 68%, respectively, among our sources.

## 2.2. SED Modeling Setup

We employ CIGALE v2022.1 (Boquien et al. 2019; Yang et al. 2020, 2022) to perform SED modeling based on the broadband photometry described in Section 2.1. The adopted CIGALE parameters are summarized in Table 1.

For star formation history (SFH), we use a standard delayed-τ module (sfhdelayed in CIGALE). We allow the e-folding time and stellar age varying from 0.5–5 Gyr to 1–5 Gyr, respectively. We employ Bruzual & Charlot (2003; bc03) for the simple stellar population (SSP) module, assuming a

Chabrier (2003) IMF with a solar metallicity (Z = 0.02). We also include the nebular module (Villa-Vélez et al. 2021) for emission from the HII regions. For the dust attenuation, we adopt dustatt_modified_starburst in CIGALE, which is primarily based on Calzetti et al. (2000) extending to short wavelengths (91.2–150 nm) with Leitherer et al. (2002). The allowed range of color excess is E(B-V) = 0–1 (see Table 1).

For galactic dust emission, we use the dl2014 module (Draine et al. 2014). The dust luminosity numerically equals the attenuated luminosity of starlight, because CIGALE strictly follows energy conservation. The Draine et al. (2014) recipe models the dust emission with two components, a diffused emission and a photodissociation region (PDR) emission associated with star formation. The two components share the same mass fraction of polycyclic aromatic hydrocarbon (PAH) compared to total dust, and we allow three different values of it (see Table 1). The major difference between the two components is the radiation field. The diffuse dust is radiated by the minimum radiation parameter (U<sub>min</sub>) set by the user, and the PDR’s radiation field ranges from U<sub>min</sub> to a maximum fixed value of U<sub>max</sub> = 10<sup>7</sup>. We allow U<sub>min</sub> varying from 0.1, 1.0, 10, and 50. The relative strength between the two components is controlled by the γ parameter (the fraction of PDR emission), and we allow it to vary in the range 0.01–0.9.

For the AGN component, we adopt the skirtor2016 module based on a clumpy torus model from Stalevski et al. (2012, 2016). The relative strength between the AGN and galaxy components are set by the frac<sub>AGN</sub> parameter, i.e., AGN fractional luminosity compared to the total over a user-defined rest-frame wavelength range (λ<sub>AGN</sub>). In this work, we allow frac<sub>AGN</sub> to vary from 0 to 0.99 and set λ<sub>AGN</sub> to a mid-IR wavelength range of 3–30 µm that contains the bulk of AGN dust emission (Table 1). We allow all available values of the 9.7 µm optical depth (i.e., 3, 5, 7, 9, and 11) in skirtor2016. The parameter of the AGN viewing angle determines the AGN type. Following Yang et al. (2021b), we
only allow type 2 AGN, because our focus in this work is to
find the obscured AGN untraceable in the UV/optical.
Including type 1 models would cause a significant model
degeneracy versus the stellar emission in the UV/optical (e.g.,
Ni et al. 2021; Zou et al. 2022). Also, type 1 AGN are rare
especially for small-area surveys like our case, and we expect
only ≲1 type 1 AGN using the empirical number density in
Yang et al. (2018). We fix the viewing angle at 70°, a typical
value for type 2 AGN (Yang et al. 2020), because the fit results
are generally insensitive to different type 2 viewing angles
(e.g., Ramos Padilla et al. 2002). Here, by using skir-
tor2016 we assume AGN obscuration is mainly caused by
dust on the torus scale (~1 pc) rather than on the galactic scale
(≥100 pc). This is because galactic-scale dust is likely heated
to a low temperature (≤100 K), and consequently the
reprocessed emission concentrates on far-IR, beyond the MIRI
wavelength coverage. We discuss the effect of possible
galactic-scale obscuration in Section 3.3.

Since most (95%) of our sources do not have spectroscopic
redshifts (spec-z) available, we allow CIGALE to fit a redshift
for each source.<fnref n="36" /> We set the redshift grid from 0.01 to 8.0
(Table 1) with 50 steps evenly spaced in terms of log(1 + z).
We note that changing the maximum allowed redshift (z = 8) to
higher values does not affect the fit results. We evaluate the
best-fit photometric redshift (photo-z) quality using the
available spec-z measurements in Section 2.3.

2.3. CIGALE Results

We run CIGALE with the settings described in Section 2.2.
We adopt the Bayesian output instead of the best-fit output.
The former is probability-weighted values considering all
models available, while the latter is from the minimum-$\chi^2$
model alone. Therefore, the Bayesian output is generally more
robust than the best-fit output. Figure 1 displays three SED-fit
examples from our run and the corresponding image cutouts.

A small subsample of our sources (29 out of 560) have secure
spec-z measurements available from the CANDEL-EGS catalog
(Stefanon et al. 2017). We compare these spec-z values versus our
photo-z estimation in Figure 2 (top). We assess the photo-z quality
by two parameters, i.e., normalized median absolute deviation
($\sigma_{\text{NMAD}}$) and outlier fraction ($\text{frac}_{\text{out}}$). The former is defined as
$\sigma_{\text{NMAD}} = 1.48 \times \text{median}\{|\Delta z - \text{median}(\Delta z)|/(1+z_{\text{spec}})\}$,
where $\Delta z$ is the difference between photo-z and spec-z; the latter
is defined as the fraction of outliers [defined as $|\Delta z|/(1+z_{\text{spec}})>0.15$] among all spec-z sources. We obtain
satisfactory photo-z quality of $\sigma_{\text{NMAD}} = 0.032$ and $\text{frac}_{\text{out}} =
0.00\%$ for our spec-z subsample. We attribute the remarkable zero
outlier fraction to the MIRI coverage of the mid-IR PAH features
and CANDELS UV-to-near-IR coverage of Lyman and/or
Balmer breaks. Indeed, previous simulations have indicated that
MIRI photometry could potentially reduce the number of photo-z
outliers significantly, thanks to the coverage of the PAH features
as a redshift indicator (e.g., Bisigello et al. 2016; Yang et al.
2021b). From Figure 2, the measured photo-z tends to be
systematically lower than the spec-z (median offset ~3%). The
exact cause is not clear, but this minor systematic does not affect
our main results. To assess the dependence of photo-z quality on
source brightness, we also plot the photo-z error (compared to

spec-z) versus the H-band magnitude. It appears that there is no
obvious trend between the error and the magnitude.

We classify each source as a star-forming galaxy, a star-
forming-AGN mixed system, or a pure AGN (hereafter SF,
mixed, and AGN), based on the best-fit frac<sub>AGN</sub> parameter (see
Section 2.2 for its definition). For an object with frac<sub>AGN</sub> ≦ 0.1,
0.1 < frac<sub>AGN</sub> < 0.5, and frac<sub>AGN</sub> ≧ 0.5, we classify it as SF,
mixed, and AGN, respectively. From these criteria, 433, 102,
and 25 objects are classified as SF, mixed, and AGN,
respectively. We note that slightly adjusting these empirical
classification criteria does not change our main results
qualitatively. Figure 3 displays the photo-z distributions for
different object types.

The CEERS field has X-ray coverage from Chandra (Nandra et al. 2015). Among our MIRI sources, there are eight X-ray-detected AGN (selected as $L_X > 10^{42.5}$ erg s⁻¹; e.g., Yang et al. 2018), and one, four, and three of these eight X-ray AGN are classified as SF, mixed, and AGN, respectively, by our SED method above. Therefore, our classifications successfully retrieved $7/8 = 87.5\%$ of the X-ray AGN. The only X-ray AGN misclassified as SF has $\text{frac}_{\text{AGN}} = 0.08 \pm 0.06$, marginally below our threshold for the mixed type ($\text{frac}_{\text{AGN}} = 0.1$).

3. Discussion

In this section, we present the mid-IR SEDs for SF, mixed,
and AGN types (Section 3.1) and discuss the cosmic evolution
history of supermassive BHs (Section 3.2). We give some
caveats for our results in Section 3.3.

3.1. Median SEDs

In Section 2.3, we have classified each source as SF, mixed,
or AGN based on the SED-fit result. We now study the
ensemble SED features for each source type from their median
SEDs in the rest frame.

To estimate the median SED, we first calculate the rest-frame
SED for each individual source by de-redshifting the observed
MIRI-band wavelengths, e.g., the F1500W band for a $z=2$
source is shifted to rest-frame 5 µm. In this procedure, we only
adopt MIRI photometry with S/N > 3 to guarantee good data
quality. We obtain each individual SED by linearly interpolat-
ing the MIRI flux density and normalizing each resulting SED
at rest-frame 3.6 µm. If an SED does not cover 3.6 µm, we
discard this individual SED in the subsequent median-SED
estimation. The numbers of SEDs without 3.6 µm coverage are
263 (SF), 79 (mixed), and 7 (AGN). The value of normalization
wavelength (3.6 µm) is chosen to optimize the AGN sample
size (our main focus in this work) as well as to avoid strong
emission features (PAHs and lines). These individual SEDs are
displayed as the gray curves in Figure 4. Finally, we obtain the
median SED at each wavelength by taking the median value of
all individual SEDs that cover this wavelength. We estimate the
associated 1σ scatter using the 16th–84th percentiles. In this
step, if the number of individual SEDs drops below 5, we do
not calculate the median SED, and this criterion leads to the
shortest and longest wavelengths of the median SED (differing
by source types). Note that we do not divide our sources into
different redshift bins and study the median SEDs for each bin,
mainly because our sample sizes (especially for AGN) are
limited (see Section 3.3). Also, we use the interpolated SEDs
directly instead of the best-fit SEDs from CIGALE, because
choosing the latter would be more strongly model dependent.

<footnote n="36">³⁶ For consistency across the sample, we also fit the redshift for spec-z sources. We note that fixing their redshifts at spec-z does not affect our main conclusions.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left: SED-fit examples from CIGALE. The object in the top/middle/bottom panel is classified as SF/mixed/AGN. Within each panel, the top subpanel displays the data points and the best-fit SED models (i.e., flux density vs. observed wavelength). The purple and red data points indicate the observed and model flux densities, respectively. The black curves represent the total model SEDs. The red and orange curves indicate the galactic dust and AGN components, respectively. The blue/yellow curves represent unattenuated/attenuated stellar. The fitted redshift and frac<sub>AGN</sub> are labeled. The bottom subpanel displays the relative residual in the flux, i.e., (observed-model)/observed. Right: the corresponding 5" × 5" cutouts from HST/F160W and MIRI.</figcaption></figure>

Figure 4 displays the resulting median SED for each type. In the left panel (SF), we also overplot the observed SEDs of two local star-forming galaxies, NGC 5992 and NGC 6090, from Brown et al. (2014), and these two SEDs appear similar to each other. The broad shape of our median SF SED is generally similar to these two local-SF SEDs: all three SEDs are

relatively flat at wavelengths below ~5 µm and rise toward longer wavelengths due to the PDR emission. However, there are also some differences. The PAH emission features appear to be significantly weaker in our median SED. This difference might be mainly caused by the broadband nature of our MIRI photometry (λ/Δλ ~ 5), i.e., we lack the spectral resolution to
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Top: comparison between our photometric redshifts (from CIGALE fits) vs. spectroscopic redshifts (when available). The red solid line indicates the 1:1 relation, while the red dashed lines indicate a 15% uncertainty (beyond which a source would be considered as a photo-z outlier). The photo-z quality, quantified by $\sigma_{NMAD}$ and $frac_{out}$, is marked. We note that the outlier fraction is zero. Bottom: relative difference between photo-z and spec-z as a function of H-band magnitude. Each red star represents the median value for each magnitude bin, while the horizontal error bar indicates the bin coverage. It appears the photo-z quality does not significantly depend on the magnitude.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Photometric redshift distribution of our MIRI sources. The blue, orange, and red bars represent SF, mixed, and AGN types, respectively. The inset panel zooms in on the $z \gtrsim 3$ part of the histogram.</figcaption></figure>

highlight the intrinsic PAH intensities. Also, the photo-z uncertainties inevitably weaken the observed PAH strength in the median SED. Detailed quantification of the PAH strength for our MIRI sources can be performed with detailed SED modeling, but this is beyond the scope of this work. Future MIRI medium-resolution spectrometer (MRS) programs

targeting distant SF galaxies will allow a robust comparison of the PAH strength versus local galaxies (approved MRS programs, e.g., Colina Robledo et al. 2017; Pope et al. 2021). At wavelengths of $\lambda \sim 2 \mu$m, the local-SF SEDs have higher fluxes than our median SED, indicating that the former have stronger stellar emission peaking at near-IR wavelengths. This is understandable, because our sources are MIRI-selected distant SF galaxies and they likely have relatively stronger mid-IR PDR emission compared to the near-IR stellar emission.

The right panel in Figure 4 compares our median AGN SED versus the median SED of hot dust-obscured galaxies (hot DOGs; Fan et al. 2016) and the average SED of Seyfert 2 galaxies (Videla et al. 2013). Hot DOGs are infrared-luminous ($L_{IR} \gtrsim 10^{13} L_{\odot}$) galaxies selected by the Wide-field Infrared Survey Explorer (WISE) typically at $z \gtrsim 2$, and they are considered to be powered by fast-growing heavily obscured AGN (e.g., Eisenhardt et al. 2012; Wu et al. 2012). Seyfert 2 galaxies are type 2 obscured AGN in the local universe, and the average SEDs are based on photometry from the nuclear regions to avoid starlight from the outskirt regions. Videla et al. (2013) published average Seyfert 2 SEDs from two populations: one with near-IR excess (i.e., potential stellar emission) and the other without. We choose the latter to focus on the AGN emission. From Figure 4, our median AGN SED is similar to the hot DOG and Seyfert 2 SEDs. This similarity indicates that the MIRI-selected AGN have comparable physical dust-obscuration structures versus hot DOGs and Seyfert 2s, at least in the ensemble sense. We stress that our MIRI-selected AGN are likely more representative in the BH growth history, because hot DOGs are too rare with a very low number density (a few in 100 deg$^2$; e.g., Assef et al. 2015) and nearby Seyfert 2s are systems in which BHs have largely finished their mass growth (e.g., Yang et al. 2018).

Figure 5 compares the median SEDs from different source types. The SF and AGN SEDs have distinctive shapes. The SF SED is significantly flatter than the AGN SED at $z \lesssim 5$ because the former has relatively strong stellar emission peaking at near-IR wavelengths. In fact, the IRAC-based color-color AGN selection methods are largely based on this difference in SF versus AGN SED shapes (e.g., Stern et al. 2005; Lacy et al. 2007; Donley et al. 2012). As expected, the mixed SED appears to fall between the SF and AGN SEDs, but it tends to be closer to the latter. We publicly release the median SEDs and their associated scatters (as displayed in Figure 5) in the format of electronic tables.

Finally, we point out that the median SEDs are not completely model independent, although they are based on the interpolated photometry instead of the best-fit models. The source classification relies on model fitting, and thus the resulting median SEDs (especially the broad shapes) mimics the corresponding models by construction. The PAH features might be less dependent on the SED-based selections, compared to the broad SED shapes.

## 3.2. Black Hole Growth History

X-ray AGN selections have several advantages of, e.g., high purity and good sensitivity (see Section 1). Thanks to these benefits, X-ray surveys have been widely employed to study the BH cosmic accretion history, often characterized by BHAD, i.e., total BH growth rate per comoving volume. However, X-ray photons can be heavily obscured if the obscuring materials are abundant (Section 1). This leads to the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Median rest-frame SEDs of SF (left), mixed (middle), and AGN (right) types. The red curves indicate the median SEDs. The gray curves indicate the MIRI SEDs (only S/N > 3 MIRI points are being used here for high quality) for individual sources. All these individual SEDs are normalized at rest-frame 3.6 µm. The median SEDs are estimated from these gray curves. In the left panel, we also display two local-SF galaxy SEDs, NGC 5992 and NGC 6090, from Brown et al. (2014). Our median SF SED is similar to the SEDs of these two local-SF galaxies. In the right panel, we also display the median SED for hot dust-obscured galaxies (hot DOG; Fan et al. 2016) and the average SED for the nuclear regions of Seyfert 2 galaxies (Videla et al. 2013). These typical SEDs are similar to our MIRI-selected AGN median SED.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison between different types of rest-frame median SEDs (same as in Figure 4). The blue, orange, and red represent SF, mixed, and AGN types, respectively. The shaded regions indicate the corresponding 1σ scatters, estimated as 16th–84th percentiles. The AGN SED is distinctive from the SF SED. These median SEDs and their associated errors are publicly available with this figure.</figcaption></figure>

(The data used to create this figure are available.)

suspicion that even the deepest X-ray surveys so far might suffer from potential incompleteness issues when characterizing the BHAD evolution (e.g., Hickox & Alexander 2018; Carroll et al. 2021; Yang et al. 2021a). It is thereby crucial to complement the X-ray results with BHAD measurements based on other wavelengths. In this section, we estimate BHAD based on our SED-fit results and compare it with X-ray measurements from the literature.

To estimate BHAD, we first calculate the BH accretion rate (BHAR) for each source classified as AGN or mixed using

$$ 
\begin{aligned}
\text{BHAR} &= \frac{L_{\text{disk}}(1 - \epsilon)}{\epsilon c^2} \\
&= \frac{1.59 L_{\text{disk}}}{10^{46} \text{ erg s}^{-1} \text{ M}_\odot \text{ yr}^{-1}},
\end{aligned}
\quad (1) $$

where c is the speed of light and ϵ is the radiative efficiency. We assume a conventional value of ϵ = 0.1 here (e.g., Brandt & Alexander 2015). L<sub>disk</sub> is the viewing-angle-averaged intrinsic accretion-disk luminosity from our SED fits (i.e., the “accretion_power” quantity in the CIGALE output; see Yang et al. 2018). It is numerically equivalent to the angle-averaged AGN bolometric luminosity due to energy conservation. The L<sub>disk</sub> (BHAR) distributions at different redshifts are displayed in Figure 6.

We then add up all the BHAR in a given redshift bin and divide the total BHAR by the comoving volume in the redshift bin sampled by the survey area (a total of 9.2 arcmin<sup>2</sup> for the four pointings). This yields our BHAD estimate for a given redshift bin, i.e.,

$$ \text{BHAD} = \frac{\sum_i \text{BHAR}_i}{V_c}, \quad (2) $$

where BHAR<sub>i</sub> is the BH accretion rate for the ith source (AGN or mixed type) in the redshift bin; V<sub>c</sub> is the comoving volume for the redshift bin sampled by the MIRI area (9.2 arcmin<sup>2</sup>), and the numerical value is calculated using ASTROPY.COSMOLOGY. Note that we do not correct for incompleteness here. The correction would require detailed knowledge of our MIRI-selection sensitivity as well as AGN IR luminosity function, and that is beyond the scope of this work. We qualitatively discuss the effects of incompleteness in Section 3.3.

The resulting BHAD versus redshift is displayed in Figure 7 as red data points. The redshift bins (same for the red and orange points) are indicated by the horizontal error bars. We estimate the uncertainties with a bootstrap method, using ASTROPY.STATS.BOOTSTRAP (Astropy Collaboration et al. 2018). The bootstrap resample number is 1000, and we adopt the 16th–84th percentile as 1σ confidence limit for each data point. The bootstrap error (δBHAD<sub>boot</sub>) only accounts for statistical fluctuations due to limited sample sizes. To also consider the errors in our SED fits, we propagate the L<sub>disk</sub> Bayesian errors (from the CIGALE output) into the BHAD error
<figure><img src="image_9.png" /><figcaption>**Figure 6.** $L_{\text{disk}}$ distributions for AGN (red, with $\text{frac}_{\text{AGN}} \ge 0.5$) and mixed (blue, with $0.1 < \text{frac}_{\text{AGN}} < 0.5$) objects. Different panels are from different redshift bins as labeled. The vertical dashed lines indicate the typical knee luminosities from Shen et al. (2020). The corresponding BHAR is also marked at the upper axis. $L_{\text{disk}}$ is numerically equivalent to angle-averaged $L_{\text{bol}}$ (see Section 3.2).</figcaption></figure>

($\delta\text{BHAD}_{\text{Bayes}}$, using standard error propagation. Finally, we estimate the final BHAD uncertainty as $\sqrt{(\delta\text{BHAD}_{\text{boot}})^2 + (\delta\text{BHAD}_{\text{Bayes}})^2}$. $\delta\text{BHAD}_{\text{Bayes}}$ is smaller than or comparable to the corresponding $\delta\text{BHAD}_{\text{boot}}$ and their ratios are 0.2–1.0.

In Figure 7, we also show some BHAD curves from the literature. The Yang et al. (2021a) curve is based on the observed SFH of bulge-dominated galaxies at $z = 0.7–1.5$ with the assumption of the bulge BHAR–SFR relation (Yang et al. 2019), and thus this BHAD measurement is an indirect method. The Yang et al. (2021a) BHAD should be considered as a lower limit, as it only accounts for the BH growth within bulge-dominated galaxies. The other literature BHAD curves are based on X-ray data (Ueda et al. 2014; Aird et al. 2015; Vito et al. 2018; Ananna et al. 2019), each including a wide range of surveys (from wide–shallow to narrow–deep) to completely sample the entire X-ray AGN population. The X-ray measurements are compiled by Yang et al. (2021a), assuming the same radiation efficiency ($\epsilon = 0.1$) as our MIRI-based BHAD.

From Figure 7, our MIRI-based BHAD roughly agrees with the X-ray results at low redshifts of $z \lesssim 3$. However, the X-ray BHAD curves drop sharply toward higher redshifts. In contrast,

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Black hole accretion density as a function of redshift. The red data points are derived based on our MIRI-selected AGN; the orange ones are from AGN+mixed sources. The vertical error bars consider both bootstrap and SED-fit (Bayesian) uncertainties; the horizontal error bars indicate the redshift bins. The orange curve representing the inferred BHAR for bulge-dominated galaxies is from Yang et al. (2021a), and is thereby a lower bound on the total accretion rate. The black curves are from X-ray measurements with references marked in the legend (see Yang et al. 2021a for details). At high redshifts of $z \gtrsim 3$, our MIRI-based BHAD values are significantly higher than the corresponding X-ray measurements, but are consistent with the Yang et al. (2021a) lower bound. The blue curve represents SFRD (Madau & Dickinson 2014) downscaled by a factor of 2000 to roughly match the X-ray-based BHAD peak at $z \sim 2$.</figcaption></figure>

our BHAD does not have significant evolution at $z \approx 1–5$. This difference in evolution leads to the fact that our BHAD is much higher than the X-ray values by $\sim 0.5$ dex at $z \gtrsim 3$.

The results above indicate that MIRI detects a great number of AGN that are missed by X-ray observations. Indeed, none of our $z > 3$ AGN/mixed objects are detected by the Chandra X-ray observations (Nandra et al. 2015). This lack of X-ray signals is likely due to heavy obscuration at high redshifts (see Section 1). In the early universe, the obscuring gas is likely more abundant and compact compared to low redshifts, resulting in a high column density that is capable of heavily damping the X-ray emission (e.g., Gilli et al. 2022). On the other hand, heavy AGN obscuration should be associated with high IR luminosities due to strong dust absorption and reemission following the energy-conservation law (Section 1). This explains why MIRI is able to detect the X-ray-missed high-$z$ AGN.

Now, we quantitatively constrain the X-ray obscuration $N_H$ via a stacking analysis. We apply the same stacking method of Yang et al. (2019) to our AGN at $z = 3–5$ (10 sources, all undetected in X-ray). In brief, this method extracts X-ray count rates at the H-band position of each source, and then converts the average count rate to the average $L_X$ (X-ray luminosity in rest-frame 2–10 keV). The details of the stacking method are described in Section 2.4 of Yang et al. (2019). For our high-$z$ AGN sample, the stacked $L_X$ is consistent with zero (likely due to our limited sample size), having a $3\sigma$ upper limit of $2.7 \times 10^{42} \text{ erg s}^{-1}$. On the hand, based on our fitted $L_{\text{disk}}$ (Figure 6), we have estimated an average intrinsic $L_X$ of
5.9 × 10<sup>43</sup> erg s<sup>-1</sup>, assuming an X-ray bolometric correction factor of 22.4 (e.g., Vasudevan & Fabian 2007). From this expected intrinsic L<sub>X</sub> and the observed 3σ upper limit of L<sub>X</sub>, we estimate a damping factor of at least ~22 (5.9 × 10<sup>43</sup> divided by 2.7 × 10<sup>42</sup>) for the AGN X-ray radiation. To reach this level of obscuration, the N<sub>H</sub> should be at least ~2.0 × 10<sup>24</sup> cm<sup>-2</sup>,<sup>37</sup> above the Compton-thick threshold 1.5 × 10<sup>24</sup> cm<sup>-2</sup>. Therefore, we conclude that our MIRI-detected high-z AGN are heavily obscured in X-ray on average. We have also applied the analysis above to the sample of AGN and mixed sources (25 objects) at z = 3–5 and reach a similar result. We note that it is not feasible to directly infer X-ray N<sub>H</sub> from the obscuration estimated from our SED fitting of the UV-to-IR photometry. This is because UV-to-IR radiation is dominantly reprocessed by dust, while X-ray is dominantly reprocessed by gas (e.g., Hickox & Alexander 2018). Actually, the skirtor2016 model we used in SED fitting does not include any gas at all (Stalevski et al. 2012, 2016), and thus we are not able to infer N<sub>H</sub> from our SED modeling.

Figure 7 also displays the redshift evolution of SFRD from Madau & Dickinson (2014). In the plot, the SFRD is numerically downscaled by a factor of 2000 to roughly match the X-ray-based BHAD peak at $z \sim 2$. Qualitatively similar to the BHAD measurements from X-ray data, the SFRD also decreases at $z \gtrsim 2$. This behavior leads to that the SFRD falls below our BHAD at $z \gtrsim 3$, indicating the BHAD/SFRD ratio rises at high redshifts. This result poses a strong question to theorists, as current simulations predict that the BHAD/SFR ratio should drop (or at most remain constant) in the early universe (e.g., Habouzit et al. 2021; Zhang et al. 2023). New theoretical models are perhaps needed for more efficient BH growth at high redshifts.

### 3.3. Caveats

Our work represents one of the first papers using JWST/ MIRI to search for obscured AGN. Currently, we are still in the early stage of JWST operation, and thereby this Letter inevitably has some limitations, as we discuss below.

A major limitation of this work is the sample size, only a total of 560 sources (433 SF, 102 mixed, and 25 AGN), over a wide redshift range of $z \approx 0-5$. This relatively small sample size is due to the survey area we have, i.e., only four pointings with a total area of 9 arcmin². The sample sizes (especially for AGN and mixed types) prevent us from studying the redshift evolution of the median SEDs (Section 3.1).

For the BHAD measurements (Section 3.2), the contributions from rare luminous AGN (such as hot DOGs) are missed, because our MIRI survey area is not sufficiently large enough to capture them. Indeed, most of our sources lie below the knee luminosity (dashed lines in Figure 6). For comparison, the least luminous hot DOGs have $L_{\text{disk}} \sim 10^{47}$ erg s⁻¹ above the knee luminosity (e.g., Fan et al. 2016). In addition, our sample does not include unobscured type 1 AGN, which is also generally luminous, due to our limited survey area (see Section 2.2). Therefore, the intrinsic BHAD might be higher than our estimated values. However, we note that this bias actually strengthens our main conclusion, i.e., the BH growth in the early universe is stronger than we thought based on X-ray detections. On the faint side, our MIRI data could miss low-luminosity AGN especially at high redshifts. For example, all

the sources at $z > 3$ have $L_{\text{disk}} \gtrsim 10^{43}$ erg s⁻¹ in Figure 6. Another AGN population we could miss is those mainly obscured by galactic-scale dust (e.g., Gilli et al. 2022), of which the reprocessed emission is beyond the MIRI wavelength coverage (see Section 2.2). Such biases again go to the direction of strengthening our main conclusion. From the discussion above, the contributions from many AGN could be missed in the MIRI-based BHAD measurements. Therefore, our estimation effectively serves as a lower-limit constraint to the intrinsic complete BHAD.

Another limitation of this work is the lack of secure spec-z for most sources in our sample, especially the high-z AGN/mixed objects that are crucial to our main results. Although the photo-z quality appears high (Section 2.3), the spec-z sources used in this quality assessment are all at low redshifts of $z < 3$ (see Figure 2). To have a rough evaluation of the quality at high redshifts, we compare the photo-z of our $z > 3$ sources versus the CANDELS-EGS photo-z (Stefanon et al. 2017). Only 1 (out of 36) source has a fractional difference >15%, indicating that our selected high-z sample is also supported by the CANDELS-EGS work. The difference between our $z > 3$ BHAD and the X-ray measurements is ~0.5 dex (Figure 7). Therefore, if this large difference were due to photo-z errors, most ($\gtrsim$70%) of our photo-z estimations for our high-z sources would have to be spurious.

Finally, our SED-based method is model dependent. It is possible that our adopted IR models (see Section 2.2) are not suitable in particular for MIRI-selected high-z objects because these sources are at a low mid-IR flux level that has not been well studied. For example, if the 3.3 µm PAH feature is exceptionally strong, the observed MIRI SED of a $z \approx 4-5$ object could be very red and thereby misclassified as an AGN (e.g., Magnelli et al. 2008). This is because the rest-frame 3.3 µm is shifted to F1800W or F2100W (our reddest MIRI bands). To address this potential issue, we reperform our analyses but discard the F1800W and F2100W data. The resulting BHAD is consistent with the value based on all MIRI bands, which is still significantly higher than the X-ray measurements at $z \approx 4-5$. Therefore, the possible model uncertainties of about 3.3 µm PAH are unlikely to affect our main conclusion qualitatively.

## 4. Summary and Future Prospects

In this work, we have performed SED analysis for a MIRI-detected sample in the CEERS survey. Our main results are summarized below.

1. Our analyses are based on the four MIRI red pointings of the CEERS survey (~9 arcmin²). From the CEERS MIRI catalog, we have selected sources with at least two MIRI-band S/N > 3 and compiled other multiwavelength (UV-to-IRAC 4) photometry from the CANDELS catalog. We have a total of 560 MIRI-detected sources.

2. We have performed SED modeling for the 560 sources with CIGALE. In this procedure, we have fitted the redshift and other source properties (e.g., stellar e-folding time and age) simultaneously. We compared our results to the available spec-z of 29 objects in the CANDELS-EGS catalog. We found $\sigma_{\text{NMAD}} = 0.032$ and no outliers with >15% uncertainties, indicating our photo-z estimate is robust. Based on the best-fit frac<sub>AGN</sub> (fractional AGN

<footnote n="37"><sup>37</sup> This is estimated with PIMMS, https://cxc.harvard.edu/toolkit/pimms.jsp.</footnote>
luminosity at rest-frame 3–30 µm), we have classified each source into SF, mixed, and AGN, respectively.

3. We have derived median rest-frame SEDs for the three types using MIRI photometry. The median SF SED has a similar broad shape compared to local-SF galaxies, but the PAH features appear to be weaker, probably due to the broad wavelength coverage of the MIRI filters and/or redshift uncertainties. The median AGN SED is very close to the typical SEDs of hot DOGs and Seyfert 2s, suggesting they are intrinsically the same type of object, i.e., actively accreting but obscured BHs. We publicly release our median SEDs along with this Letter.

4. We have studied cosmic BH accretion history by estimating the BHAD as a function of redshift based on our classified AGN (and mixed) types of sources. The result agrees with the X-ray measurements at $z \lesssim 3$. However, our MIRI-based BHAD becomes significantly higher than the latter toward higher redshifts, and the difference is $\sim 0.5$ dex at $z \gtrsim 3$. We interpret this difference as that MIRI is able to detect many heavily obscured AGN in the early universe.

With the accumulation of total JWST observing time, there will be more and more archival MIRI imaging data suitable for such studies as this work. This will effectively address the potential issues related to the small sample size and also allows studying the redshift evolution of the median SEDs (Section 3.3). Future MIRI/MRS follow-up observations of the imaging-selected SF/AGN targets will also be useful: MRS can reveal the redshift evolution of PAH strength among SF galaxies (Section 3.1); MRS can provide secure spectroscopic redshifts for our high-$z$ AGN candidates, essentially testing our conclusion of the BHAD excess compared to X-ray measurements (Section 3.2).

An alternative approach to probe AGN at high redshifts is mining the existing JWST spectroscopic database (e.g., NIRSpec Multi-Object Spectroscopy and NIRCam grism slitless spectroscopy), many of which target high-$z$ galaxy candidates. Based on the spectra, the presence of AGN can be identified by searching for broad emission lines and/or narrow high-ionization lines (e.g., Cleri et al. 2023; Kocevski et al. 2023). The spectroscopic selections can potentially probe AGN at the highest redshifts. For example, the recent work of Larson et al. (2023) reported a broad-line AGN at $z = 8.679$, the earliest BH discovered so far.

## Acknowledgments

We thank the referee for helpful feedback that improved this work. We acknowledge the hard work of our colleagues in the CEERS collaboration and everyone involved in the JWST mission. G.Y., K.I.C., and E.I. acknowledge funding from the Netherlands Research School for Astronomy (NOVA). K.I.C. and V.K. acknowledge funding from the Dutch Research Council (NWO) through the award of the Vici Grant VI. C.212.036. C.P. thanks Marsha and Ralph Schilling for generous support of this research. T.A.H. is supported by an appointment to the NASA Postdoctoral Program (NPP) at NASA Goddard Space Flight Center, administered by Oak Ridge Associated Universities under contract with NASA. This work acknowledges support from the NASA/ESA/CSA James Webb Space Telescope through the Space Telescope Science Institute, which is operated by the Association of Universities

for Research in Astronomy, Incorporated, under NASA contract NAS5-03127. Support for program No. JWST-ERS01345 was provided through a grant from the STScI under NASA contract NAS5-03127. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific observations analyzed can be accessed via doi:10.17909/agda-2w34. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS526555. Support to MAST for these data is provided by the NASA Office of Space Science via grant NAG57584 and by other grants and contracts.

Software: ASTROPY (Astropy Collaboration et al. 2018), CIGALE (Boquien et al. 2019; Yang et al. 2020, 2022).

## ORCID iDs

G. Yang (杨光) https://orcid.org/0000-0001-8835-7722

K. I. Caputi https://orcid.org/0000-0001-8183-1460

C. Papovich https://orcid.org/0000-0001-7503-8482

P. Arrabal Haro https://orcid.org/0000-0002-7959-8783

M. B. Bagley https://orcid.org/0000-0002-9921-9218

P. Behroozi https://orcid.org/0000-0002-2517-6446

E. F. Bell https://orcid.org/0000-0002-5564-9873

L. Bisigello https://orcid.org/0000-0003-0492-4924

V. Buat https://orcid.org/0000-0003-3441-903X

D. Burgarella https://orcid.org/0000-0002-4193-2539

Y. Cheng https://orcid.org/0000-0001-8551-071X

N. J. Cleri https://orcid.org/0000-0001-7151-009X

R. Davé https://orcid.org/0000-0003-2842-9434

M. Dickinson https://orcid.org/0000-0001-5414-5131

D. Elbaz https://orcid.org/0000-0002-7631-647X

H. C. Ferguson https://orcid.org/0000-0001-7113-2738

S. L. Finkelstein https://orcid.org/0000-0001-8519-1130

N. A. Grogin https://orcid.org/0000-0001-9440-8872

N. P. Hathi https://orcid.org/0000-0001-6145-5090

M. Hirschmann https://orcid.org/0000-0002-3301-3321

B. W. Holwerda https://orcid.org/0000-0002-4884-6756

M. Huertas-Company https://orcid.org/0000-0002-1416-8483

T. A. Hutchison https://orcid.org/0000-0001-6251-4988

E. Iani https://orcid.org/0000-0001-8386-3546

J. S. Kartaltepe https://orcid.org/0000-0001-9187-3605

A. Kirkpatrick https://orcid.org/0000-0002-1306-1545

D. D. Kocevski https://orcid.org/0000-0002-8360-3880

A. M. Koekemoer https://orcid.org/0000-0002-6610-2048

V. Kokorev https://orcid.org/0000-0002-5588-9156

R. L. Larson https://orcid.org/0000-0003-2366-8858

R. A. Lucas https://orcid.org/0000-0003-1581-7825

P. G. Pérez-González https://orcid.org/0000-0003-4528-5639

P. Rinaldi https://orcid.org/0000-0002-5104-8245

L. Shen https://orcid.org/0000-0001-9495-7759

J. R. Trump https://orcid.org/0000-0002-1410-0470

A. de la Vega https://orcid.org/0000-0002-6219-5558

L. Y. A. Yung https://orcid.org/0000-0003-3466-035X

J. A. Zavala https://orcid.org/0000-0002-7051-1100

## References

Aird, J., Coil, A. L., & Georgakakis, A. 2019, MNRAS, 484, 4360

Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892

Aird, J., Coil, A. L., & Kocevski, D. D. 2022, MNRAS, 515, 4860
Alonso-Herrero, A., Pérez-González, P. G., Alexander, D. M., et al. 2006, ApJ, 640, 167

Ananna, T. T., Treister, E., Urry, C. M., et al. 2019, ApJ, 871, 240

Assef, R. J., Eisenhardt, P. R. M., Stern, D., et al. 2015, ApJ, 804, 27

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Bisigello, L., Caputi, K. I., Colina, L., et al. 2016, ApJS, 227, 19

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Brandt, W. N., & Alexander, D. M. 2015, A&ARv, 23, 1

Brandt, W. N., & Yang, G. 2022, in Handbook of X-ray and Gamma-ray Astrophysics, ed. C. Bambi & A. Santangelo (Singapore: Springer), doi:10.1007/978-981-16-4544-0_130-1

Brown, M. J. I., Moustakas, J., Smith, J. D. T., et al. 2014, ApJS, 212, 18

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Caputi, K. I. 2013, ApJ, 768, 103

Carroll, C. M., Hickox, R. C., Masini, A., et al. 2021, ApJ, 908, 185

Chabrier, G. 2003, ApJL, 586, L133

Cleri, N. J., Yang, G., Papovich, C., et al. 2023, ApJ, 948, 112

Colina Robledo, L., Ferruit, P., Alvarez-Marquez, J., et al. 2017, NIRSpec and MIRI IFS of SMGs, JWST Proposal Cycle 1, ID, #1264

Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Eisenhardt, P. R. M., Wu, J., Tsai, C.-W., et al. 2012, ApJ, 755, 173

Fan, L., Han, Y., Nikutta, R., Drouart, G., & Knudsen, K. K. 2016, ApJ, 823, 107

Finkelstein, S. L., Dickinson, M., Ferguson, H. C., et al. 2017, The Cosmic Evolution Early Release Science (CEERS) Survey, JWST Proposal ID 1345. Cycle 0 Early Release Science

Gilli, R., Norman, C., Calura, F., et al. 2022, A&A, 666, A17

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Habouzit, M., Li, Y., Somerville, R. S., et al. 2021, MNRAS, 503, 1940

Hickox, R. C., & Alexander, D. M. 2018, ARA&A, 56, 625

Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055

Kirkpatrick, A., Alberts, S., Pope, A., et al. 2017, ApJ, 849, 111

Kocevski, D. D., Onoue, M., Inayoshi, K., et al. 2023, arXiv:2302.00012

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Lacy, M., Petric, A. O., Sajina, A., et al. 2007, AJ, 133, 186

Larson, R. L., Finkelstein, S. L., Kocevski, D. D., et al. 2023, arXiv:2303.08918

Leitherer, C., Li, I. H., Calzetti, D., & Heckman, T. M. 2002, ApJS, 140, 303

Lyu, J., Alberts, S., Rieke, G. H., & Rujopakarn, W. 2022, ApJ, 941, 191

Lyu, J., & Rieke, G. 2022, Univ, 8, 304

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Magnelli, B., Chary, R. R., Pope, A., et al. 2008, ApJ, 681, 258

Mendez, A. J., Coil, A. L., Aird, J., et al. 2013, ApJ, 770, 40

Merlin, E., Fontana, A., Ferguson, H. C., et al. 2015, A&A, 582, A15

Merloni, A., & Heinz, S. 2008, MNRAS, 388, 1011

Nandra, K., Laird, E. S., Aird, J. A., et al. 2015, ApJS, 220, 10

Ni, Q., Brandt, W. N., Yang, G., et al. 2021, MNRAS, 500, 4989

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&ARv, 25, 2

Papovich, C., Cole, J., Yang, G., et al. 2023, ApJL, 949, L18

Pope, A., Sajina, A., Yan, L., et al. 2021, Halfway to the Peak: A Bridge Program to Map Coeval Star Formation and Supermassive Black Hole Growth, JWST Proposal Cycle 1, ID. #1762

Ramos Padilla, A. F., Wang, L., Malek, K., Efstathiou, A., & Yang(杨光), G. 2002, AJ, 510, 687

Richards, G. T., Fan, X., Newberg, H. J., et al. 2002, AJ, 123, 2945

Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 495, 3252

Stalevski, M., Fritz, J., Baes, M., Nakos, T., & Popović, L. Č. 2012, MNRAS, 420, 2756

Stalevski, M., Ricci, C., Ueda, Y., et al. 2016, MNRAS, 458, 2288

Stefanon, M., Yan, H., Mobasher, B., et al. 2017, ApJS, 229, 32

Stern, D., Eisenhardt, P., Gorjian, V., et al. 2005, ApJ, 631, 163

Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104

Vasudevan, R. V., & Fabian, A. C. 2007, MNRAS, 381, 1235

Videla, L., Lira, P., Andrews, H., et al. 2013, ApJS, 204, 23

Villa-Vélez, J. A., Buat, V., Theulé, P., Boquien, M., & Burgarella, D. 2021, A&A, 654, A153

Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS, 473, 2378

Wu, J., Tsai, C.-W., Sayers, J., et al. 2012, ApJ, 756, 96

Yang, G., Boquien, M., Brandt, W. N., et al. 2022, ApJ, 927, 192

Yang, G., Boquien, M., Buat, V., et al. 2020, MNRAS, 491, 740

Yang, G., Brandt, W. N., Alexander, D. M., et al. 2019, MNRAS, 485, 3721

Yang, G., Brandt, W. N., Vito, F., et al. 2018, MNRAS, 475, 1887

Yang, G., Estrada-Carpenter, V., Papovich, C., et al. 2021a, ApJ, 921, 170

Yang, G., Papovich, C., Bagley, M. B., et al. 2021b, ApJ, 908, 144

Zhang, H., Behroozi, P., Volonteri, M., et al. 2023, MNRAS, 518, 2123

Zou, F., Brandt, W. N., Chen, C.-T., et al. 2022, ApJS, 262, 15