CEERS Key Paper. IV. A Triality in the Nature of HST-dark Galaxies

Pablo G. Pérez-González, Guillermo Barro, Marianna Annunziatella, Luca Costantin, Ángela García-Argumánez, Elizabeth J. Mcgrath, Rosa M. Mérida, Jorge A. Zavala, Pablo Arrabal Haro, Micaela B. Bagley, et al.

► To cite this version:

Pablo G. Pérez-González, Guillermo Barro, Marianna Annunziatella, Luca Costantin, Ángela García-Argumánez, et al. CEERS Key Paper. IV. A Triality in the Nature of HST-dark Galaxies. The Astrophysical journal letters, 2023, 946, 10.3847/2041-8213/acb3a5. insu-04479097

HAL Id: insu-04479097

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

Submitted on 4 Mar 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
CEERS Key Paper. IV. A Triality in the Nature of HST-dark Galaxies

Pablo G. Pérez-González, Guillermo Barro, Marianna Annunziatella, Luca Costantin, Ángela García-Argumánez, Elizabeth J. Mcgrath, Rosa M. Mérida, Jorge A. Zavala, Pablo Arrabal Haro, Micaela B. Bagley, et al.

► To cite this version:

Pablo G. Pérez-González, Guillermo Barro, Marianna Annunziatella, Luca Costantin, Ángela García-Argumánez, et al. CEERS Key Paper. IV. A Triality in the Nature of HST-dark Galaxies. The Astrophysical Journal Letters, 2023, 946, 10.3847/2041-8213/acb3a5. insu-04479097

HAL Id: insu-04479097

https://insu.hal.science/insu-04479097v1

Submitted on 4 Mar 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
<figure><img src="image_2.png" /></figure>

CEERS Key Paper. IV. A Triality in the Nature of HST-dark Galaxies

Pablo G. Pérez-González<fnref n="1" /> ID, Guillermo Barro<fnref n="2" /> ID, Marianna Annunziatella<fnref n="1" /> ID, Luca Costantin<fnref n="1" /> ID, Ángela García-Argumánez<fnref n="3" />,<fnref n="4" /> ID,
Elizabeth J. McGrath<fnref n="5" /> ID, Rosa M. Mérida<fnref n="1" /> ID, Jorge A. Zavala<fnref n="6" /> ID, Pablo Arrabal Haro<fnref n="7" /> ID, Micaela B. Bagley<fnref n="8" /> ID,

Bren E. Backhaus<fnref n="9" /> ID, Peter Behroozi<fnref n="10" />,<fnref n="11" /> ID, Eric F. Bell<fnref n="12" /> ID, Laura Bisigello<fnref n="13" />,<fnref n="14" /> ID, Véronique Buat<fnref n="15" /> ID, Antonello Calabro<fnref n="16" /> ID,
Caitlin M. Casey<fnref n="8" /> ID, Nikko J. Cleri<fnref n="17" />,<fnref n="18" /> ID, Rosemary T. Coogan<fnref n="19" /> ID, M. C. Cooper<fnref n="20" /> ID, Asantha R. Cooray<fnref n="20" /> ID,

Avishai Dekel<fnref n="21" /> ID, Mark Dickinson<fnref n="7" /> ID, David Elbaz<fnref n="19" /> ID, Henry C. Ferguson<fnref n="22" /> ID, Steven L. Finkelstein<fnref n="8" /> ID, Adriano Fontana<fnref n="23" /> ID,
Maximilien Franco<fnref n="8" /> ID, Jonathan P. Gardner<fnref n="24" /> ID, Mauro Giavalisco<fnref n="25" /> ID, Carlos Gómez-Guijarro<fnref n="19" /> ID, Andrea Grazian<fnref n="14" /> ID,
Norman A. Grogin<fnref n="22" /> ID, Yuchen Guo<fnref n="8" /> ID, Marc Huertas-Company<fnref n="26" />,<fnref n="27" />,<fnref n="28" /> ID, Shardha Jogee<fnref n="8" /> ID, Jeyhan S. Kartaltepe<fnref n="29" /> ID,
Lisa J. Kewley<fnref n="30" /> ID, Allison Kirkpatrick<fnref n="31" /> ID, Dale D. Kocevski<fnref n="5" /> ID, Anton M. Koekemoer<fnref n="32" /> ID, Arianna S. Long<fnref n="8" />,<fnref n="44" /> ID,
Jennifer M. Lotz<fnref n="33" /> ID, Ray A. Lucas<fnref n="32" /> ID, Casey Papovich<fnref n="17" />,<fnref n="18" /> ID, Nor Pirzkal<fnref n="34" /> ID, Swara Ravindranath<fnref n="32" /> ID,

Rachel S. Somerville<fnref n="35" /> ID, Sandro Tacchella<fnref n="36" />,<fnref n="37" /> ID, Jonathan R. Trump<fnref n="9" /> ID, Weichen Wang<fnref n="38" /> ID, Stephen M. Wilkins³⁹,<fnref n="40" /> ID,
Stijn Wuyts<fnref n="41" /> ID, Guang Yang<fnref n="42" />,<fnref n="43" /> ID, and L. Y. Aaron Yung<fnref n="24" />,<fnref n="45" /> ID

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

<footnote n="1">¹Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain; pgperez@cab.inta-csic.es</footnote>

<footnote n="2">²Department of Physics, University of the Pacific, Stockton, CA 90340 USA</footnote>

<footnote n="3">³Departamento de Física de la Tierra y Astrofísica, Facultad de CC Físicas, Universidad Complutense de Madrid, E-28040, Madrid, Spain</footnote>

<footnote n="4">⁴Instituto de Física de Partículas y del Cosmos IPARCOS, Facultad de CC Físicas, Universidad Complutense de Madrid, E-28040 Madrid, Spain</footnote>

<footnote n="5">⁵Department of Physics and Astronomy, Colby College, Waterville, ME 04901, USA</footnote>

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

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

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

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

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

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

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

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

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

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

<footnote n="16">¹⁶Osservatorio Astronomico di Roma, via Frascati 33, Monte Porzio Catone, Italy</footnote>

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

<footnote n="18">¹⁸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="19">¹⁹Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France</footnote>

<footnote n="20">²⁰Department of Physics & Astronomy, University of California, Irvine, 4129 Reines Hall, Irvine, CA 92697, USA</footnote>

<footnote n="21">²¹Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel</footnote>

<footnote n="22">²²Space Telescope Science Institute, Baltimore, MD, USA</footnote>

<footnote n="23">²³INAF-Osservatorio Astronomico di Roma, via di Frascati 33, I-00078 Monte Porzio Catone, Italy</footnote>

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

<footnote n="25">²⁵University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA</footnote>

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

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

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

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

<footnote n="30">³⁰Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="31">³¹Department of Physics and Astronomy, University of Kansas, Lawrence, KS 66045, USA</footnote>

<footnote n="32">³²Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="33">³³Gemini Observatory/NSF's National Optical-IR Astronomy Research Laboratory, 950 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="34">³⁴ESA/AURA Space Telescope Science Institute, USA</footnote>

<footnote n="35">³⁵Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY, 10010, USA</footnote>

<footnote n="36">³⁶Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK</footnote>

<footnote n="37">³⁷Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK</footnote>

<footnote n="38">³⁸Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA
Astronomy Centre, University of Sussex, Falmer, Brighton BN1 9QH, UK</footnote>

<footnote n="40">⁴⁰Institute of Space Sciences and Astronomy, University of Malta, Msida MSD 2080, Malta</footnote>

<footnote n="41">⁴¹Department of Physics, University of Bath, Claverton Down, Bath BA2 7AY, UK</footnote>

<footnote n="42">⁴²Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands</footnote>

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

<footnote>Received 2022 October 31; revised 2022 December 20; accepted 2023 January 1; published 2023 March 31</footnote>

<footnote n="44">⁴⁴NASA Hubble Fellow.</footnote>

<footnote n="45">⁴⁵ NASA Postdoctoral Fellow.</footnote>
<figure><img src="image_2.png" /></figure>

CEERS Key Paper. IV. A Triality in the Nature of HST-dark Galaxies

Pablo G. Pérez-González<fnref n="1" />, Guillermo Barro<fnref n="2" />, Marianna Annunziatella<fnref n="1" />, Luca Costantin<fnref n="1" />, Ángela García-Argumánez<fnref n="3" />,<fnref n="4" />, Elizabeth J. McGrath<fnref n="5" />, Rosa M. Mérida<fnref n="1" />, Jorge A. Zavala<fnref n="6" />, Pablo Arrabal Haro<fnref n="7" />, Micaela B. Bagley<fnref n="8" />,

Bren E. Backhaus<fnref n="9" />, Peter Behroozi<fnref n="10" />,<fnref n="11" />, Eric F. Bell<fnref n="12" />, Laura Bisigello<fnref n="13" />,<fnref n="14" />, Véronique Buat<fnref n="15" />, Antonello Calabro<fnref n="16" />, Caitlin M. Casey<fnref n="8" />, Nikko J. Cleri<fnref n="17" />,<fnref n="18" />, Rosemary T. Coogan<fnref n="19" />, M. C. Cooper<fnref n="20" />, Asantha R. Cooray<fnref n="20" />,

Avishai Dekel<fnref n="21" />, Mark Dickinson<fnref n="7" />, David Elbaz<fnref n="19" />, Henry C. Ferguson<fnref n="22" />, Steven L. Finkelstein<fnref n="8" />, Adriano Fontana<fnref n="23" />, Maximilien Franco<fnref n="8" />, Jonathan P. Gardner<fnref n="24" />, Mauro Giavalisco<fnref n="25" />, Carlos Gómez-Guijarro<fnref n="19" />, Andrea Grazian<fnref n="14" />, Norman A. Grogin<fnref n="22" />, Yuchen Guo<fnref n="8" />, Marc Huertas-Company<fnref n="26" />,<fnref n="27" />,<fnref n="28" />, Shardha Jogee<fnref n="8" />, Jeyhan S. Kartaltepe<fnref n="29" />, Lisa J. Kewley<fnref n="30" />, Allison Kirkpatrick<fnref n="31" />, Dale D. Kocevski<fnref n="5" />, Anton M. Koekemoer<fnref n="32" />, Arianna S. Long<fnref n="8" />,<fnref n="44" />, Jennifer M. Lotz<fnref n="33" />, Ray A. Lucas<fnref n="32" />, Casey Papovich<fnref n="17" />,<fnref n="18" />, Nor Pirzkal<fnref n="34" />, Swara Ravindranath<fnref n="32" />,

Rachel S. Somerville<fnref n="35" />, Sandro Tacchella<fnref n="36" />,<fnref n="37" />, Jonathan R. Trump<fnref n="9" />, Weichen Wang<fnref n="38" />, Stephen M. Wilkins³⁹,<fnref n="40" />, Stijn Wuyts<fnref n="41" />, Guang Yang<fnref n="42" />,<fnref n="43" />, and L. Y. Aaron Yung<fnref n="24" />,⁵⁰

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

<footnote n="1">¹Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain; pgperez@cab.inta-csic.es</footnote>

<footnote n="2">²Department of Physics, University of the Pacific, Stockton, CA 90340 USA</footnote>

<footnote n="3">³Departamento de Física de la Tierra y Astrofísica, Facultad de CC Físicas, Universidad Complutense de Madrid, E-28040, Madrid, Spain</footnote>

<footnote n="4">⁴Instituto de Física de Partículas y del Cosmos IPARCOS, Facultad de CC Físicas, Universidad Complutense de Madrid, E-28040 Madrid, Spain</footnote>

<footnote n="5">⁵Department of Physics and Astronomy, Colby College, Waterville, ME 04901, USA</footnote>

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

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

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

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

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

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

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

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

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

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

<footnote n="16">¹⁶Osservatorio Astronomico di Roma, via Frascati 33, Monte Porzio Catone, Italy</footnote>

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

<footnote n="18">¹⁸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="19">¹⁹Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France</footnote>

<footnote n="20">²⁰Department of Physics & Astronomy, University of California, Irvine, 4129 Reines Hall, Irvine, CA 92697, USA</footnote>

<footnote n="21">²¹Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 91904, Israel</footnote>

<footnote n="22">²²Space Telescope Science Institute, Baltimore, MD, USA</footnote>

<footnote n="23">²³INAF-Osservatorio Astronomico di Roma, via di Frascati 33, I-00078 Monte Porzio Catone, Italy</footnote>

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

<footnote n="25">²⁵University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003-9305, USA</footnote>

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

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

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

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

<footnote n="30">³⁰Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="31">³¹Department of Physics and Astronomy, University of Kansas, Lawrence, KS 66045, USA</footnote>

<footnote n="32">³²Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="33">³³Gemini Observatory/NSF's National Optical-IR Astronomy Research Laboratory, 950 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="34">³⁴ESA/AURA Space Telescope Science Institute, USA</footnote>

<footnote n="35">³⁵Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY, 10010, USA</footnote>

<footnote n="36">³⁶Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK</footnote>

<footnote n="37">³⁷Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK</footnote>

<footnote n="38">³⁸Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA
Astronomy Centre, University of Sussex, Falmer, Brighton BN1 9QH, UK</footnote>

<footnote n="40">⁴⁰Institute of Space Sciences and Astronomy, University of Malta, Msida MSD 2080, Malta</footnote>

<footnote n="41">⁴¹Department of Physics, University of Bath, Claverton Down, Bath BA2 7AY, UK</footnote>

<footnote n="42">⁴²Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands</footnote>

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

<footnote>Received 2022 October 31; revised 2022 December 20; accepted 2023 January 1; published 2023 March 31</footnote>

<footnote n="44">⁴⁴NASA Hubble Fellow.</footnote>

<footnote n="45">⁴⁵ NASA Postdoctoral Fellow.</footnote>
Abstract

The new capabilities that JWST offers in the near- and mid-infrared (IR) are used to investigate in unprecedented detail the nature of optical/near-IR-faint, mid-IR-bright sources, with HST-dark galaxies among them. We gather JWST data from the CEERS survey in the Extended Groth Strip, jointly with HST data, and analyze spatially resolved optical-to-mid-IR spectral energy distributions to estimate photometric redshifts in two dimensions and stellar population properties on a pixel-by-pixel basis for red galaxies detected by NIRCam. We select 138 galaxies with $F150W - F356W > 1.5$ mag and $F356W < 27.5$ mag. The nature of these sources is threefold: (1) 71% are dusty star-forming galaxies (SFGs) at $2 < z < 6$ with $9 < \log M_*/M_{\odot} < 11$ and a variety of specific SFRs ($<1$ to $>100$ Gyr$^{-1}$); (2) 18% are quiescent/dormant (i.e., subject to reignition/rejuvenation) galaxies (QGs) at $3 < z < 5$, with $\log M_*/M_{\odot} \sim 10$ and poststarburst mass-weighted ages (0.5-1.0 Gyr); and (3) 11% are strong young starbursts with indications of high equivalent width emission lines (typically, [O III]+Hβ) at $6 < z < 7$ (XELG-z6) and $\log M_*/M_{\odot} \sim 9.5$. The sample is dominated by disk-like galaxies with remarkable compactness for XELG-z6 (effective radii smaller than 0.4 kpc). Large attenuations in SFGs, $2 < A(V) < 5$ mag, are found within 1.5 times the effective radius, approximately 2 kpc, while QGs present $A(V) \sim 0.2$ mag. Our SED-fitting technique reproduces the expected dust emission luminosities of IR-bright and submillimeter galaxies. This study implies high levels of star formation activity between $z \sim 20$ and $z \sim 10$, where virtually 100% of our galaxies had already formed $10^8 M_{\odot}$, 60% had assembled $10^9 M_{\odot}$, and 10% up to $10^{10} M_{\odot}$ (in situ or ex situ).

*Unified Astronomy Thesaurus concepts:* Galaxy formation (595); Galaxy evolution (594); High-redshift galaxies (734); Stellar populations (1622); Broad band photometry (184); Galaxy ages (576); James Webb Space Telescope (2291); Galaxy quenching (2040); Galaxy processes (614); Quenched galaxies (2016); Starburst galaxies (1570); Galaxy photometry (611)

*Supporting material:* machine-readable table

1. Introduction

Every time a new space observatory is launched with
improved capabilities at longer and longer wavelengths, the
universe reveals the existence of surprisingly massive and/or
highly metal-enriched and/or dusty galaxies at higher and
higher redshifts, where some of us naïve astronomers are
surprised to discover such evolved systems. This was the case
for IRAS (e.g., Sanders & Mirabel 1996), ISO (e.g., Chary &
Elbaz 2001), Spitzer (e.g., Pérez-González et al. 2005; Reddy
et al. 2006), and Herschel (e.g., Elbaz et al. 2011; Gruppioni
et al. 2013), but also from ground infrared observatories (e.g.,
Elston et al. 1988; Franx et al. 2003; Daddi et al. 2004;
Chapman et al. 2005; Dannerbauer et al. 2008).

One of the last episodes of this recursive story happened
with the comparison of data coming from the most powerful
mid-infrared (mid-IR) space telescope existing before 2022,
Spitzer, and the most advanced optical space telescope to date,
Hubble. The combination of their capabilities revealed a
population of red galaxies, nearly undetected or even escaping
the sensitivity limits in the near-infrared (near-IR) provided by
the WFC3 instrument on board HST (and even dimmer in the
optical), but relatively bright in the mid-IR, well detected by
the IRAC camera on board Spitzer (Huang et al. 2011; Caputi
et al. 2012; Wang et al. 2016; Alcalde Pampliega et al. 2019;
Sun et al. 2021; Kokorev et al. 2022).

Among the population of mid-IR-bright, near-IR/optically faint galaxies, some are missed by HST due to depth limitations, and they are thus called HST-dark sources, although the name is misleading due to the variety of depths of HST/WFC3 observations on the sky. The nature of most mid-IR-selected systems, and especially HST-dark sources, is still uncertain. This is mainly because their stellar emission could only be detected robustly by IRAC (and their dust content by ALMA, NOEMA, or SCUBA; see references below), their spectral energy distributions are thus poorly constrained with just a few flux data points, and the vast majority of them are too faint for spectrographs (operating at

any wavelength from the UV to the submillimeter range; but
see Swinbank et al. 2012; Williams et al. 2019, and Zhou et al.
2020) to analyze their nature in detail.

However, the discovery of these intrinsically red sources
directly points to an elusive high-redshift (probably $z > 2-3$)
population of either massive galaxies experiencing a very dusty
star-forming event or quiescently evolving from very early
cosmic times. This is supported by the large fraction of
detections by Spitzer/MIPS, Herschel/PACS, and SPIRE
(Huang et al. 2011; Caputi et al. 2012; Alcalde Pampliega
et al. 2019), and at submillimeter and radio wavelengths
(Simpson et al. 2014; Franco et al. 2018; Yamaguchi et
2019; Williams et al. 2019; Wang et al. 2021; Talia et al. 2021;
Manning et al. 2022; Xiao et al. 2022). In addition, a minor
fraction of HST-dark galaxies could be linked to emission-line
galaxies presenting high equivalent widths and thus being
prominent in the mid-IR when the lines enter the IRAC bands
(see Figure 12 in Alcalde Pampliega et al. 2019 and references
therein).

The cosmological importance of HST-dark galaxies to
improve our knowledge about galaxy evolution is not
negligible. The existence of a numerous population of massive
galaxies already in place or even evolving passively at $z > 2-3$,
and even up to $z \sim 6$ (and beyond; see Barrufet et al. 2022;
Endsley et al. 2022; Labbe et al. 2022; Nelson et al. 2022;
Tacchella et al. 2022a), is very difficult to reproduce by state-
of-the-art galaxy evolution models (see discussion in Alcalde
Pampliega et al. 2019 and the recent discussion on JWST
results in Lovell et al. 2023), which also suffer from their
limitations in area to understand the nature of these samples
presenting relatively small volume densities.

In this paper, we want to benefit from the huge jump offered
by the brand-new JWST and its near- and mid-infrared unique
capabilities to understand the nature of this interesting
population of galaxies—HST-dark galaxies among them,
near-infrared-faint and mid-infrared sources in general—whose
nature and relevance are still debatable.
Abstract

The new capabilities that JWST offers in the near- and mid-infrared (IR) are used to investigate in unprecedented detail the nature of optical/near-IR-faint, mid-IR-bright sources, with HST-dark galaxies among them. We gather JWST data from the CEERS survey in the Extended Groth Strip, jointly with HST data, and analyze spatially resolved optical-to-mid-IR spectral energy distributions to estimate photometric redshifts in two dimensions and stellar population properties on a pixel-by-pixel basis for red galaxies detected by NIRCam. We select 138 galaxies with $F150W - F356W > 1.5$ mag and $F356W < 27.5$ mag. The nature of these sources is threefold: (1) 71% are dusty star-forming galaxies (SFGs) at $2 < z < 6$ with $9 < \log M_*/M_{\odot} < 11$ and a variety of specific SFRs ($<1$ to $>100$ Gyr$^{-1}$); (2) 18% are quiescent/dormant (i.e., subject to reignition/rejuvenation) galaxies (QGs) at $3 < z < 5$, with $\log M_*/M_{\odot} \sim 10$ and poststarburst mass-weighted ages (0.5-1.0 Gyr); and (3) 11% are strong young starbursts with indications of high equivalent width emission lines (typically, [O III]+Hβ) at $6 < z < 7$ (XELG-z6) and $\log M_*/M_{\odot} \sim 9.5$. The sample is dominated by disk-like galaxies with remarkable compactness for XELG-z6 (effective radii smaller than 0.4 kpc). Large attenuations in SFGs, $2 < A(V) < 5$ mag, are found within 1.5 times the effective radius, approximately 2 kpc, while QGs present $A(V) \sim 0.2$ mag. Our SED-fitting technique reproduces the expected dust emission luminosities of IR-bright and submillimeter galaxies. This study implies high levels of star formation activity between $z \sim 20$ and $z \sim 10$, where virtually 100% of our galaxies had already formed $10^8 M_{\odot}$, 60% had assembled $10^9 M_{\odot}$, and 10% up to $10^{10} M_{\odot}$ (in situ or ex situ).

*Unified Astronomy Thesaurus concepts:* Galaxy formation (595); Galaxy evolution (594); High-redshift galaxies (734); Stellar populations (1622); Broad band photometry (184); Galaxy ages (576); James Webb Space Telescope (2291); Galaxy quenching (2040); Galaxy processes (614); Quenched galaxies (2016); Starburst galaxies (1570); Galaxy photometry (611)

*Supporting material:* machine-readable table

1. Introduction

Every time a new space observatory is launched with
improved capabilities at longer and longer wavelengths, the
universe reveals the existence of surprisingly massive and/or
highly metal-enriched and/or dusty galaxies at higher and
higher redshifts, where some of us naïve astronomers are
surprised to discover such evolved systems. This was the case
for IRAS (e.g., Sanders & Mirabel 1996), ISO (e.g., Chary &
Elbaz 2001), Spitzer (e.g., Pérez-González et al. 2005; Reddy
et al. 2006), and Herschel (e.g., Elbaz et al. 2011; Gruppioni
et al. 2013), but also from ground infrared observatories (e.g.,
Elston et al. 1988; Franx et al. 2003; Daddi et al. 2004;
Chapman et al. 2005; Dannerbauer et al. 2008).

One of the last episodes of this recursive story happened
with the comparison of data coming from the most powerful
mid-infrared (mid-IR) space telescope existing before 2022,
Spitzer, and the most advanced optical space telescope to date,
Hubble. The combination of their capabilities revealed a
population of red galaxies, nearly undetected or even escaping
the sensitivity limits in the near-infrared (near-IR) provided by
the WFC3 instrument on board HST (and even dimmer in the
optical), but relatively bright in the mid-IR, well detected by
the IRAC camera on board Spitzer (Huang et al. 2011; Caputi
et al. 2012; Wang et al. 2016; Alcalde Pampliega et al. 2019;
Sun et al. 2021; Kokorev et al. 2022).

Among the population of mid-IR-bright, near-IR/optically faint galaxies, some are missed by HST due to depth limitations, and they are thus called HST-dark sources, although the name is misleading due to the variety of depths of HST/WFC3 observations on the sky. The nature of most mid-IR-selected systems, and especially HST-dark sources, is still uncertain. This is mainly because their stellar emission could only be detected robustly by IRAC (and their dust content by ALMA, NOEMA, or SCUBA; see references below), their spectral energy distributions are thus poorly constrained with just a few flux data points, and the vast majority of them are too faint for spectrographs (operating at

any wavelength from the UV to the submillimeter range; but
see Swinbank et al. 2012; Williams et al. 2019, and Zhou et al.
2020) to analyze their nature in detail.

However, the discovery of these intrinsically red sources
directly points to an elusive high-redshift (probably $z > 2-3$)
population of either massive galaxies experiencing a very dusty
star-forming event or quiescently evolving from very early
cosmic times. This is supported by the large fraction of
detections by Spitzer/MIPS, Herschel/PACS, and SPIRE
(Huang et al. 2011; Caputi et al. 2012; Alcalde Pampliega
et al. 2019), and at submillimeter and radio wavelengths
(Simpson et al. 2014; Franco et al. 2018; Yamaguchi et
2019; Williams et al. 2019; Wang et al. 2021; Talia et al. 2021;
Manning et al. 2022; Xiao et al. 2022). In addition, a minor
fraction of HST-dark galaxies could be linked to emission-line
galaxies presenting high equivalent widths and thus being
prominent in the mid-IR when the lines enter the IRAC bands
(see Figure 12 in Alcalde Pampliega et al. 2019 and references
therein).

The cosmological importance of HST-dark galaxies to
improve our knowledge about galaxy evolution is not
negligible. The existence of a numerous population of massive
galaxies already in place or even evolving passively at $z > 2-3$,
and even up to $z \sim 6$ (and beyond; see Barrufet et al. 2022;
Endsley et al. 2022; Labbe et al. 2022; Nelson et al. 2022;
Tacchella et al. 2022a), is very difficult to reproduce by state-
of-the-art galaxy evolution models (see discussion in Alcalde
Pampliega et al. 2019 and the recent discussion on JWST
results in Lovell et al. 2023), which also suffer from their
limitations in area to understand the nature of these samples
presenting relatively small volume densities.

In this paper, we want to benefit from the huge jump offered
by the brand-new JWST and its near- and mid-infrared unique
capabilities to understand the nature of this interesting
population of galaxies—HST-dark galaxies among them,
near-infrared-faint and mid-infrared sources in general—whose
nature and relevance are still debatable.
Indeed, with the advent of JWST, new mid-IR-bright, near-IR-faint galaxies have been identified and started to be characterized (Barrufet et al. 2022; Carnall et al. 2023; Labbe et al. 2022; Nelson et al. 2022). These works have confirmed the high-redshift nature of HST-dark galaxies, their relevance for compiling a complete stellar-mass census in the first 1–2 Gyr of the lifetime of the universe, and their challenging numbers and properties for galaxy evolution models. Disagreement with theoretical predictions has also been (preliminary) found with the large numbers of very high-redshift ($z > 10$) galaxy candidates being found in relatively shallow JWST/NIRCam data (e.g., Adams et al. 2023; Castellano et al. 2022; Donnan et al. 2023; Finkelstein et al. 2022; Harikane et al. 2022; Naidu et al. 2022; Ono et al. 2022; Rodighiero et al. 2023), implying high levels of star formation activity in the very early universe (even though these candidates must still be confirmed and separated from contaminants such as dusty galaxies at $2 < z < 5$, see Zavala et al. 2023).

In this paper, we take advantage of the new JWST data for studying galaxies that are red in the near- to mid-IR color (some of them HST-dark), benefiting from two different aspects. First, we exploit the deeper flux limits provided by JWST and extend the selection of red mid-IR-detected galaxies carried out by previously published works, mainly using Spitzer/IRAC and HST/WFC3 data. And second, we also exploit the exquisite spatial resolution provided by JWST and HST in an unprecedentedly wide spectral range covering from the observed optical wavelengths (around $0.6~\mu$m) to the observed mid-IR ($\sim 4.4~\mu$m). We remark that this spectral range has been used extensively to study the integrated light of galaxies at cosmological distances up to $z \sim 10$, and we are now able to use the same type of spectral energy distributions (SEDs) but with much better spatial resolution in the entire (and very common in the last 20 yr since Spitzer was launched) optical-to-mid-IR spectral window thanks to JWST.

With this methodology, our goals are twofold. First, we aim to robustly characterize the redshift, stellar mass, and evolutionary stage of HST-dark galaxies, and mid-to-near-IR red galaxies in general. Within this goal, we are especially interested in identifying the first massive galaxies ever formed and, particularly, evolved galaxies that quenched $0.5-1.0$ Gyr after a starburst. Such systems are classically called quiescent galaxies. We note that rather than the term quiescent, it may be more adequate at early epochs to talk about dormant galaxies; since the universe itself is not much older than those ages, there is no time to be quiescent (in the same way used to talk about other galaxy populations at lower redshifts, such as red nuggets; Costantin et al. 2020; Tortora et al. 2020; ; Lisiecki et al. 2023) but these early galaxies may reignite and rejuvenate in their future evolution (Pérez-González et al. 2008a; Fumagalli et al. 2014; Tacchella et al. 2022b; Woodrum et al. 2022).

Our second goal is to benefit from the $8 \times$ jump in spatial resolution provided by JWST in the mid-IR (with respect to Spitzer), and $1 \times -2 \times$ in the near-infrared (with respect to HST). With those enhanced capabilities, we will study the two-dimensional distribution of the stellar populations in HST-dark galaxies, searching for the location of dusty starbursts, clues into their mass-growth and structural assembly channels, and, more generally, information about the first episodes in the star formation history of massive galaxies occurring at the earliest cosmic epochs.

This paper is organized as follows. Section 2 presents the observations used in our study. Section 3 discusses the selection of near-IR-faint, mid-IR-bright galaxies and the expectations for their nature. In Section 4, we describe our photometric measurements. Next, Section 5 presents the methodology to analyze our sample using spatially resolved spectral energy distributions to obtain photometric redshifts and stellar population properties. Section 6 presents our results about the integrated properties of our sample as well as the radial gradients of the stellar populations they harbor and the morphologies of our galaxies. Finally, Section 7 summarizes our results and conclusions.

Throughout the paper, we assume a flat cosmology with $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$, and a Hubble constant $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. We use AB magnitudes (Oke & Gunn 1983). All stellar-mass and SFR estimations assume a universal Chabrier (2003) IMF.

## 2. Data

In this paper, we use the JWST imaging data acquired within the Cosmic Evolution Early Released Science (CEERS) project (Finkelstein et al. 2017) in 2022 June. This data set comprises four NIRCam pointings (1, 2, 3, and 6, as named in the APT observing strategy file) in the Extended Groth Strip (EGS, Davis et al. 2007), observed in filters F115W, F150W, F200W, F277W, F356W, F410M, and F444W. The total surveyed area is $38'8^2$. The $5\sigma$ depths in each filter for point-like sources using $0.1''$15 circular apertures are 29.24, 29.08, 29.23, 30.08, 30.13, 29.74, and 29.50 mag; respectively, and 28.96, 28.84, 29.02, 29.74, 29.77, and 29.34, 29.09 mag after applying aperture corrections. The $5\sigma$ limits for individual pixels (important for our two-dimensional study of the stellar populations in high-redshift galaxies) are 31.61, 31.45, 31.61, 32.45, 32.50, 32.11, and 31.85 mag for the filters mentioned above, respectively.

The NIRCam data were calibrated with JWST pipeline version 1.7.2, reference files in pmap version 0989 (which includes detector-to-detector matched, improved absolute photometric calibration; see M. Bagley et al. 2023, in preparation). Apart from the standard pipeline stages, we performed careful alignment of the three images available for each pointing and filter before stacking them. We also applied a background homogenization algorithm prior to obtaining the final mosaics. The whole data set was registered to the same World Coordinate System reference frame (based on Gaia DR1.2; Gaia Collaboration et al. 2016a, 2016b) and put to the same pixel scale, namely, $0.03 \text{ pixel}^{-1}$.

Complementary to the JWST observations, we also used the Hubble Space Telescope HST data provided by the AEGIS (Davis et al. 2007) and CANDELS (Grogin et al. 2011; Koekemoer et al. 2011) collaborations in the EGS region. The data were reduced again and drizzled to the same pixel size of the JWST observations, $0.03 \text{ pixel}^{-1}$, by the CEERS team (version v1.9; see M. Bagley et al. 2023, in preparation). We added to our analysis those bands with better spatial resolution than the worst obtained by JWST n filter F444W (see Section 4), namely F606W and F814W from the Advanced Camera for Surveys, i.e., we exclude WFC3 observations from our analysis. Depths for these images are 29.29 and 29.05 mag, respectively, for point-like sources in a $0.1''$15 circular aperture, accounting for negligible aperture corrections. The $5\sigma$ limits for individual pixels are 31.66 and 31.42 mag, respectively.
Indeed, with the advent of JWST, new mid-IR-bright, near-IR-faint galaxies have been identified and started to be characterized (Barrufet et al. 2022; Carnall et al. 2023; Labbe et al. 2022; Nelson et al. 2022). These works have confirmed the high-redshift nature of HST-dark galaxies, their relevance for compiling a complete stellar-mass census in the first 1–2 Gyr of the lifetime of the universe, and their challenging numbers and properties for galaxy evolution models. Disagreement with theoretical predictions has also been (preliminary) found with the large numbers of very high-redshift ($z > 10$) galaxy candidates being found in relatively shallow JWST/NIRCam data (e.g., Adams et al. 2023; Castellano et al. 2022; Donnan et al. 2023; Finkelstein et al. 2022; Harikane et al. 2022; Naidu et al. 2022; Ono et al. 2022; Rodighiero et al. 2023), implying high levels of star formation activity in the very early universe (even though these candidates must still be confirmed and separated from contaminants such as dusty galaxies at $2 < z < 5$, see Zavala et al. 2023).

In this paper, we take advantage of the new JWST data for studying galaxies that are red in the near- to mid-IR color (some of them HST-dark), benefiting from two different aspects. First, we exploit the deeper flux limits provided by JWST and extend the selection of red mid-IR-detected galaxies carried out by previously published works, mainly using Spitzer/IRAC and HST/WFC3 data. And second, we also exploit the exquisite spatial resolution provided by JWST and HST in an unprecedentedly wide spectral range covering from the observed optical wavelengths (around $0.6~\mu$m) to the observed mid-IR ($\sim 4.4~\mu$m). We remark that this spectral range has been used extensively to study the integrated light of galaxies at cosmological distances up to $z \sim 10$, and we are now able to use the same type of spectral energy distributions (SEDs) but with much better spatial resolution in the entire (and very common in the last 20 yr since Spitzer was launched) optical-to-mid-IR spectral window thanks to JWST.

With this methodology, our goals are twofold. First, we aim to robustly characterize the redshift, stellar mass, and evolutionary stage of HST-dark galaxies, and mid-to-near-IR red galaxies in general. Within this goal, we are especially interested in identifying the first massive galaxies ever formed and, particularly, evolved galaxies that quenched $0.5–1.0$ Gyr after a starburst. Such systems are classically called quiescent galaxies. We note that rather than the term quiescent, it may be more adequate at early epochs to talk about dormant galaxies; since the universe itself is not much older than those ages, there is no time to be quiescent (in the same way used to talk about other galaxy populations at lower redshifts, such as red nuggets; Costantin et al. 2020; Tortora et al. 2020; ; Lisiecki et al. 2023) but these early galaxies may reignite and rejuvenate in their future evolution (Pérez-González et al. 2008a; Fumagalli et al. 2014; Tacchella et al. 2022b; Woodrum et al. 2022).

Our second goal is to benefit from the $8 \times$ jump in spatial resolution provided by JWST in the mid-IR (with respect to Spitzer), and $1 \times -2 \times$ in the near-infrared (with respect to HST). With those enhanced capabilities, we will study the two-dimensional distribution of the stellar populations in HST-dark galaxies, searching for the location of dusty starbursts, clues into their mass-growth and structural assembly channels, and, more generally, information about the first episodes in the star formation history of massive galaxies occurring at the earliest cosmic epochs.

This paper is organized as follows. Section 2 presents the observations used in our study. Section 3 discusses the selection of near-IR-faint, mid-IR-bright galaxies and the expectations for their nature. In Section 4, we describe our photometric measurements. Next, Section 5 presents the methodology to analyze our sample using spatially resolved spectral energy distributions to obtain photometric redshifts and stellar population properties. Section 6 presents our results about the integrated properties of our sample as well as the radial gradients of the stellar populations they harbor and the morphologies of our galaxies. Finally, Section 7 summarizes our results and conclusions.

Throughout the paper, we assume a flat cosmology with $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$, and a Hubble constant $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. We use AB magnitudes (Oke & Gunn 1983). All stellar-mass and SFR estimations assume a universal Chabrier (2003) IMF.

## 2. Data

In this paper, we use the JWST imaging data acquired within the Cosmic Evolution Early Released Science (CEERS) project (Finkelstein et al. 2017) in 2022 June. This data set comprises four NIRCam pointings (1, 2, 3, and 6, as named in the APT observing strategy file) in the Extended Groth Strip (EGS, Davis et al. 2007), observed in filters F115W, F150W, F200W, F277W, F356W, F410M, and F444W. The total surveyed area is $38'8^2$. The $5\sigma$ depths in each filter for point-like sources using $0.1''$15 circular apertures are 29.24, 29.08, 29.23, 30.08, 30.13, 29.74, and 29.50 mag; respectively, and 28.96, 28.84, 29.02, 29.74, 29.77, and 29.34, 29.09 mag after applying aperture corrections. The $5\sigma$ limits for individual pixels (important for our two-dimensional study of the stellar populations in high-redshift galaxies) are 31.61, 31.45, 31.61, 32.45, 32.50, 32.11, and 31.85 mag for the filters mentioned above, respectively.

The NIRCam data were calibrated with JWST pipeline version 1.7.2, reference files in pmap version 0989 (which includes detector-to-detector matched, improved absolute photometric calibration; see M. Bagley et al. 2023, in preparation). Apart from the standard pipeline stages, we performed careful alignment of the three images available for each pointing and filter before stacking them. We also applied a background homogenization algorithm prior to obtaining the final mosaics. The whole data set was registered to the same World Coordinate System reference frame (based on Gaia DR1.2; Gaia Collaboration et al. 2016a, 2016b) and put to the same pixel scale, namely, $0.03 \text{ pixel}^{-1}$.

Complementary to the JWST observations, we also used the Hubble Space Telescope HST data provided by the AEGIS (Davis et al. 2007) and CANDELS (Grogin et al. 2011; Koekemoer et al. 2011) collaborations in the EGS region. The data were reduced again and drizzled to the same pixel size of the JWST observations, $0.03 \text{ pixel}^{-1}$, by the CEERS team (version v1.9; see M. Bagley et al. 2023, in preparation). We added to our analysis those bands with better spatial resolution than the worst obtained by JWST n filter F444W (see Section 4), namely F606W and F814W from the Advanced Camera for Surveys, i.e., we exclude WFC3 observations from our analysis. Depths for these images are 29.29 and 29.05 mag, respectively, for point-like sources in a $0.1''$15 circular aperture, accounting for negligible aperture corrections. The $5\sigma$ limits for individual pixels are 31.66 and 31.42 mag, respectively.
We cross-correlated our sample with the IR catalogs for Spitzer MIPS and Herschel PACS and SPIRE presented in Barro et al. (2019), constructed with a deconvolution methodology described in Pérez-González et al. (2010, see also Rawle et al. 2016; Rodríguez-Muñoz et al. 2019). The $5\sigma$ depths for these data are 45 µJy and 3.5 mJy for MIPS 24 µm and 70 µm, 8.7 and 13.1 mJy for PACS 100 and 160 µm, and 14.7, 17.3, and 17.9 mJy for SPIRE 250, 350, and 500 µm.

At longer wavelengths, we also used the submillimeter catalog presented in Zavala et al. (2017, 2018), which reaches a $5\sigma$ central depth of 0.2 mJy at 850 µm with SCUBA-2 as well as our NOEMA 1.1 mm follow-up project on a subsample of these SCUBA-2 galaxies (proposal W20CK; see also Zavala et al. 2023, and L. Ciesla et al. 2023, in preparation). To search for radio galaxies, we used the VLA/EGS 20 cm (1.4 GHz) survey described by Ivison et al. (2007). In addition, we searched for counterparts for our galaxies in the X-ray catalog presented in Nandra et al. (2015).

Finally, we used the spectroscopic redshifts compiled by Stefanon et al. (2017), adding new sources published by the MOSFIRE Deep Evolution Field (MOSDEF) team (Reddy et al. 2015; Kriek et al. 2015) in their 2021 data release. The spectroscopic sample in the four NIRCam fields includes 466 galaxies at a median and quartiles redshift $z = 0.81_{0.55}^{1.41}$, with 17% at $z > 2$, all below $z = 3.7$.

### 3. Sample Selection

In this section, we describe how, by using an F150W–F356W versus F356W color–magnitude diagram, we select a sample of red galaxies in the near-to-mid-infrared color, HST-dark systems among them. To motivate the discussion in the following sections, we also describe in some detail the redshifts and types of galaxies expected for a selection based on different thresholds in the color–magnitude diagram.

#### 3.1. Pre-JWST Searches for Mid-IR-bright, Optical/Near-IR-faint Galaxies

Previous works have shown that searching for massive, dusty, and quiescent galaxies at $z \gtrsim 3$ using color–magnitude diagrams is a very effective technique (Mancini et al. 2009; Huang et al. 2011; Caputi et al. 2012; Labbé et al. 2013; Nayyeri et al. 2014; Stefanon et al. 2015; Wang et al. 2016; Alcalde Pampliega et al. 2019). Briefly, these methods use a long-baseline optical/NIR color to identify galaxies with very red SEDs, which would be caused either by strong Balmer breaks/D4000 spectral features or large, dust-driven UV attenuations, which are more common in evolved, massive galaxies. Mid-to-near-IR large colors might also be present in galaxies with strong (high equivalent width) emission lines (see, e.g., Smit et al. 2014, 2015; Roberts-Borsani et al. 2016; Alcalde Pampliega et al. 2019). Prior to JWST, these methods typically relied on the HST/F160W filter for the blue band and either Spitzer/IRAC [3.6] or [4.5] channels for the red (and selection) band. The resulting samples find galaxies with relatively bright IRAC magnitudes ([3.6] and or [4.5] $\le 24.5$ mag) that are often very faint, even undetected, in the H band, implying F160W $> 26–27$ mag (depending on the observed field). Thus, they are commonly referred to as IRAC-bright, HST-dark (or -drops) galaxies (Caputi et al. 2012; Stefanon et al. 2015; Alcalde Pampliega et al. 2019).

The new and extremely deep near- and mid-IR first data sets from JWST reach more than 2 mag deeper in F150W and F356W than the previous HST and Spitzer photometry and, consequently, it has become easier to identify and characterize both IR-bright and much fainter HST-dark galaxies (Barrufet et al. 2022; Endsley et al. 2022; Nelson et al. 2022; Labbe et al. 2022).

Here we aim to study the properties of optical/near-IR faint galaxies selected with a similar method to those in previous papers but taking advantage of the new JWST capabilities. Therefore, we use the following thresholds in color–magnitude: F150W – F356W $> 1.5$ mag, F356W $< 27.5$ mag. The red color is similar to the values used in previous works (e.g., Wang et al. 2016; Alcalde Pampliega et al. 2019), while the much deeper limiting magnitude in F356W, compared to IRAC depths, [3.6] $\sim 24.5$ mag, implies that we can potentially find massive galaxies ($\log M_*/M_{\odot} > 11$) up to $z=8$ without including too many lower-redshift galaxies with $\log M_*/M_{\odot} < 10$, therefore keeping the overall sample focused on massive galaxies.

In order to understand better its heterogeneous nature, we divide the color-selected sample into two subsets, with magnitudes brighter and fainter than F150W = 26 mag, labeled HST-faint and HST-dark, respectively. The latter are a proxy for HST/WFC3 dropouts that were previously only detected in IRAC and/or IR/radio surveys, or completely unknown. The HST-faint sample, which consists mainly of HST/F160W-detected galaxies identified before in CANDELS and 3D-HST catalogs (Skelton et al. 2014; Stefanon et al. 2017), will be used for completeness, reference, and comparison to illustrate how the formal definition of an HST/WFC3-dropout (which depends on the depth of the survey and changes from paper to paper) affects the median redshifts and stellar population properties of an HST-dark sample (see next section). We do not consider here galaxies brighter than F150W $< 23$ mag, which are typically massive galaxies at $z \lesssim 2$ and have been well-characterized by HST and ground-based telescopes, even with spectroscopy.

#### 3.2. Understanding the Selection of Different Galaxy Types in the Color–Magnitude Diagram

Figure 1 shows our selection method. The left panel depicts the F150W–F356W versus F356W color–magnitude plot for the galaxies in the CEERS sample using a gray-scale hexdiagram to indicate the regions with a higher number density. The selection of sources and integrated photometry used in this plot have been performed with a hot+cold Sextractor run (Bertin & Arnouts 1996), using Kron (1980) apertures. The details of the catalog will be presented in S. Finkelstein et al. (2023, in preparation).

The color threshold used to select the sample is indicated with a dashed gray line, and the resulting sample of HST-faint and -dark galaxies is shown with circles and hexagons (light and dark colors, respectively). For comparison purposes, we show the sample of HST-dark galaxies from Alcalde Pampliega et al. (2019), selected in the CANDELS/GOODS region, with green circles. The overall distribution of those galaxies agrees very well with the location of our HST-dark sample at relatively bright F356W $\lesssim 25$ mag, except for the three bright and very red galaxies in the Alcalde Pampliega et al. (2019) sample, which have no equivalent in our JWST-based selection.
We cross-correlated our sample with the IR catalogs for Spitzer MIPS and Herschel PACS and SPIRE presented in Barro et al. (2019), constructed with a deconvolution methodology described in Pérez-González et al. (2010, see also Rawle et al. 2016; Rodríguez-Muñoz et al. 2019). The $5\sigma$ depths for these data are 45 µJy and 3.5 mJy for MIPS 24 µm and 70 µm, 8.7 and 13.1 mJy for PACS 100 and 160 µm, and 14.7, 17.3, and 17.9 mJy for SPIRE 250, 350, and 500 µm.

At longer wavelengths, we also used the submillimeter catalog presented in Zavala et al. (2017, 2018), which reaches a $5\sigma$ central depth of 0.2 mJy at 850 µm with SCUBA-2 as well as our NOEMA 1.1 mm follow-up project on a subsample of these SCUBA-2 galaxies (proposal W20CK; see also Zavala et al. 2023, and L. Ciesla et al. 2023, in preparation). To search for radio galaxies, we used the VLA/EGS 20 cm (1.4 GHz) survey described by Ivison et al. (2007). In addition, we searched for counterparts for our galaxies in the X-ray catalog presented in Nandra et al. (2015).

Finally, we used the spectroscopic redshifts compiled by Stefanon et al. (2017), adding new sources published by the MOSFIRE Deep Evolution Field (MOSDEF) team (Reddy et al. 2015; Kriek et al. 2015) in their 2021 data release. The spectroscopic sample in the four NIRCam fields includes 466 galaxies at a median and quartiles redshift $z = 0.81_{0.55}^{1.41}$, with 17% at $z > 2$, all below $z = 3.7$.

### 3. Sample Selection

In this section, we describe how, by using an F150W–F356W versus F356W color–magnitude diagram, we select a sample of red galaxies in the near-to-mid-infrared color, HST-dark systems among them. To motivate the discussion in the following sections, we also describe in some detail the redshifts and types of galaxies expected for a selection based on different thresholds in the color–magnitude diagram.

#### 3.1. Pre-JWST Searches for Mid-IR-bright, Optical/Near-IR-faint Galaxies

Previous works have shown that searching for massive, dusty, and quiescent galaxies at $z \gtrsim 3$ using color–magnitude diagrams is a very effective technique (Mancini et al. 2009; Huang et al. 2011; Caputi et al. 2012; Labbé et al. 2013; Nayyeri et al. 2014; Stefanon et al. 2015; Wang et al. 2016; Alcalde Pampliega et al. 2019). Briefly, these methods use a long-baseline optical/NIR color to identify galaxies with very red SEDs, which would be caused either by strong Balmer breaks/D4000 spectral features or large, dust-driven UV attenuations, which are more common in evolved, massive galaxies. Mid-to-near-IR large colors might also be present in galaxies with strong (high equivalent width) emission lines (see, e.g., Smit et al. 2014, 2015; Roberts-Borsani et al. 2016; Alcalde Pampliega et al. 2019). Prior to JWST, these methods typically relied on the HST/F160W filter for the blue band and either Spitzer/IRAC [3.6] or [4.5] channels for the red (and selection) band. The resulting samples find galaxies with relatively bright IRAC magnitudes ([3.6] and or [4.5] $\le 24.5$ mag) that are often very faint, even undetected, in the H band, implying F160W $> 26–27$ mag (depending on the observed field). Thus, they are commonly referred to as IRAC-bright, HST-dark (or -drops) galaxies (Caputi et al. 2012; Stefanon et al. 2015; Alcalde Pampliega et al. 2019).

The new and extremely deep near- and mid-IR first data sets from JWST reach more than 2 mag deeper in F150W and F356W than the previous HST and Spitzer photometry and, consequently, it has become easier to identify and characterize both IR-bright and much fainter HST-dark galaxies (Barrufet et al. 2022; Endsley et al. 2022; Nelson et al. 2022; Labbe et al. 2022).

Here we aim to study the properties of optical/near-IR faint galaxies selected with a similar method to those in previous papers but taking advantage of the new JWST capabilities. Therefore, we use the following thresholds in color–magnitude: F150W–F356W $> 1.5$ mag, F356W $< 27.5$ mag. The red color is similar to the values used in previous works (e.g., Wang et al. 2016; Alcalde Pampliega et al. 2019), while the much deeper limiting magnitude in F356W, compared to IRAC depths, [3.6] $\sim 24.5$ mag, implies that we can potentially find massive galaxies ($\log M_*/M_{\odot} > 11$) up to $z=8$ without including too many lower-redshift galaxies with $\log M_*/M_{\odot} < 10$, therefore keeping the overall sample focused on massive galaxies.

In order to understand better its heterogeneous nature, we divide the color-selected sample into two subsets, with magnitudes brighter and fainter than F150W = 26 mag, labeled HST-faint and HST-dark, respectively. The latter are a proxy for HST/WFC3 dropouts that were previously only detected in IRAC and/or IR/radio surveys, or completely unknown. The HST-faint sample, which consists mainly of HST/F160W-detected galaxies identified before in CANDELS and 3D-HST catalogs (Skelton et al. 2014; Stefanon et al. 2017), will be used for completeness, reference, and comparison to illustrate how the formal definition of an HST/WFC3-dropout (which depends on the depth of the survey and changes from paper to paper) affects the median redshifts and stellar population properties of an HST-dark sample (see next section). We do not consider here galaxies brighter than F150W $< 23$ mag, which are typically massive galaxies at $z \lesssim 2$ and have been well-characterized by HST and ground-based telescopes, even with spectroscopy.

#### 3.2. Understanding the Selection of Different Galaxy Types in the Color–Magnitude Diagram

Figure 1 shows our selection method. The left panel depicts the F150W–F356W versus F356W color–magnitude plot for the galaxies in the CEERS sample using a gray-scale hexdiagram to indicate the regions with a higher number density. The selection of sources and integrated photometry used in this plot have been performed with a hot+cold Sextractor run (Bertin & Arnouts 1996), using Kron (1980) apertures. The details of the catalog will be presented in S. Finkelstein et al. (2023, in preparation).

The color threshold used to select the sample is indicated with a dashed gray line, and the resulting sample of HST-faint and -dark galaxies is shown with circles and hexagons (light and dark colors, respectively). For comparison purposes, we show the sample of HST-dark galaxies from Alcalde Pampliega et al. (2019), selected in the CANDELS/GOODS region, with green circles. The overall distribution of those galaxies agrees very well with the location of our HST-dark sample at relatively bright F356W $\lesssim 25$ mag, except for the three bright and very red galaxies in the Alcalde Pampliega et al. (2019) sample, which have no equivalent in our JWST-based selection.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The left panel shows an F150W–F356W vs. F356W color–magnitude diagram. The overall population of galaxies in the CEERS sample is shown in a gray-scale density map (cut to F356W < 28 mag). The markers show the color-selected, F150W–F356W > 1.5 mag (horizontal dashed gray line) sample divided into HST-faint (light circles) and HST-dark (dark hexagons) galaxies with magnitudes brighter and fainter than F150W = 26 mag, respectively. Diagonal dashed gray lines at 26, 28, 30, and 32 mag are shown. Nearly 70% of the HST-dark galaxies are undetected in previous HST/WFC3-based surveys (~90% of those with F150W > 27 mag). Cyan circles and blue stars indicate galaxies within our sample detected in the far-IR and radio/submillimeter wavelengths, and green circles show the sample of HST-dark galaxies from Alcalde Pampliega et al. (2019). The middle panel shows the same plot with colored lines presenting the redshift dependence (marked with vertical lines) in color and magnitude for six different galaxy types with stellar masses $log M_*/M_⊙ = 11$ and 10, thin and thick lines, respectively. The typical SEDs and rest-frame UVJ colors of each type are indicated in the bottom and top panels on the right. The templates correspond to star-forming galaxies with low (blue tones) and high (magenta/purple) levels of obscuration, a quiescent galaxy template (red), and an extreme starburst spectrum with high equivalent width lines (black).</figcaption></figure>

To illustrate in more detail the types of galaxies and redshift ranges selected by the color–magnitude diagram and also highlight the advantages and shortcomings of this method, the middle panel of Figure 1 shows the tracks in the color–magnitude diagram of six different galaxy templates as a function of redshift. These galaxies are chosen to be (broadly) representative of the typical types in different regions of the UVJ diagram (top-right panel; Williams et al. 2009; Whitaker et al. 2010; Patel et al. 2012; Leja et al. 2019b; Zuckerman et al. 2021), i.e., they range from blue, low-mass, and low-attenuation galaxies (bottom-left region of the UVJ diagram; blue colors) to redder, more massive, dusty star-forming galaxies (SFGs; top right; purple colors) and quiescent galaxies (top left region of the UVJ diagram; red/orange colors). The average SEDs for these galaxies are shown in the bottom-right panel of Figure 1. The SEDs are normalized at ~1.6 µm rest frame to highlight the reddening of the UV/optical emission with increasing age and/or dust attenuation (from blue to purple templates).

In Figure 1, the thin/thick line-width color-redshift tracks show the effect of scaling the templates to low and high stellar masses, $log M_*/M_⊙ = 11$ and 10, respectively, which are roughly traced by the F356W magnitude (with some degeneracies and scatter). The low-mass/extinction templates are shown only for $log M_*/M_⊙ = 10$, since high dust contents typically exist at the high-mass end.

Overall, the color-redshift tracks illustrate that red colors F150W–F356W ≳ 1.5 mag and bright F356W magnitude provide an excellent proxy to identify massive, dusty SFGs at redshifts above $z \gtrsim 2$ and massive, quiescent galaxies at slightly higher redshifts above $z \gtrsim 3$. We note also that some of the bluest quiescent galaxies ($U-V \sim 1.3$ mag, orange color track) are still slightly bluer than F150W–F356W = 1.5 mag at $z = 3$, and were missed by color thresholds of previous HST

+Spitzer-based papers. For this paper, considering composite stellar population models with some additional residual recent star formation and not adding much mass (and thus not affecting the F356W flux significantly), we use a color cut $F150W-F356W = 1.5$ mag.

The color-redshift tracks also show that using a limiting magnitude down to F356W ~ 27 mag includes many more galaxies at higher redshifts and/or lower masses than previous IRAC-bright searches, which are capped at [3.6] ~ 24.5 mag (and in relatively large sky areas), often picked the most massive galaxies ($log M_*/M_⊙ \sim 11$) at redshifts of $z \gtrsim 3-4$ (Wang et al. 2016; Alcalde Pampliega et al. 2019). Indeed, the $log M_*/M_⊙ = 10$ color tracks (thick lines) indicate that the selection would be roughly complete down to that mass for quiescent galaxies up to $z \lesssim 8$ (red/orange) and even for the dustiest galaxies up to $z \sim 6$ (purple/magenta). Similarly, at the high-mass end ($log M_*/M_⊙ > 11$; thin lines), the color-magnitude selection would identify all the quiescent or dusty galaxies up to $z \sim 6$. We note that the top-right corner of the figure, where we would expect the most massive dusty galaxies at $z > 5$, is nearly empty. This is likely due to the relatively low number density of such sources, and very massive galaxies in general (e.g., Muzzin et al. 2013; Stefanon et al. 2015) combined with the small area surveyed by the current CEERS pointings (~5 × smaller than for example the CANDELS/GOODS regions). We further discuss these implications in Section 6.2.

One potential concern for mass completeness would be the existence of massive galaxies at $z > 5-7$ with less canonical (composite) SEDs, such as the ones identified by Labbe et al. (2022). Those galaxies have mirror-flipped-L-shaped SEDs (V-shaped in $f_λ$) with red colors in the rest-frame near-IR ($V-J \sim 2$ mag) and blue colors in the rest-frame UV similar to those of our blue templates (black tracks and marker in
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The left panel shows an F150W–F356W vs. F356W color–magnitude diagram. The overall population of galaxies in the CEERS sample is shown in a gray-scale density map (cut to F356W < 28 mag). The markers show the color-selected, F150W–F356W > 1.5 mag (horizontal dashed gray line) sample divided into HST-faint (light circles) and HST-dark (dark hexagons) galaxies with magnitudes brighter and fainter than F150W = 26 mag, respectively. Diagonal dashed gray lines at 26, 28, 30, and 32 mag are shown. Nearly 70% of the HST-dark galaxies are undetected in previous HST/WFC3-based surveys (~90% of those with F150W > 27 mag). Cyan circles and blue stars indicate galaxies within our sample detected in the far-IR and radio/submillimeter wavelengths, and green circles show the sample of HST-dark galaxies from Alcalde Pampliega et al. (2019). The middle panel shows the same plot with colored lines presenting the redshift dependence (marked with vertical lines) in color and magnitude for six different galaxy types with stellar masses $log M_*/M_⊙ = 11$ and 10, thin and thick lines, respectively. The typical SEDs and rest-frame UVJ colors of each type are indicated in the bottom and top panels on the right. The templates correspond to star-forming galaxies with low (blue tones) and high (magenta/purple) levels of obscuration, a quiescent galaxy template (red), and an extreme starburst spectrum with high equivalent width lines (black).</figcaption></figure>

To illustrate in more detail the types of galaxies and redshift ranges selected by the color–magnitude diagram and also highlight the advantages and shortcomings of this method, the middle panel of Figure 1 shows the tracks in the color–magnitude diagram of six different galaxy templates as a function of redshift. These galaxies are chosen to be (broadly) representative of the typical types in different regions of the UVJ diagram (top-right panel; Williams et al. 2009; Whitaker et al. 2010; Patel et al. 2012; Leja et al. 2019b; Zuckerman et al. 2021), i.e., they range from blue, low-mass, and low-attenuation galaxies (bottom-left region of the UVJ diagram; blue colors) to redder, more massive, dusty star-forming galaxies (SFGs; top right; purple colors) and quiescent galaxies (top left region of the UVJ diagram; red/orange colors). The average SEDs for these galaxies are shown in the bottom-right panel of Figure 1. The SEDs are normalized at ~1.6 µm rest frame to highlight the reddening of the UV/optical emission with increasing age and/or dust attenuation (from blue to purple templates).

In Figure 1, the thin/thick line-width color-redshift tracks show the effect of scaling the templates to low and high stellar masses, $log M_*/M_⊙ = 11$ and 10, respectively, which are roughly traced by the F356W magnitude (with some degeneracies and scatter). The low-mass/extinction templates are shown only for $log M_*/M_⊙ = 10$, since high dust contents typically exist at the high-mass end.

Overall, the color-redshift tracks illustrate that red colors F150W–F356W ≳ 1.5 mag and bright F356W magnitude provide an excellent proxy to identify massive, dusty SFGs at redshifts above $z \gtrsim 2$ and massive, quiescent galaxies at slightly higher redshifts above $z \gtrsim 3$. We note also that some of the bluest quiescent galaxies ($U-V \sim 1.3$ mag, orange color track) are still slightly bluer than F150W–F356W = 1.5 mag at $z = 3$, and were missed by color thresholds of previous HST

+Spitzer-based papers. For this paper, considering composite stellar population models with some additional residual recent star formation and not adding much mass (and thus not affecting the F356W flux significantly), we use a color cut $F150W-F356W = 1.5$ mag.

The color-redshift tracks also show that using a limiting magnitude down to F356W ~27 mag includes many more galaxies at higher redshifts and/or lower masses than previous IRAC-bright searches, which are capped at [3.6] ~24.5 mag (and in relatively large sky areas), often picked the most massive galaxies ($log M_*/M_⊙ \sim 11$) at redshifts of $z \gtrsim 3-4$ (Wang et al. 2016; Alcalde Pampliega et al. 2019). Indeed, the $log M_*/M_⊙ = 10$ color tracks (thick lines) indicate that the selection would be roughly complete down to that mass for quiescent galaxies up to $z \lesssim 8$ (red/orange) and even for the dustiest galaxies up to $z \sim 6$ (purple/magenta). Similarly, at the high-mass end ($log M_*/M_⊙ > 11$; thin lines), the color-magnitude selection would identify all the quiescent or dusty galaxies up to $z \sim 6$. We note that the top-right corner of the figure, where we would expect the most massive dusty galaxies at $z > 5$, is nearly empty. This is likely due to the relatively low number density of such sources, and very massive galaxies in general (e.g., Muzzin et al. 2013; Stefanon et al. 2015) combined with the small area surveyed by the current CEERS pointings (~5 × smaller than for example the CANDELS/GOODS regions). We further discuss these implications in Section 6.2.

One potential concern for mass completeness would be the existence of massive galaxies at $z > 5-7$ with less canonical (composite) SEDs, such as the ones identified by Labbe et al. (2022). Those galaxies have mirror-flipped-L-shaped SEDs (V-shaped in $f_λ$) with red colors in the rest-frame near-IR ($V-J \sim 2$ mag) and blue colors in the rest-frame UV similar to those of our blue templates (black tracks and marker in
<table><caption>**Table 1**
Sample of Mid-IR Bright Near-IR Faint Galaxies Presented in this Paper</caption><thead><tr><th>ID</th><th>R.A. (J2000)<br>(degrees)</th><th>Decl. (J2000)<br>(degrees)</th><th>F356W<br>(mag)</th><th>F150W - F356W<br>(mag)</th><th>z</th><th>log M<sub>x</sub><br>(M<sub>☉</sub>)</th><th>logSFR(10 Myr)<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>Comment</th></tr></thead><tbody><tr><td>nircam1-349</td><td>214.95571202</td><td>+52.98342646</td><td>27.31 ± 0.06</td><td>>2.77</td><td>6.24<sup>+0.03</sup><sub>-0.02</sub></td><td>9.01 ± 0.09</td><td>-0.29 ± 0.10</td><td>XELG-z6</td></tr><tr><td>nircam1-1085</td><td>214.95787499</td><td>+52.98030148</td><td>23.31 ± 0.03</td><td>1.91 ± 0.04</td><td>3.52<sup>+0.04</sup><sub>-0.06</sub></td><td>10.46 ± 0.03</td><td>-0.76 ± 0.03</td><td>QG</td></tr><tr><td>nircam1-1623</td><td>214.92577403</td><td>+52.95443735</td><td>24.79 ± 0.03</td><td>2.11 ± 0.05</td><td>3.07<sup>+0.06</sup><sub>-0.13</sub></td><td>9.82 ± 0.03</td><td>1.47 ± 0.06</td><td>SFG</td></tr><tr><td>nircam1-1744</td><td>215.01127115</td><td>+53.01358970</td><td>24.79 ± 0.03</td><td>1.78 ± 0.04</td><td>5.17<sup>+0.00</sup><sub>-0.00</sub></td><td>9.96 ± 0.06</td><td>2.79 ± 0.05</td><td>SFG</td></tr><tr><td>nircam1-1752</td><td>215.01221476</td><td>+53.01469376</td><td>24.72 ± 0.03</td><td>1.56 ± 0.04</td><td>3.62<sup>+0.03</sup><sub>-0.13</sub></td><td>9.80 ± 0.03</td><td>1.53 ± 0.03</td><td>SFG</td></tr><tr><td>nircam1-1867</td><td>214.99840491</td><td>+53.00461910</td><td>25.95 ± 0.03</td><td>2.03 ± 0.06</td><td>6.52<sup>+0.03</sup><sub>-0.08</sub></td><td>9.63 ± 0.10</td><td>2.07 ± 0.04</td><td>XELG-z6</td></tr><tr><td>nircam1-2080</td><td>214.98181546</td><td>+52.99123394</td><td>22.32 ± 0.03</td><td>1.72 ± 0.04</td><td>3.17<sup>+0.18</sup><sub>-0.24</sub></td><td>10.69 ± 0.01</td><td>0.52 ± 0.06</td><td>QG</td></tr></tbody></table>

**Note.** Basic information about the sample of galaxies in this paper: ID, coordinates, magnitude and colors used in the selection, redshift, stellar mass, SFR (averaged in the last 10 Myr), and comments (including galaxy activity type based on SED-based sSFR estimations (see text for details; we also mark quiescence results based on the UVJ diagram), and detection by MIPS, submillimeter, and/or X-ray surveys).

(This table is available in its entirety in machine-readable form.)

Figure 1). The strong UV upturn and relatively red optical SED suggest that these are experiencing a recent burst of unobscured star formation but also have a preexisting more evolved or dust-enshrouded stellar population. However, the optical emission could be partially contaminated by the presence of high-equivalent-width emission lines, which makes the estimate of the stellar mass more uncertain (Endsley et al. 2022). The color tracks for these galaxies decline very quickly down to F150W - F356W ~ 1 mag at $z \gtrsim 5$ instead of plateauing at around F150W - F356W ~ 2.5 mag slightly under the quiescent and dusty tracks. As we will demonstrate in the following sections, we have several of these types of galaxies in our sample thanks to our color cut at F150W - F356W = 1.5 mag, but still, some more active systems could exist and be missed by our selection. If those galaxies are common at the highest redshifts, they will be systematically missed by the typical color-magnitude selections unless we decrease the color cut (at the expense of sample contamination by less extreme systems).

### 3.3. Overall Properties of the Color-selected Sample: HST-faint and -dark

Our color-magnitude-selected sample is composed of 138 galaxies, whose median/quartiles magnitude and colors are $\langle F356W \rangle = 24.2_{23.0}^{25.6}$ mag and $\langle F150W - F356W \rangle = 1.9_{1.6}^{2.3}$ mag. Half of the sample, 50% (36%), 69 (49) galaxies, are HST-dark, $F150W > 26(27)$ mag, the rest would qualify as HST-faint.

Among the HST-dark subsample, the fraction of sources included in HST catalogs of the EGS area (Skelton et al. 2014; Stefanon et al. 2017) is less than 28%, and all of these galaxies are brighter than $F160W = 27$ mag. This number goes up to 35% if we include an IRAC selection (Barro et al. 2011). Among the HST-faint subsample, 85% of galaxies were previously cataloged by HST-selected and/or IRAC-selected surveys.

### 3.4. Counterparts in Far-IR, Submillimeter, and X-Ray Catalogs

Out of the 138 galaxies in our sample, 32 of them (23%) can be identified with MIPS 24 µm sources in the catalog presented in Barro et al. (2019), searching in a $1''5$ radius region, 75% of them being in the HST-faint sample. One-sixth of them have MIPS 70 µm measurements, all with a signal-to-noise ratio (S/N) < 2. A visual inspection of the MIPS24 sources reveals that

11 of them (75% in the HST-faint sample) are very bright and isolated MIPS emitters, with an average flux of 67 µJy detected at the $9\sigma$ level and average redshift $\langle z \rangle = 2.8$. The rest of the counterparts, 21 galaxies (65% in the HST-faint sample), are quite faint and/or possibly blended with nearby objects, presenting an average flux of 50 µJy detected at the $7\sigma$ level and average redshift $\langle z \rangle = 3.2$.

Herschel detections (all above $4\sigma$) by PACS at 100 µm were found for seven galaxies (5% of the sample), average flux 3.1 mJy detected on average at the $7\sigma$ level, with two of them in the HST-dark sample. For PACS 160 µm, we found six galaxies with $S/N > 4$ measurements, average flux 7.5 mJy detected on average at the $7\sigma$ level, and one in the HST-dark sample. For SPIRE bands, we found three possible counterparts with $S/N > 4$, average flux 13 mJy detected at, on average, $4\sigma$ level, and one in the HST-dark sample.

A total of 10 galaxies in our sample (7% of the total) are likely associated with SCUBA-2 sources reported in Zavala et al. (2017, we note that this survey does not cover our whole area). Nine of them are unequivocally associated with the dusty galaxies via 1.1 mm high-angular-resolution observations obtained with NOEMA (Zavala et al. 2023; L. Ciesla et al. in preparation), while the remaining four sources lie within the $1''5$ of the 850 µm position. An additional three sources were found in the VLA 20 cm maps from Ivison et al. (2007), within a $1''5$ radius.

Finally, we also cross-correlated our sources with the X-ray catalog presented in Nandra et al. (2015) finding two (three) association(s) within $0''5$ ($1''6$).

The IDs, coordinates, and basic information about our sample are given in Table 1.

### 4. Photometry in Two Dimensions

Taking advantage of the unprecedented spatial resolution combined with a wide wavelength coverage and depth provided by JWST and HST, in this paper, we base our analysis on two-dimensional spectral energy distributions (2D-SEDs) rather than the more common study of the integrated emission of galaxies. We identify three benefits from this approach:

1. The spectral analysis is similar to what has been extensively done before in the literature based on HST and Spitzer/IRAC data for complete galaxies, i.e.,
<table><caption>**Table 1**
Sample of Mid-IR Bright Near-IR Faint Galaxies Presented in this Paper</caption><thead><tr><th>ID</th><th>R.A. (J2000)<br>(degrees)</th><th>Decl. (J2000)<br>(degrees)</th><th>F356W<br>(mag)</th><th>F150W - F356W<br>(mag)</th><th>z</th><th>log M<sub>x</sub><br>(M<sub>☉</sub>)</th><th>logSFR(10 Myr)<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>Comment</th></tr></thead><tbody><tr><td>nircam1-349</td><td>214.95571202</td><td>+52.98342646</td><td>27.31 ± 0.06</td><td>>2.77</td><td>6.24<sup>+0.03</sup><sub>-0.02</sub></td><td>9.01 ± 0.09</td><td>-0.29 ± 0.10</td><td>XELG-z6</td></tr><tr><td>nircam1-1085</td><td>214.95787499</td><td>+52.98030148</td><td>23.31 ± 0.03</td><td>1.91 ± 0.04</td><td>3.52<sup>+0.04</sup><sub>-0.06</sub></td><td>10.46 ± 0.03</td><td>-0.76 ± 0.03</td><td>QG</td></tr><tr><td>nircam1-1623</td><td>214.92577403</td><td>+52.95443735</td><td>24.79 ± 0.03</td><td>2.11 ± 0.05</td><td>3.07<sup>+0.06</sup><sub>-0.13</sub></td><td>9.82 ± 0.03</td><td>1.47 ± 0.06</td><td>SFG</td></tr><tr><td>nircam1-1744</td><td>215.01127115</td><td>+53.01358970</td><td>24.79 ± 0.03</td><td>1.78 ± 0.04</td><td>5.17<sup>+0.00</sup><sub>-0.00</sub></td><td>9.96 ± 0.06</td><td>2.79 ± 0.05</td><td>SFG</td></tr><tr><td>nircam1-1752</td><td>215.01221476</td><td>+53.01469376</td><td>24.72 ± 0.03</td><td>1.56 ± 0.04</td><td>3.62<sup>+0.03</sup><sub>-0.13</sub></td><td>9.80 ± 0.03</td><td>1.53 ± 0.03</td><td>SFG</td></tr><tr><td>nircam1-1867</td><td>214.99840491</td><td>+53.00461910</td><td>25.95 ± 0.03</td><td>2.03 ± 0.06</td><td>6.52<sup>+0.03</sup><sub>-0.08</sub></td><td>9.63 ± 0.10</td><td>2.07 ± 0.04</td><td>XELG-z6</td></tr><tr><td>nircam1-2080</td><td>214.98181546</td><td>+52.99123394</td><td>22.32 ± 0.03</td><td>1.72 ± 0.04</td><td>3.17<sup>+0.18</sup><sub>-0.24</sub></td><td>10.69 ± 0.01</td><td>0.52 ± 0.06</td><td>QG</td></tr></tbody></table>

**Note.** Basic information about the sample of galaxies in this paper: ID, coordinates, magnitude and colors used in the selection, redshift, stellar mass, SFR (averaged in the last 10 Myr), and comments (including galaxy activity type based on SED-based sSFR estimations (see text for details; we also mark quiescence results based on the UVJ diagram), and detection by MIPS, submillimeter, and/or X-ray surveys).

(This table is available in its entirety in machine-readable form.)

Figure 1). The strong UV upturn and relatively red optical SED suggest that these are experiencing a recent burst of unobscured star formation but also have a preexisting more evolved or dust-enshrouded stellar population. However, the optical emission could be partially contaminated by the presence of high-equivalent-width emission lines, which makes the estimate of the stellar mass more uncertain (Endsley et al. 2022). The color tracks for these galaxies decline very quickly down to F150W - F356W ~ 1 mag at $z \gtrsim 5$ instead of plateauing at around F150W - F356W ~ 2.5 mag slightly under the quiescent and dusty tracks. As we will demonstrate in the following sections, we have several of these types of galaxies in our sample thanks to our color cut at F150W - F356W = 1.5 mag, but still, some more active systems could exist and be missed by our selection. If those galaxies are common at the highest redshifts, they will be systematically missed by the typical color-magnitude selections unless we decrease the color cut (at the expense of sample contamination by less extreme systems).

### 3.3. Overall Properties of the Color-selected Sample: HST-faint and -dark

Our color-magnitude-selected sample is composed of 138 galaxies, whose median/quartiles magnitude and colors are $\langle F356W \rangle = 24.2_{23.0}^{25.6}$ mag and $\langle F150W - F356W \rangle = 1.9_{1.6}^{2.3}$ mag. Half of the sample, 50% (36%), 69 (49) galaxies, are HST-dark, $F150W > 26(27)$ mag, the rest would qualify as HST-faint.

Among the HST-dark subsample, the fraction of sources included in HST catalogs of the EGS area (Skelton et al. 2014; Stefanon et al. 2017) is less than 28%, and all of these galaxies are brighter than $F160W = 27$ mag. This number goes up to 35% if we include an IRAC selection (Barro et al. 2011). Among the HST-faint subsample, 85% of galaxies were previously cataloged by HST-selected and/or IRAC-selected surveys.

### 3.4. Counterparts in Far-IR, Submillimeter, and X-Ray Catalogs

Out of the 138 galaxies in our sample, 32 of them (23%) can be identified with MIPS 24 µm sources in the catalog presented in Barro et al. (2019), searching in a $1''5$ radius region, 75% of them being in the HST-faint sample. One-sixth of them have MIPS 70 µm measurements, all with a signal-to-noise ratio (S/N) < 2. A visual inspection of the MIPS24 sources reveals that

11 of them (75% in the HST-faint sample) are very bright and isolated MIPS emitters, with an average flux of 67 µJy detected at the $9\sigma$ level and average redshift $\langle z \rangle = 2.8$. The rest of the counterparts, 21 galaxies (65% in the HST-faint sample), are quite faint and/or possibly blended with nearby objects, presenting an average flux of 50 µJy detected at the $7\sigma$ level and average redshift $\langle z \rangle = 3.2$.

Herschel detections (all above $4\sigma$) by PACS at 100 µm were found for seven galaxies (5% of the sample), average flux 3.1 mJy detected on average at the $7\sigma$ level, with two of them in the HST-dark sample. For PACS 160 µm, we found six galaxies with $S/N > 4$ measurements, average flux 7.5 mJy detected on average at the $7\sigma$ level, and one in the HST-dark sample. For SPIRE bands, we found three possible counterparts with $S/N > 4$, average flux 13 mJy detected at, on average, $4\sigma$ level, and one in the HST-dark sample.

A total of 10 galaxies in our sample (7% of the total) are likely associated with SCUBA-2 sources reported in Zavala et al. (2017, we note that this survey does not cover our whole area). Nine of them are unequivocally associated with the dusty galaxies via 1.1 mm high-angular-resolution observations obtained with NOEMA (Zavala et al. 2023; L. Ciesla et al. in preparation), while the remaining four sources lie within the $1''5$ of the 850 µm position. An additional three sources were found in the VLA 20 cm maps from Ivison et al. (2007), within a $1''5$ radius.

Finally, we also cross-correlated our sources with the X-ray catalog presented in Nandra et al. (2015) finding two (three) association(s) within $0''5$ ($1''6$).

The IDs, coordinates, and basic information about our sample are given in Table 1.

### 4. Photometry in Two Dimensions

Taking advantage of the unprecedented spatial resolution combined with a wide wavelength coverage and depth provided by JWST and HST, in this paper, we base our analysis on two-dimensional spectral energy distributions (2D-SEDs) rather than the more common study of the integrated emission of galaxies. We identify three benefits from this approach:

1. The spectral analysis is similar to what has been extensively done before in the literature based on HST and Spitzer/IRAC data for complete galaxies, i.e.,
choosing a 2D approach does not imply losing photo-
metric bands now that we have JWST data.

2. Smaller parts of a galaxy (especially if it is significantly extended) should have less complex star formation histories (SFHs), so a 2D approach is expected to help in diminishing the effect of stellar population synthesis degeneracies, especially for high-redshift sources where the stellar population age is relatively more limited. More specifically, our 2D approach can produce relatively complex SFH for an entire galaxy (similar to “nonpara- metric” approaches; see, e.g., Pacifici et al. 2012; Leja et al. 2019a; Johnson et al. 2021; Ji & Giavalisco 2023) based on simple parameterizations (e.g., an exponentially decaying SFH) for smaller parts of it.

3. Independent estimations of the photometric redshift should allow us to obtain more accurate results, especially if the galaxies present stellar population gradients, as well as help with the segmentation of images.

On the con side of our 2D approach, we can mention:

1. Squeezing the 2D stellar population synthesis procedure to very small regions of a galaxy can be detrimental since the derived properties can be affected by stochasticity in the IMF and/or correlations of properties in physically connected regions (e.g., linked to the escape of ionizing photons from one region of the galaxy to another).

2. The low-surface-brightness component of galaxies cannot be studied with the same spatial resolution as the “core” of the galaxy.

3. The combination of results from individual pixels or resolution elements to recover the properties of the galaxy as a whole might not be straightforward (e.g., in the determination of the redshift of the galaxy based on regions with very different S/N photometry).

4. The faintest and smallest galaxies would now allow for a 2D approach.

5. The pixel-by-pixel approach needs well-determined PSFs and matching kernels as well as accurately aligned images, details that might not be as important when working with photometry using large apertures.

Overall, the advantages in the analysis of our galaxies are significant, so we follow this route and perform a two-dimensional analysis of the emission of galaxies, which we describe next. We invite the reader to see an evaluation of the method for *z* > 1 galaxies in García-Argumánez et al. (2022) and similar approaches based on HST data (Abraham et al. 1999; Wuyts et al. 2012). This paper is focused on presenting an interesting population of galaxies that JWST has revealed and allowed us to study in detail and will discuss the optimization of procedures for 2D SED-fitting in more detail in future work, when more comprehensive data sets and data calibrations are available.

The JWST and HST images described in Section 2 were
registered to the same World Coordinate System including a
local solution for each galaxy in the sample, in order to ensure
an accurate alignment on a pixel-by-pixel basis. To do that, we
considered a square region around each source and measured
centroids for all galaxies in the CANDELS catalog within that
region for all images in the different filters. Based on the cross-
correlation of those sources with respect to a reference band,
F444W in our case, we realigned and remapped all images to a

common grid. The accuracy of this alignment is better than
0.′′001, on average, with a scatter smaller than 0.′′006 (one-fifth
of the pixel). These statistics are independent of the band and
are not affected by more complex morphologies of galaxies
observed in the bluest bands, in part because the relative WCS
alignment of different JWST (and HST) bands is extremely
accurate (see M. Bagley et al. 2023, in preparation).

After aligning the images, we PSF-matched them using F444W as the reference (see example in Figure 2). This is the band with the most extended PSF in our data set, FWHM 0.′′16 (based on the actual PSF constructed with CEERS data). Using all known stars provided by the CANDELS catalog and fainter stars identified with JWST colors, we built PSFs for each band and kernels using photutils to match them to the F444W PSF.

We then measured photometry on a pixel-by-pixel basis using the locally registered, PSF-matched JWST and HST images. We built SEDs for each pixel and kept for analysis those which had more than three bands with fluxes measured at the 3σ level or better. For bands with low-significance fluxes, we assumed 5σ upper limits in all subsequent analyses, including photometric redshift estimation and stellar population synthesis modeling. Considering 3σ upper limits instead had a very small impact on photometric redshifts, with an average difference between estimations $\langle \Delta z/(1+z) \rangle = 0.02$ and 5 galaxies out of the full sample of 138 sources changing redshift by more than $\Delta z = 0.2$. Fixing the redshift, the effect of using 3σ upper limits in the stellar population synthesis is negligible since the fits are mostly dominated by actually measured fluxes.

We also measured integrated photometry based on a Kron (1980) aperture (marked in red in Figure 2). We note that the latter covers an area that is larger than the sum of the individual pixels for which we were able to build SEDs. The typical aperture correction from the sum of pixels to the Kron aperture is 7% for our sample.

The median number of pixels with usable SEDs per galaxy in our sample is 472 (approximately 20 resolution elements, defined as the area of the FWHM), with 85% (65%) of the sample having more than 100 (150) pixels (~4/7 resolution elements). We note that the pixels in our images correspond to physical sizes between 250 and 170 pc, and the F444W FWHM to 1.3 and 1 kpc, for the redshift range covered by our sample (2 ≲ z ≲ 7).

**5. Photometric Redshifts and Stellar Population Synthesis**
**Modeling in Two Dimensions**

The SEDs described in Section 4 were used to estimate photometric redshift and stellar populations properties on a pixel-by-pixel basis and also for the integrated aperture. We describe the procedures in the following subsections.

5.1. *Estimation of Photometric Redshifts*

Photometric redshifts for each pixel were estimated with a
modified version of the EAZY code (Brammer et al. 2008). The
modification consisted in allowing the template-fitting algo-
rithm to use (5σ) upper limits as an input, not allowing any fit
to provide brighter fluxes (achieved by penalizing the χ²
calculation) than those limits, and excluding the band in the χ²
calculation for templates providing lower fluxes (Mérida et al.
2023). The code was run with the template error feature
disabled (which penalized IRAC bands as a default; i.e.,
choosing a 2D approach does not imply losing photo-
metric bands now that we have JWST data.

2. Smaller parts of a galaxy (especially if it is significantly extended) should have less complex star formation histories (SFHs), so a 2D approach is expected to help in diminishing the effect of stellar population synthesis degeneracies, especially for high-redshift sources where the stellar population age is relatively more limited. More specifically, our 2D approach can produce relatively complex SFH for an entire galaxy (similar to “nonpara- metric” approaches; see, e.g., Pacifici et al. 2012; Leja et al. 2019a; Johnson et al. 2021; Ji & Giavalisco 2023) based on simple parameterizations (e.g., an exponentially decaying SFH) for smaller parts of it.

3. Independent estimations of the photometric redshift should allow us to obtain more accurate results, especially if the galaxies present stellar population gradients, as well as help with the segmentation of images.

On the con side of our 2D approach, we can mention:

1. Squeezing the 2D stellar population synthesis procedure to very small regions of a galaxy can be detrimental since the derived properties can be affected by stochasticity in the IMF and/or correlations of properties in physically connected regions (e.g., linked to the escape of ionizing photons from one region of the galaxy to another).

2. The low-surface-brightness component of galaxies cannot be studied with the same spatial resolution as the “core” of the galaxy.

3. The combination of results from individual pixels or resolution elements to recover the properties of the galaxy as a whole might not be straightforward (e.g., in the determination of the redshift of the galaxy based on regions with very different S/N photometry).

4. The faintest and smallest galaxies would now allow for a 2D approach.

5. The pixel-by-pixel approach needs well-determined PSFs and matching kernels as well as accurately aligned images, details that might not be as important when working with photometry using large apertures.

Overall, the advantages in the analysis of our galaxies are significant, so we follow this route and perform a two-dimensional analysis of the emission of galaxies, which we describe next. We invite the reader to see an evaluation of the method for *z* > 1 galaxies in García-Argumánez et al. (2022) and similar approaches based on HST data (Abraham et al. 1999; Wuyts et al. 2012). This paper is focused on presenting an interesting population of galaxies that JWST has revealed and allowed us to study in detail and will discuss the optimization of procedures for 2D SED-fitting in more detail in future work, when more comprehensive data sets and data calibrations are available.

The JWST and HST images described in Section 2 were
registered to the same World Coordinate System including a
local solution for each galaxy in the sample, in order to ensure
an accurate alignment on a pixel-by-pixel basis. To do that, we
considered a square region around each source and measured
centroids for all galaxies in the CANDELS catalog within that
region for all images in the different filters. Based on the cross-
correlation of those sources with respect to a reference band,
F444W in our case, we realigned and remapped all images to a

common grid. The accuracy of this alignment is better than
0.′′001, on average, with a scatter smaller than 0.′′006 (one-fifth
of the pixel). These statistics are independent of the band and
are not affected by more complex morphologies of galaxies
observed in the bluest bands, in part because the relative WCS
alignment of different JWST (and HST) bands is extremely
accurate (see M. Bagley et al. 2023, in preparation).

After aligning the images, we PSF-matched them using F444W as the reference (see example in Figure 2). This is the band with the most extended PSF in our data set, FWHM 0.′′16 (based on the actual PSF constructed with CEERS data). Using all known stars provided by the CANDELS catalog and fainter stars identified with JWST colors, we built PSFs for each band and kernels using photutils to match them to the F444W PSF.

We then measured photometry on a pixel-by-pixel basis using the locally registered, PSF-matched JWST and HST images. We built SEDs for each pixel and kept for analysis those which had more than three bands with fluxes measured at the 3σ level or better. For bands with low-significance fluxes, we assumed 5σ upper limits in all subsequent analyses, including photometric redshift estimation and stellar population synthesis modeling. Considering 3σ upper limits instead had a very small impact on photometric redshifts, with an average difference between estimations $\langle \Delta z/(1+z) \rangle = 0.02$ and 5 galaxies out of the full sample of 138 sources changing redshift by more than $\Delta z = 0.2$. Fixing the redshift, the effect of using 3σ upper limits in the stellar population synthesis is negligible since the fits are mostly dominated by actually measured fluxes.

We also measured integrated photometry based on a Kron (1980) aperture (marked in red in Figure 2). We note that the latter covers an area that is larger than the sum of the individual pixels for which we were able to build SEDs. The typical aperture correction from the sum of pixels to the Kron aperture is 7% for our sample.

The median number of pixels with usable SEDs per galaxy in our sample is 472 (approximately 20 resolution elements, defined as the area of the FWHM), with 85% (65%) of the sample having more than 100 (150) pixels (~4/7 resolution elements). We note that the pixels in our images correspond to physical sizes between 250 and 170 pc, and the F444W FWHM to 1.3 and 1 kpc, for the redshift range covered by our sample (2 ≲ z ≲ 7).

**5. Photometric Redshifts and Stellar Population Synthesis**
**Modeling in Two Dimensions**

The SEDs described in Section 4 were used to estimate photometric redshift and stellar populations properties on a pixel-by-pixel basis and also for the integrated aperture. We describe the procedures in the following subsections.

5.1. *Estimation of Photometric Redshifts*

Photometric redshifts for each pixel were estimated with a
modified version of the EAZY code (Brammer et al. 2008). The
modification consisted in allowing the template-fitting algo-
rithm to use (5σ) upper limits as an input, not allowing any fit
to provide brighter fluxes (achieved by penalizing the χ²
calculation) than those limits, and excluding the band in the χ²
calculation for templates providing lower fluxes (Mérida et al.
2023). The code was run with the template error feature
disabled (which penalized IRAC bands as a default; i.e.,
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Example of a two-dimensional photometric redshift determination for one of our extended sources (identified as two galaxies at different redshifts; see text for more details). The upper-left panel shows an RGB 8" × 8" postage stamp of nircam2-7122, using NIRCam filters F115W, F150W, and F200W. The top-right thumbnails show individual images of the source in all HST/ACS and JWST/NIRCam filters used in this work. The images have been PSF-matched to the F444W image. The red ellipse marks the aperture used for the integrated photometric measurements, and the blue ellipse depicts the aperture containing half of the total flux in F356W. On the bottom left, we show the histogram of the most probable photometric redshifts for all pixels with enough S/N (see text for details) within the integrated aperture. The combined redshift probability distribution function (zPDF) for all these pixels is shown in red, while the most probable redshift and its uncertainty for this source are marked and written in red. The inset shows the photometric redshift map. On the bottom right, we show the SEDs of two circular apertures positioned on the central region of the source and the northern part, respectively (see circular regions marked in the RGB image insets). In red, we show the sum of all photometric redshift templates fitting the corresponding pixels (offset to the average redshift of the region for clarity). The bottom-right insets in these SED plots show the individual probability distribution functions for all pixels within the aperture, plotted with a transparency parameter.</figcaption></figure>

NIRCam LW fluxes would also be affected) and not using any prior. We used v1.3 templates, which include a dusty galaxy and a high-EW emission-line galaxy spectrum.

After obtaining photometric redshifts for each pixel, we combined the probability distribution functions (zPDF) for all pixels belonging to a given galaxy using the procedure described in Dahlen et al. (2013). Briefly, different zPDFs are combined allowing for catastrophic fits as well as some degree of interdependency among the different estimations, parameterized with a quantity called $\alpha$, which can take values from 1 to the number of combined zPDFs. We note that, in our case, the estimation of photometric redshifts for nearby pixels should be correlated, since the FWHM of the PSF-matched data is 0.16 (5 pixels), but for most galaxies, we counted with completely independent estimations given their large extension. We made some tests on the parameter used to account for this interdependence of estimations and found no significant variations in the final most probable redshifts provided for each galaxy for $\alpha > 2$.

The procedure of photometric redshift estimation in two dimensions is exemplified in Figure 2 with a peculiar source, as explained next. We show one of the selected galaxies, nircam2-7122, which seems to have two stellar clumps in the JWST

data. The northern part, orange in the RGB image, was included in HST catalogs (e.g., in Stefanon et al. 2017). At wavelengths redder than 2 µm, another clump starts to be visible to the south, a region that gets brighter at longer wavelengths while the northern knot gets relatively dimmer (and almost disappears by 4 µm). Our 2D analysis indicates that the northern (already known) galaxy lies at a lower redshift, as revealed by the photo-z map, which counts with 4–5 resolution elements for the northern galaxy and 6–7 times more for the new JWST source. SEDs and photo-z fits are shown in Figure 2 for the centers of the two distinct galaxies. A similar analysis was performed in the rest of the sample, just keeping the regions that present consistent photometric redshifts in the study of the stellar populations presented in the following sections.

Our approach to analyzing the photo-z maps was that in order to consider possible contamination from an interloper at a different redshift, the region must have a size of at least two resolution elements (defined in terms of the FWHM). This allows us to get rid of noise, visible in some zones of the photo-z map presented in Figure 2. In that plot, we clearly identify a different galaxy or knot with this criterion to the north (but not in other zones). We must also consider that the map is
<figure><img src="image_5.png" /></figure>

**Figure 2.** Example of a two-dimensional photometric redshift determination for one of our extended sources (identified as two galaxies at different redshifts; see text for more details). The upper-left panel shows an RGB 8" × 8" postage stamp of nircam2-7122, using NIRCam filters F115W, F150W, and F200W. The top-right thumbnails show individual images of the source in all HST/ACS and JWST/NIRCam filters used in this work. The images have been PSF-matched to the F444W image. The red ellipse marks the aperture used for the integrated photometric measurements, and the blue ellipse depicts the aperture containing half of the total flux in F356W. On the bottom left, we show the histogram of the most probable photometric redshifts for all pixels with enough S/N (see text for details) within the integrated aperture. The combined redshift probability distribution function (zPDF) for all these pixels is shown in red, while the most probable redshift and its uncertainty for this source are marked and written in red. The inset shows the photometric redshift map. On the bottom right, we show the SEDs of two circular apertures positioned on the central region of the source and the northern part, respectively (see circular regions marked in the RGB image insets). In red, we show the sum of all photometric redshift templates fitting the corresponding pixels (offset to the average redshift of the region for clarity). The bottom-right insets in these SED plots show the individual probability distribution functions for all pixels within the aperture, plotted with a transparency parameter.

NIRCam LW fluxes would also be affected) and not using any prior. We used v1.3 templates, which include a dusty galaxy and a high-EW emission-line galaxy spectrum.

After obtaining photometric redshifts for each pixel, we combined the probability distribution functions (zPDF) for all pixels belonging to a given galaxy using the procedure described in Dahlen et al. (2013). Briefly, different zPDFs are combined allowing for catastrophic fits as well as some degree of interdependency among the different estimations, parameterized with a quantity called $\alpha$, which can take values from 1 to the number of combined zPDFs. We note that, in our case, the estimation of photometric redshifts for nearby pixels should be correlated, since the FWHM of the PSF-matched data is 0.16 (5 pixels), but for most galaxies, we counted with completely independent estimations given their large extension. We made some tests on the parameter used to account for this interdependence of estimations and found no significant variations in the final most probable redshifts provided for each galaxy for $\alpha > 2$.

The procedure of photometric redshift estimation in two dimensions is exemplified in Figure 2 with a peculiar source, as explained next. We show one of the selected galaxies, nircam2-7122, which seems to have two stellar clumps in the JWST

data. The northern part, orange in the RGB image, was included in HST catalogs (e.g., in Stefanon et al. 2017). At wavelengths redder than 2 µm, another clump starts to be visible to the south, a region that gets brighter at longer wavelengths while the northern knot gets relatively dimmer (and almost disappears by 4 µm). Our 2D analysis indicates that the northern (already known) galaxy lies at a lower redshift, as revealed by the photo-z map, which counts with 4–5 resolution elements for the northern galaxy and 6–7 times more for the new JWST source. SEDs and photo-z fits are shown in Figure 2 for the centers of the two distinct galaxies. A similar analysis was performed in the rest of the sample, just keeping the regions that present consistent photometric redshifts in the study of the stellar populations presented in the following sections.

Our approach to analyzing the photo-z maps was that in order to consider possible contamination from an interloper at a different redshift, the region must have a size of at least two resolution elements (defined in terms of the FWHM). This allows us to get rid of noise, visible in some zones of the photo-z map presented in Figure 2. In that plot, we clearly identify a different galaxy or knot with this criterion to the north (but not in other zones). We must also consider that the map is
presented in terms of the most probable photo-z, but each pixel
counts with a complete zPDF.

With this in mind, in order to assign redshifts to different
regions of a given galaxy, we considered a probability
threshold to decide when two zPDFs (constructed with several
pixels as explained above) provide incompatible redshifts. This
assumes that a zPDF would catch degeneracies in the photo-z
determination linked to stellar population gradients throughout
the galaxy (or belonging to different galaxies). We considered
as separate galaxies those regions where the combined zPDF
implied a probability smaller than 10% of the region being at
the redshift of the rest of the galaxy.

We remark that the identification of redshift interlopers in
this work implies that some regions of a given source were
removed from the subsequent analysis (as in Figure 11). This
means that the total masses could be underestimated if the
removed zone is indeed part of the galaxy, but the SFHs (or
other properties such as SFRs, sSFR, A(V), etc...) would be
robust for the retained galaxy area.

The evaluation and calibration of the 2D photometric
redshift determination, i.e., how to combine the information
from different pixels to obtain the best redshift for a given
galaxy as well as help with the segmentation of objects, would
need more extended data sets in terms of area and availability
of accurate measurements at higher redshifts than what current
spectroscopic samples provide (and would be extended in the
coming months with JWST observations).

As a first evaluation, we run the procedure through all the
 galaxies with spectroscopic redshifts in the four available
 CEERS pointings (see Section 2), just using the nine filters
 selected for our study of HST-dark/faint galaxies. These
 galaxies sum up 1.1 million pixels, for which we calculated
 photometric redshifts. We obtained σ<sub>NMAD</sub> values (as defined
 in Brammer et al. 2008) a factor of 2 better (0.035 versus
 0.018) than using the integrated SEDs.

5.2. *Stellar Population Synthesis Modeling*

Once we obtained a most probable redshift by combining the estimations for all pixels belonging to each galaxy, we fixed its value and analyzed the stellar populations. We fitted the pixel-by-pixel and integrated SEDs with the synthesizer code (Pérez-González et al. 2003; Pérez-González et al. 2008b) assuming a delayed exponential as the SFH, with timescale values τ between 100 Myr and 5 Gyr, ages between 1 Myr and the age of the universe at the redshift of the source, all discrete metallicities provided by the Bruzual & Charlot (2003) models, and a Calzetti et al. (2000) attenuation law with V-band extinction values, A(V) between 0 and 5 mag. We assume a universal initial mass function (IMF) that follows the Chabrier (2003) functional form. Nebular continuum and emission lines were added to the models as described in Pérez-González et al. (2008b). A Monte Carlo method was carried out in order to obtain uncertainties and account for degeneracies (see Domínguez Sánchez et al. 2016). We warn the reader that this kind of methodology would not account for systematic uncertainties such as those arising from the parameterization of the SFH, metallicity evolution, or the complexity of dust effects (in general, and as a function of position in the galaxy).

Figure 3 shows the results of the stellar population synthesis (SPS) modeling for the galaxy presented in Figure 2 after removing the northern component, which is revealed to lie at a different redshift by our 2D photo-z method. All relevant

quantities have been corrected by the aperture correction
mentioned in Section 4. The galaxy presents a knot with
significant amounts of dust, A(V) ~ 3 mag, separated from the
stellar-mass centroid, and (integrated) SFR around 10 M⊙ yr⁻¹.
The rest of the galaxy presents a less active and low-dust-
content disky structure (see also the UVJ diagram per pixel).
The SFH of the galaxy shows three main episodes of star
formation, one beyond z ~ 10, another (extended) at z ~ 5.5,
and the recent burst in the center. The galaxy is not detected by
MIPS, PACS, SPIRE, or (sub)millimeter wavelengths (nor in
X-ray) at the 3σ level or better. This galaxy is transiting
through a poststarburst phase with some residual dusty star
formation near the nuclear region, presenting mass-weighted
ages of up to a few hundred Myr.

The 2D UVJ diagram on the bottom right of Figure 3 shows that the full galaxy would be located in an intermediate region between blue, low-dust content, SFGs and dusty starbursts. The pixels near the mass-weighted center of the galaxy lie toward the dust-enshrouded zone, while outer pixels beyond the effective radius tend to lie in the blue cloud and poststarburst/quiescent wedge. In order to demonstrate the reliability of this spatially resolved UVJ diagram, we show color maps directly constructed with the JWST bands that roughly probe the UVJ bands, namely, F150W, F277W, and F444W. For the latter, given that we would need some extrapolation to probe the J band, we multiply it by a factor of 1.5, so we roughly match in the color map the V–J values in the color–color plot.

Our 2D-SPS methodology was compared in García-Argu-
mánez et al. (2022) with fits of the integrated emission of
simulated galaxies in Illustris using SFHs including a single
burst and two star formation events, each following a delayed
exponentially decaying SFH. We demonstrated that this more
classical approach could underestimate the mass-weighted ages
by factors of a few (see Figure 8 in García-Argumánez et al.
2022). The paper also showed that our 2D-SPS approach could
reproduce the Illustris SFHs with 20%–30% accuracies in
physical properties such as mass-weighted ages and formation
times. For this paper, we compared our 2D-SPS-based stellar
masses with those obtained from single-burst delayed expo-
nentially decaying SFHs fitting the integrated emission,
obtaining a very small systematic offset, Δlog M<sub>⋆</sub> = 0.02
dex, and a scatter that ranges from 0.05 dex at
log M<sub>⋆</sub>/M<sub>☉</sub> > 10 to 0.3 dex for smaller masses. A detailed
comparison of modeling procedures including 2D-SPS and
integrated photometry with parametric and nonparametric
SFHs will be presented in other papers (e.g., A. García-
Argumánez et al. 2023, in preparation).

6. Analysis of the Integrated and Bidimensional Stellar Population Properties

In this section, we analyze the photometric redshifts, stellar
masses, morphology, and both the integrated and spatially
resolved current versus past star formation activity of the 138
galaxies selected as “HST-dark” and “HST-faint” using the
mid-to-near-IR color–magnitude diagram.

6.1. Redshifts and Stellar Masses

The left panel of Figure 4 shows the photometric redshift
versus stellar-mass diagram for galaxies in the HST WFC3/
F160W-selected CANDELS catalog (gray-scale density map)
overlapping with the CEERS region, and the galaxy sample
presented in terms of the most probable photo-z, but each pixel counts with a complete zPDF.

With this in mind, in order to assign redshifts to different regions of a given galaxy, we considered a probability threshold to decide when two zPDFs (constructed with several pixels as explained above) provide incompatible redshifts. This assumes that a zPDF would catch degeneracies in the photo-z determination linked to stellar population gradients throughout the galaxy (or belonging to different galaxies). We considered as separate galaxies those regions where the combined zPDF implied a probability smaller than 10% of the region being at the redshift of the rest of the galaxy.

We remark that the identification of redshift interlopers in this work implies that some regions of a given source were removed from the subsequent analysis (as in Figure 11). This means that the total masses could be underestimated if the removed zone is indeed part of the galaxy, but the SFHs (or other properties such as SFRs, sSFR, $A(V)$, etc...) would be robust for the retained galaxy area.

The evaluation and calibration of the 2D photometric redshift determination, i.e., how to combine the information from different pixels to obtain the best redshift for a given galaxy as well as help with the segmentation of objects, would need more extended data sets in terms of area and availability of accurate measurements at higher redshifts than what current spectroscopic samples provide (and would be extended in the coming months with JWST observations).

As a first evaluation, we run the procedure through all the galaxies with spectroscopic redshifts in the four available CEERS pointings (see Section 2), just using the nine filters selected for our study of HST-dark/faint galaxies. These galaxies sum up 1.1 million pixels, for which we calculated photometric redshifts. We obtained $\sigma_{\text{NMAD}}$ values (as defined in Brammer et al. 2008) a factor of 2 better (0.035 versus 0.018) than using the integrated SEDs.

## 5.2. Stellar Population Synthesis Modeling

Once we obtained a most probable redshift by combining the estimations for all pixels belonging to each galaxy, we fixed its value and analyzed the stellar populations. We fitted the pixel-by-pixel and integrated SEDs with the synthesizer code (Pérez-González et al. 2003; Pérez-González et al. 2008b) assuming a delayed exponential as the SFH, with timescale values $\tau$ between 100 Myr and 5 Gyr, ages between 1 Myr and the age of the universe at the redshift of the source, all discrete metallicities provided by the Bruzual & Charlot (2003) models, and a Calzetti et al. (2000) attenuation law with V-band extinction values, $A(V)$ between 0 and 5 mag. We assume a universal initial mass function (IMF) that follows the Chabrier (2003) functional form. Nebular continuum and emission lines were added to the models as described in Pérez-González et al. (2008b). A Monte Carlo method was carried out in order to obtain uncertainties and account for degeneracies (see Domínguez Sánchez et al. 2016). We warn the reader that this kind of methodology would not account for systematic uncertainties such as those arising from the parameterization of the SFH, metallicity evolution, or the complexity of dust effects (in general, and as a function of position in the galaxy).

Figure 3 shows the results of the stellar population synthesis (SPS) modeling for the galaxy presented in Figure 2 after removing the northern component, which is revealed to lie at a different redshift by our 2D photo-z method. All relevant

quantities have been corrected by the aperture correction mentioned in Section 4. The galaxy presents a knot with significant amounts of dust, $A(V) \sim 3$ mag, separated from the stellar-mass centroid, and (integrated) SFR around $10 M_{\odot} \text{yr}^{-1}$. The rest of the galaxy presents a less active and low-dust-content disky structure (see also the UVJ diagram per pixel). The SFH of the galaxy shows three main episodes of star formation, one beyond $z \sim 10$, another (extended) at $z \sim 5.5$, and the recent burst in the center. The galaxy is not detected by MIPS, PACS, SPIRE, or (sub)millimeter wavelengths (nor in X-ray) at the $3\sigma$ level or better. This galaxy is transiting through a poststarburst phase with some residual dusty star formation near the nuclear region, presenting mass-weighted ages of up to a few hundred Myr.

The 2D UVJ diagram on the bottom right of Figure 3 shows that the full galaxy would be located in an intermediate region between blue, low-dust content, SFGs and dusty starbursts. The pixels near the mass-weighted center of the galaxy lie toward the dust-enshrouded zone, while outer pixels beyond the effective radius tend to lie in the blue cloud and poststarburst/quiescent wedge. In order to demonstrate the reliability of this spatially resolved UVJ diagram, we show color maps directly constructed with the JWST bands that roughly probe the UVJ bands, namely, F150W, F277W, and F444W. For the latter, given that we would need some extrapolation to probe the J band, we multiply it by a factor of 1.5, so we roughly match in the color map the V-J values in the color-color plot.

Our 2D-SPS methodology was compared in García-Argumánez et al. (2022) with fits of the integrated emission of simulated galaxies in Illustris using SFHs including a single burst and two star formation events, each following a delayed exponentially decaying SFH. We demonstrated that this more classical approach could underestimate the mass-weighted ages by factors of a few (see Figure 8 in García-Argumánez et al. 2022). The paper also showed that our 2D-SPS approach could reproduce the Illustris SFHs with 20%-30% accuracies in physical properties such as mass-weighted ages and formation times. For this paper, we compared our 2D-SPS-based stellar masses with those obtained from single-burst delayed exponentially decaying SFHs fitting the integrated emission, obtaining a very small systematic offset, $\Delta \log M_{\star} = 0.02$ dex, and a scatter that ranges from 0.05 dex at $\log M_{\star}/M_{\odot} > 10$ to 0.3 dex for smaller masses. A detailed comparison of modeling procedures including 2D-SPS and integrated photometry with parametric and nonparametric SFHs will be presented in other papers (e.g., A. García-Argumánez et al. 2023, in preparation).

# 6. Analysis of the Integrated and Bidimensional Stellar Population Properties

In this section, we analyze the photometric redshifts, stellar masses, morphology, and both the integrated and spatially resolved current versus past star formation activity of the 138 galaxies selected as “HST-dark” and “HST-faint” using the mid-to-near-IR color-magnitude diagram.

## 6.1. Redshifts and Stellar Masses

The left panel of Figure 4 shows the photometric redshift versus stellar-mass diagram for galaxies in the HST WFC3/F160W-selected CANDELS catalog (gray-scale density map) overlapping with the CEERS region, and the galaxy sample
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Results of our pixel-by-pixel SED-modeling procedure for the example source presented in Figure 2. In this figure, we just keep the southern component of the source, which we identify as a distinct galaxy that meets our selection criteria and lies at a different redshift than the northern region of the source. The top four panels show maps of the stellar-mass surface density, the mass-weighted age, the SFR averaged over the last 25 Myr of the SFH of the galaxy, and the dust attenuation in the V band. The centroid of the galaxy calculated with the stellar-mass map is marked with a gray star. The striped circle depicts the FWHM of the data. On the bottom left, we show the integrated SFH and its uncertainty obtained by adding the results for all pixels, after applying an aperture correction based on the average flux ratio between the sum of all pixels used in the 2D analysis and the integrated aperture. On the bottom-center panel, we show a UVJ diagram for all pixels used in our analysis, color-coded by distance to the galaxy centroid. Average uncertainties in both axis as a function of radius are plotted in the top-left (uncertainties in *V* − *J*) and bottom-right (in *U* − *V*) corners of the panel. The position of the integrated colors of the galaxy is marked with a black circle. On the bottom-right panel, we show color maps using the closest observed filters to the UVJ rest-frame colors considering the redshift of the galaxy (z = 3.7), namely F150W − F277W for *U* − *V* and F277W − 1.5 × F444W for *V* − *J* (the multiplicative factor being applied to account for the extrapolation to probe the *J* band).</figcaption></figure>

selected by our JWST color–magnitude criterion. We also distinguish between JWST galaxies that are brighter and fainter than F150W = 26 mag, namely, HST-faint and -dark galaxies (light and dark colors, hexagons and circles), respectively. At first glance, the main difference in the distribution of galaxies selected with HST and the new JWST-selected sample is that the latter reveals a population of red, massive (log M*/M⊙ = 9–11) galaxies at z ≳ 2, many of which were previously undetected (HST-dark) even in the deepest HST/WFC3 surveys such as CANDELS.

In more detail, the HST-faint sample selects predominantly massive galaxies (log M*/M☉ ≥ 10) over a redshift range around z ~ 3. This is consistent with the distribution of points for HST-faint galaxies in the color–magnitude diagram presented in Figure 1. As illustrated by the template color tracks, similarly massive galaxies at redshifts z ≲ 2 are not in the HST-faint sample because they have bluer colors

(F150W − F356W < 1.5 mag) and/or brighter F150W mag, while massive galaxies at higher redshifts tend to be redder and fainter, and thus fall into the selection region of the HST-dark sample.

Consequently, the HST-dark sample exhibits a higher median redshift *z* ~ 4.4 and a broader distribution that extends up to *z* ~ 7. The redshift overlap between HST-faint and -dark samples at *z* ~ 3.5 (central panel of Figure 4) is a direct consequence of the magnitude limit used to define the two sets (F150W = 26 mag). If we choose a fainter magnitude, many of the HST-dark galaxies at higher-z will move to the HST-faint sample, increasing its median redshift and lowering the median redshift of the HST-dark sample. We note that while F150W = 26 mag is an adequate limit to identify WFC3 dropouts in the shallower CANDELS mosaics in EGS, the deeper WFC3 imaging in GOODS/UDF fields (see, e.g., Guo et al. 2013; Barro et al. 2019) means that HST-dark galaxies in
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Results of our pixel-by-pixel SED-modeling procedure for the example source presented in Figure 2. In this figure, we just keep the southern component of the source, which we identify as a distinct galaxy that meets our selection criteria and lies at a different redshift than the northern region of the source. The top four panels show maps of the stellar-mass surface density, the mass-weighted age, the SFR averaged over the last 25 Myr of the SFH of the galaxy, and the dust attenuation in the V band. The centroid of the galaxy calculated with the stellar-mass map is marked with a gray star. The striped circle depicts the FWHM of the data. On the bottom left, we show the integrated SFH and its uncertainty obtained by adding the results for all pixels, after applying an aperture correction based on the average flux ratio between the sum of all pixels used in the 2D analysis and the integrated aperture. On the bottom-center panel, we show a UVJ diagram for all pixels used in our analysis, color-coded by distance to the galaxy centroid. Average uncertainties in both axis as a function of radius are plotted in the top-left (uncertainties in *V* − *J*) and bottom-right (in *U* − *V*) corners of the panel. The position of the integrated colors of the galaxy is marked with a black circle. On the bottom-right panel, we show color maps using the closest observed filters to the UVJ rest-frame colors considering the redshift of the galaxy (z = 3.7), namely F150W − F277W for *U* − *V* and F277W − 1.5 × F444W for *V* − *J* (the multiplicative factor being applied to account for the extrapolation to probe the *J* band).</figcaption></figure>

selected by our JWST color–magnitude criterion. We also distinguish between JWST galaxies that are brighter and fainter than F150W = 26 mag, namely, HST-faint and -dark galaxies (light and dark colors, hexagons and circles), respectively. At first glance, the main difference in the distribution of galaxies selected with HST and the new JWST-selected sample is that the latter reveals a population of red, massive (log M*/M⊙ = 9–11) galaxies at z ≳ 2, many of which were previously undetected (HST-dark) even in the deepest HST/WFC3 surveys such as CANDELS.

In more detail, the HST-faint sample selects predominantly massive galaxies (log M*/M☉ ≥ 10) over a redshift range around z ~ 3. This is consistent with the distribution of points for HST-faint galaxies in the color–magnitude diagram presented in Figure 1. As illustrated by the template color tracks, similarly massive galaxies at redshifts z < ∼ 2 are not in the HST-faint sample because they have bluer colors

(F150W − F356W < 1.5 mag) and/or brighter F150W mag, while massive galaxies at higher redshifts tend to be redder and fainter, and thus fall into the selection region of the HST-dark sample.

Consequently, the HST-dark sample exhibits a higher median redshift *z* ~ 4.4 and a broader distribution that extends up to *z* ~ 7. The redshift overlap between HST-faint and -dark samples at *z* ~ 3.5 (central panel of Figure 4) is a direct consequence of the magnitude limit used to define the two sets (F150W = 26 mag). If we choose a fainter magnitude, many of the HST-dark galaxies at higher-z will move to the HST-faint sample, increasing its median redshift and lowering the median redshift of the HST-dark sample. We note that while F150W = 26 mag is an adequate limit to identify WFC3 dropouts in the shallower CANDELS mosaics in EGS, the deeper WFC3 imaging in GOODS/UDF fields (see, e.g., Guo et al. 2013; Barro et al. 2019) means that HST-dark galaxies in
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The left panel shows the stellar-mass vs. photometric redshift diagram for all the galaxies in the CANDELS catalog (gray-scale density map) and the color-magnitude-selected sample divided into HST-dark (dark hexagons) and HST-faint (light circles). MIPS24 emitters are marked with a cyan circle, submillimeter and radio galaxies with a blue star. The central and right panels show the photometric redshift and stellar-mass distributions for the HST-dark and HST-faint galaxies (with the same color scheme) and the whole sample (gray histogram). Some statistical information (median and quartiles) for the HST-faint and HST-dark subsamples are provided.</figcaption></figure>

those fields will be fainter in F160W (by 0.5 to 1 mag) and will have, on average, higher redshifts (as presented exactly in those fields in, e.g., Alcalde Pampliega et al. 2019).

The redshift histogram for HST-dark galaxies reveals a peak around $z \sim 6.5$. The sources in this peak are peculiar, as we will also show in the coming sections. Here we emphasize that these sources typically present a blue F356W–F410M color, on average $-0.3 \pm 0.2$ mag, compared with an average of $+0.2 \pm 0.2$ mag for the rest of the sample. Some of them are also brighter in the F410M compared to both F356W and F444W.

In terms of the stellar-mass distribution (right panel of Figure 4), the HST-dark sample has a lower average stellar mass than the HST-faint sample, $\log M_*/M_{\odot} \sim 9.8$ versus $\log M_*/M_{\odot} \sim 10.4$, although both distributions have galaxies with masses up to $\log M_*/M_{\odot} \sim 11.0$. The HST-dark sample also extends to lower stellar masses than the HST-faint sample. This is primarily due to the faint limiting magnitude of the selection in F356W, which, to first order, correlates with the stellar mass of the galaxies. As illustrated by the template color tracks in Figure 1, some of the galaxies in the HST-dark sample with magnitudes around F356W $\sim 27$ mag will have $\log M_*/M_{\odot} \lesssim 10$ at redshifts ranging from $z=3$ to 7.

Another consequence of having a fainter limiting magnitude is that the average stellar mass of our HST-dark sample is also lower than those of pre-JWST HST-dark studies, which selected galaxies using Spitzer/IRAC and thus were restricted to much brighter limiting magnitudes of $[3.6] \sim 24.0–24.5$ mag. Therefore, these samples had virtually no galaxies with $\log M_*/M_{\odot} \lesssim 10$. For example, Wang et al. (2016) or Alcalde Pampliega et al. (2019) find average stellar masses of the order of $\log M_*/M_{\odot} = 10.5$ and $10.8$, respectively, for their HST-dark, IRAC-bright samples. A different effect that can contribute to the larger stellar masses of these IRAC-based studies is that they cover a larger area than the current CEERS survey by up to a factor of 4–5 and, therefore, they are more likely to find a higher fraction of rare galaxies with very large stellar masses. Furthermore, the number of those galaxies detected over the relatively small area surveyed so far by CEERS is more susceptible to cosmic variance.

## 6.2. Recent/Past Star Formation and Colors

In this section, we study the star formation properties and rest-frame colors of our sample in the context of the star formation main sequence and the UVJ diagram. Based on this analysis, we will divide the sample into star formation activity subsets that will be explored in the upcoming sections to analyze their stacked properties (e.g., star formation or mass-assembly histories). In particular, we will classify the sample into three classes: (1) quiescent sources; (2) SFGs, many of them detected at far-infrared and submillimeter wavelengths (note that the coverage of the CEERS footprint by submillimeter surveys is very inhomogeneous, so this is not a complete sample); and (3) peculiar galaxies at $z \sim 6$. We refer the reader to Section 6.2.3 for more details.

### 6.2.1. Recent Star Formation Properties and the Star-forming Main Sequence

Figure 5 shows the SFR versus stellar-mass plot for our sample, known as the star formation main-sequence (SFMS) diagram. In order to further understand the properties of our sample, which derive directly from the selection method, we divide the sample into several redshift bins in the following plots. The selected redshift bins are: one dominated by HST-faint galaxies at $2 < z < 3$, one dominated by HST-dark galaxies at $4 < z < 7$, and an intermediate range with a similar number of galaxies from the two subsamples, $3 < z < 4$. SFRs have been derived by calculating the average value from the SFHs of each galaxy in different time intervals, ranging from the last 5 to the last 100 Myr (namely, 5, 10, 25, 50, 75, and 100 Myr). In this figure, we use the SFR defined as that averaged over the last 10 Myr. Similar trends are observed in plots using other definitions of the SFR.

In Figure 5, we see our sample forming a band in the SFR-stellar-mass plane around the SFMS; most galaxies would qualify as SFGs and a small fraction presents very small specific SFRs (a $0.1 \text{ Gyr}^{-1}$ limit is marked in the plot), expediting the classification as quiescent/dormant galaxy. Although our SFR and stellar-mass estimations do follow a sequence, many of the galaxies lie below the “classical” SFMS (this is especially clear when using SFRs averaged on the last 75–100 Myr), typically obtained by estimating stellar masses
<figure><img src="image_7.png" /></figure>

**Figure 4.** The left panel shows the stellar-mass vs. photometric redshift diagram for all the galaxies in the CANDELS catalog (gray-scale density map) and the color-magnitude-selected sample divided into HST-dark (dark hexagons) and HST-faint (light circles). MIPS24 emitters are marked with a cyan circle, submillimeter and radio galaxies with a blue star. The central and right panels show the photometric redshift and stellar-mass distributions for the HST-dark and HST-faint galaxies (with the same color scheme) and the whole sample (gray histogram). Some statistical information (median and quartiles) for the HST-faint and HST-dark subsamples are provided.

those fields will be fainter in F160W (by 0.5 to 1 mag) and will have, on average, higher redshifts (as presented exactly in those fields in, e.g., Alcalde Pampliega et al. 2019).

The redshift histogram for HST-dark galaxies reveals a peak around $z \sim 6.5$. The sources in this peak are peculiar, as we will also show in the coming sections. Here we emphasize that these sources typically present a blue F356W–F410M color, on average $-0.3 \pm 0.2$ mag, compared with an average of $+0.2 \pm 0.2$ mag for the rest of the sample. Some of them are also brighter in the F410M compared to both F356W and F444W.

In terms of the stellar-mass distribution (right panel of Figure 4), the HST-dark sample has a lower average stellar mass than the HST-faint sample, $\log M_*/M_{\odot} \sim 9.8$ versus $\log M_*/M_{\odot} \sim 10.4$, although both distributions have galaxies with masses up to $\log M_*/M_{\odot} \sim 11.0$. The HST-dark sample also extends to lower stellar masses than the HST-faint sample. This is primarily due to the faint limiting magnitude of the selection in F356W, which, to first order, correlates with the stellar mass of the galaxies. As illustrated by the template color tracks in Figure 1, some of the galaxies in the HST-dark sample with magnitudes around F356W $\sim 27$ mag will have $\log M_*/M_{\odot} \lesssim 10$ at redshifts ranging from $z=3$ to 7.

Another consequence of having a fainter limiting magnitude is that the average stellar mass of our HST-dark sample is also lower than those of pre-JWST HST-dark studies, which selected galaxies using Spitzer/IRAC and thus were restricted to much brighter limiting magnitudes of $[3.6] \sim 24.0–24.5$ mag. Therefore, these samples had virtually no galaxies with $\log M_*/M_{\odot} \lesssim 10$. For example, Wang et al. (2016) or Alcalde Pampliega et al. (2019) find average stellar masses of the order of $\log M_*/M_{\odot} = 10.5$ and $10.8$, respectively, for their HST-dark, IRAC-bright samples. A different effect that can contribute to the larger stellar masses of these IRAC-based studies is that they cover a larger area than the current CEERS survey by up to a factor of 4–5 and, therefore, they are more likely to find a higher fraction of rare galaxies with very large stellar masses. Furthermore, the number of those galaxies detected over the relatively small area surveyed so far by CEERS is more susceptible to cosmic variance.

## 6.2. Recent/Past Star Formation and Colors

In this section, we study the star formation properties and rest-frame colors of our sample in the context of the star formation main sequence and the UVJ diagram. Based on this analysis, we will divide the sample into star formation activity subsets that will be explored in the upcoming sections to analyze their stacked properties (e.g., star formation or mass-assembly histories). In particular, we will classify the sample into three classes: (1) quiescent sources; (2) SFGs, many of them detected at far-infrared and submillimeter wavelengths (note that the coverage of the CEERS footprint by submillimeter surveys is very inhomogeneous, so this is not a complete sample); and (3) peculiar galaxies at $z \sim 6$. We refer the reader to Section 6.2.3 for more details.

### 6.2.1. Recent Star Formation Properties and the Star-forming Main Sequence

Figure 5 shows the SFR versus stellar-mass plot for our sample, known as the star formation main-sequence (SFMS) diagram. In order to further understand the properties of our sample, which derive directly from the selection method, we divide the sample into several redshift bins in the following plots. The selected redshift bins are: one dominated by HST-faint galaxies at $2 < z < 3$, one dominated by HST-dark galaxies at $4 < z < 7$, and an intermediate range with a similar number of galaxies from the two subsamples, $3 < z < 4$. SFRs have been derived by calculating the average value from the SFHs of each galaxy in different time intervals, ranging from the last 5 to the last 100 Myr (namely, 5, 10, 25, 50, 75, and 100 Myr). In this figure, we use the SFR defined as that averaged over the last 10 Myr. Similar trends are observed in plots using other definitions of the SFR.

In Figure 5, we see our sample forming a band in the SFR-stellar-mass plane around the SFMS; most galaxies would qualify as SFGs and a small fraction presents very small specific SFRs (a $0.1 \text{ Gyr}^{-1}$ limit is marked in the plot), expediting the classification as quiescent/dormant galaxy. Although our SFR and stellar-mass estimations do follow a sequence, many of the galaxies lie below the “classical” SFMS (this is especially clear when using SFRs averaged on the last 75–100 Myr), typically obtained by estimating stellar masses
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Main-sequence (SFR vs. stellar-mass, SFMS) plots for our sample of galaxies with red near-to-mid-IR colors, compared to the CANDELS sample. From left to right, each panel shows the SFMS in different redshift bins, selected to exemplify the range dominated by HST-faint galaxies (F150W < 26 mag), 2 < z < 3; HST-dark sources (F150W > 26 mag), 4 < z < 7; and a transition range. In each panel, our sample is shown with hexagons (HST-dark) and circles (HST-faint). MIPS24 emitters are marked with a cyan circle, submillimeter and radio galaxies with a blue star, quiescent galaxies with a red square, and XELG-z6 with green triangles. The hexmap shows the distribution of galaxies in the CANDELS catalog (Stefanon et al. 2017), with SFRs estimated with the method explained in Barro et al. (2019). The colored lines depict the SFMS at z = 2–7 from Schreiber et al. (2015), Santini et al. (2017), and Barro et al. (2019). The cutouts on the bottom show the 5" × 5" images of the quiescent galaxies marked with letters in the SFMS, also providing IDs and redshifts. From bottom to top, the images are color composites of three HST bands and three JWST bands, short and long wavelengths, respectively.</figcaption></figure>

from integrated-aperture SED fitting and SFRs from observed rest-frame UV luminosities corrected for attenuation with recipes based on UV slopes, or with SFR tracers linked to dust emission. At first sight, this kind of behavior (which has also been observed when comparing the observed SFMS with galaxy formation simulations such as Illustris; see Sparre et al. 2015; Donnari et al. 2019) might be interpreted as the SED-fitting analysis underestimating the SFRs due to high dust contents.

To probe the possible existence of highly embedded, optically thick star-forming knots that cannot be characterized with rest-frame UV-to-near-IR SEDs and to test how the 2D approach handles this classical problem, we discuss the effects of dust on the SFMS plot in Figure 6. For this purpose, we calculate the total energy absorbed by dust according to our pixel-by-pixel SPS analysis, adding up the contributions of all pixels for each galaxy to calculate the stellar dust-absorbed luminosity. Assuming an energy balance, this luminosity should be converted into dust emission in the mid- to far-IR. We then translate the dust-absorbed stellar luminosity to an IR-based SFR, applying the calibration in Kennicutt (1998). We remark that this exercise should provide results that are completely comparable to the regular technique used to construct the SFMS plot at the high-mass end. For low-mass galaxies, with lower dust contents, SFRs are typically based on rest-frame UV luminosities, maybe added to absorbed emission with hybrid SFR tracers (see Buat et al. 2005; Calzetti et al. 2007) or dust-corrected using UV slopes (with more or less complicated recipes, see Barro et al. 2019). Therefore, for our galaxies, which span a relatively wide range of masses, the

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Star formation main-sequence plot for our sample of mid-IR-bright, near-IR-faint galaxies constructed with SFRs derived from our 2D stellar population synthesis analysis through the calculation of the dust-absorbed stellar emission. The energy extincted by dust in the UV/optical is assumed to be reemitted in the mid- to far-IR, from 8 to 1000 µm. The luminosity in this wavelength range is then converted to SFR using the Kennicutt (1998) factor. Symbols in gray show the SFRs obtained from the average of the integrated SFHs in the last 20 Myr. Blue curves depict literature main-sequence relationships (same as in Figure 5).</figcaption></figure>

SFRs calculated with an energy-budget technique might be underestimated at the low stellar-mass end if significant amounts of nonobscured star formation exist.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Main-sequence (SFR vs. stellar-mass, SFMS) plots for our sample of galaxies with red near-to-mid-IR colors, compared to the CANDELS sample. From left to right, each panel shows the SFMS in different redshift bins, selected to exemplify the range dominated by HST-faint galaxies (F150W < 26 mag), 2 < z < 3; HST-dark sources (F150W > 26 mag), 4 < z < 7; and a transition range. In each panel, our sample is shown with hexagons (HST-dark) and circles (HST-faint). MIPS24 emitters are marked with a cyan circle, submillimeter and radio galaxies with a blue star, quiescent galaxies with a red square, and XELG-z6 with green triangles. The hexmap shows the distribution of galaxies in the CANDELS catalog (Stefanon et al. 2017), with SFRs estimated with the method explained in Barro et al. (2019). The colored lines depict the SFMS at z = 2–7 from Schreiber et al. (2015), Santini et al. (2017), and Barro et al. (2019). The cutouts on the bottom show the 5" × 5" images of the quiescent galaxies marked with letters in the SFMS, also providing IDs and redshifts. From bottom to top, the images are color composites of three HST bands and three JWST bands, short and long wavelengths, respectively.</figcaption></figure>

from integrated-aperture SED fitting and SFRs from observed rest-frame UV luminosities corrected for attenuation with recipes based on UV slopes, or with SFR tracers linked to dust emission. At first sight, this kind of behavior (which has also been observed when comparing the observed SFMS with galaxy formation simulations such as Illustris; see Sparre et al. 2015; Donnari et al. 2019) might be interpreted as the SED-fitting analysis underestimating the SFRs due to high dust contents.

To probe the possible existence of highly embedded, optically thick star-forming knots that cannot be characterized with rest-frame UV-to-near-IR SEDs and to test how the 2D approach handles this classical problem, we discuss the effects of dust on the SFMS plot in Figure 6. For this purpose, we calculate the total energy absorbed by dust according to our pixel-by-pixel SPS analysis, adding up the contributions of all pixels for each galaxy to calculate the stellar dust-absorbed luminosity. Assuming an energy balance, this luminosity should be converted into dust emission in the mid- to far-IR. We then translate the dust-absorbed stellar luminosity to an IR-based SFR, applying the calibration in Kennicutt (1998). We remark that this exercise should provide results that are completely comparable to the regular technique used to construct the SFMS plot at the high-mass end. For low-mass galaxies, with lower dust contents, SFRs are typically based on rest-frame UV luminosities, maybe added to absorbed emission with hybrid SFR tracers (see Buat et al. 2005; Calzetti et al. 2007) or dust-corrected using UV slopes (with more or less complicated recipes, see Barro et al. 2019). Therefore, for our galaxies, which span a relatively wide range of masses, the

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Star formation main-sequence plot for our sample of mid-IR-bright, near-IR-faint galaxies constructed with SFRs derived from our 2D stellar population synthesis analysis through the calculation of the dust-absorbed stellar emission. The energy extincted by dust in the UV/optical is assumed to be reemitted in the mid- to far-IR, from 8 to 1000 µm. The luminosity in this wavelength range is then converted to SFR using the Kennicutt (1998) factor. Symbols in gray show the SFRs obtained from the average of the integrated SFHs in the last 20 Myr. Blue curves depict literature main-sequence relationships (same as in Figure 5).</figcaption></figure>

SFRs calculated with an energy-budget technique might be underestimated at the low stellar-mass end if significant amounts of nonobscured star formation exist.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Rest-frame *U* − *V* vs. *V* − *J* diagram at different redshifts for the galaxies in the CEERS sample (color gray-scale) and the color-selected sample divided into HST-faint and -dark (light and dark colors). MIPS24 emitters are marked with a cyan circle, submillimeter and radio galaxies with a blue star, quiescent galaxies by sSFR with a red square, and XELG-z6 with green triangles. The green circles show quiescent galaxies at *z* > 3 from Carnall et al. (2023). HST-dark galaxies make up the bulk of the sample at *z* > 4 while HST-faint ones dominate in number at ⟨*z*⟩ = 2.5, and there is an even split between them at ⟨*z*⟩ = 3.5. The average color, extinction, and stellar mass are indicated in the top-left corner of each plot. Overall color-selected galaxies are predominantly massive, dusty, or quiescent galaxies in agreement with the goal of the selection criteria in Figure 1. Furthermore, there is a clear downward trend in extinction and stellar mass ranging from ⟨*A*<sub>V</sub>⟩ ∼ 2.2 mag and ⟨log *M*<sub>*</sub>/M<sub>☉</sub>⟩ ∼ 10.5 at *z* < 3 to ⟨*A*<sub>V</sub>⟩ ∼ 1.2 mag and ⟨log *M*<sub>*</sub>/M<sub>☉</sub>⟩ = 10 at *z* > 4. The cutouts on the bottom show the 5″ × 5″ images of the MIPS and submillimeter/radio galaxies marked with letters in the UVJ diagram, also providing IDs and redshifts. From bottom to top, the images are color composites of three HST bands and three JWST bands, short and long wavelengths, respectively.</figcaption></figure>

Figure 6 demonstrates that the 2D-based, global dust attenuations for the galaxies in our sample do reproduce the dust content and emission observed in galaxies and are used for the construction of the SFMS. Even though the 5–100 Myr time-averaged SFRs fall below the MS, the fair comparison with observations behaves much better. On average, and down to log M*/M☉ ~ 9.7, our energy-budget-based SFRs reproduce the SFMS, with the median sSFR lying within 0.1 dex of the literature curves. However, the scatter of our SFMS is 2× larger than the typical 0.2–0.3 dex reported in the literature. Whether this behavior is intrinsic or related to observational uncertainties is beyond the scope of this paper.

In Figure 6 (and the right panel of Figure 5), we mark the galaxies in the redshift peak at z ~ 6 reported in the previous section (green triangles). All of these sources present very large SFRs for their stellar mass, well above the SFMS, consequently qualifying as starburst galaxies. The SFRs derived with the energy-budget calculation are also large, indicating the presence of significant amounts of dust. We remark here that the starburst galaxies seem to follow a separate sequence, very similar to the bimodality reported in Rinaldi et al. (2022). We cannot claim any robust interpretation for this effect due to the small number of galaxies we have in our sample, but our modeling points to systematic effects linked to age (rather, for example, than for extinction) since most of the galaxies located far away from the SFMS present very young bursts, many of them belonging to the XELG-z6 type. This would be consistent with the presence of strong emission lines, most probably [O III]+Hβ, in these starburst systems, similar to the sources reported in Endsley et al. (2022), Matthee et al. (2022), and

Laporte et al. (2022), but they also present relatively high dust contents, with average attenuations around A(V) = 2 mag. Indeed, Matthee et al. (2022) report not negligible amounts of dust in some of the (NIRCam-grism-selected) spectroscopically confirmed [O III] emitters.

We conclude that our 2D-SPS method based on JWST+HST data does provide robust dust estimations at spatial resolutions around 200 pc and that SFRs from SED fitting should be compared with SFRs from classical tracers with caution.

### 6.2.2. Rest-frame UVJ Colors

Figure 7 shows the overall distribution of all the CEERS galaxies (gray-scale density map) and the color-selected sample, divided into HST-faint and -dark in the rest-frame UVJ diagrams at different redshifts. The distribution of rest-frame colors is roughly consistent with what we see in the SFMS plot: The bulk of the sample is massive, dusty SFGs with a small fraction of quiescent galaxies starting at z ≳ 3. The majority of the galaxies identified as quiescent by sSFR (red squares) are located in the quiescent region of the UVJ diagram, showing that the two methods are quite consistent.

At 2 < *z* < 3, nearly all the color-selected galaxies are HST-faint, dusty and star-forming, with high average extinctions ⟨*A*<sub>V</sub>⟩ = 2.2 mag and relatively large stellar masses of ⟨log *M*<sub>*</sub>/M<sub>☉</sub>⟩ = 10.4. This is consistent with the color tracks in Figure 1 that indicate that only dusty galaxies are redder than the selection threshold at *z* < 3. That is, quiescent galaxies at these redshifts do not make the color cut. Similarly, since these galaxies are also relatively massive, they tend to be on the
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Rest-frame *U* − *V* vs. *V* − *J* diagram at different redshifts for the galaxies in the CEERS sample (color gray-scale) and the color-selected sample divided into HST-faint and -dark (light and dark colors). MIPS24 emitters are marked with a cyan circle, submillimeter and radio galaxies with a blue star, quiescent galaxies by sSFR with a red square, and XELG-z6 with green triangles. The green circles show quiescent galaxies at *z* > 3 from Carnall et al. (2023). HST-dark galaxies make up the bulk of the sample at *z* > 4 while HST-faint ones dominate in number at ⟨*z*⟩ = 2.5, and there is an even split between them at ⟨*z*⟩ = 3.5. The average color, extinction, and stellar mass are indicated in the top-left corner of each plot. Overall color-selected galaxies are predominantly massive, dusty, or quiescent galaxies in agreement with the goal of the selection criteria in Figure 1. Furthermore, there is a clear downward trend in extinction and stellar mass ranging from ⟨*A*<sub>V</sub>⟩ ∼ 2.2 mag and ⟨log *M*<sub>*</sub>/M<sub>☉</sub>⟩ ∼ 10.5 at *z* < 3 to ⟨*A*<sub>V</sub>⟩ ∼ 1.2 mag and ⟨log *M*<sub>*</sub>/M<sub>☉</sub>⟩ = 10 at *z* > 4. The cutouts on the bottom show the 5″ × 5″ images of the MIPS and submillimeter/radio galaxies marked with letters in the UVJ diagram, also providing IDs and redshifts. From bottom to top, the images are color composites of three HST bands and three JWST bands, short and long wavelengths, respectively.</figcaption></figure>

Figure 6 demonstrates that the 2D-based, global dust attenuations for the galaxies in our sample do reproduce the dust content and emission observed in galaxies and are used for the construction of the SFMS. Even though the 5–100 Myr time-averaged SFRs fall below the MS, the fair comparison with observations behaves much better. On average, and down to log M*/M☉ ~ 9.7, our energy-budget-based SFRs reproduce the SFMS, with the median sSFR lying within 0.1 dex of the literature curves. However, the scatter of our SFMS is 2× larger than the typical 0.2–0.3 dex reported in the literature. Whether this behavior is intrinsic or related to observational uncertainties is beyond the scope of this paper.

In Figure 6 (and the right panel of Figure 5), we mark the galaxies in the redshift peak at z ~ 6 reported in the previous section (green triangles). All of these sources present very large SFRs for their stellar mass, well above the SFMS, consequently qualifying as starburst galaxies. The SFRs derived with the energy-budget calculation are also large, indicating the presence of significant amounts of dust. We remark here that the starburst galaxies seem to follow a separate sequence, very similar to the bimodality reported in Rinaldi et al. (2022). We cannot claim any robust interpretation for this effect due to the small number of galaxies we have in our sample, but our modeling points to systematic effects linked to age (rather, for example, than for extinction) since most of the galaxies located far away from the SFMS present very young bursts, many of them belonging to the XELG-z6 type. This would be consistent with the presence of strong emission lines, most probably [O III]+Hβ, in these starburst systems, similar to the sources reported in Endsley et al. (2022), Matthee et al. (2022), and

Laporte et al. (2022), but they also present relatively high dust contents, with average attenuations around A(V) = 2 mag. Indeed, Matthee et al. (2022) report not negligible amounts of dust in some of the (NIRCam-grism-selected) spectroscopically confirmed [O III] emitters.

We conclude that our 2D-SPS method based on JWST+HST data does provide robust dust estimations at spatial resolutions around 200 pc and that SFRs from SED fitting should be compared with SFRs from classical tracers with caution.

6.2.2. Rest-frame UVJ Colors

Figure 7 shows the overall distribution of all the CEERS galaxies (gray-scale density map) and the color-selected sample, divided into HST-faint and -dark in the rest-frame UVJ diagrams at different redshifts. The distribution of rest-frame colors is roughly consistent with what we see in the SFMS plot: The bulk of the sample is massive, dusty SFGs with a small fraction of quiescent galaxies starting at z ≳ 3. The majority of the galaxies identified as quiescent by sSFR (red squares) are located in the quiescent region of the UVJ diagram, showing that the two methods are quite consistent.

At 2 < *z* < 3, nearly all the color-selected galaxies are HST-faint, dusty and star-forming, with high average extinctions ⟨*A*<sub>V</sub>⟩ = 2.2 mag and relatively large stellar masses of ⟨log *M*<sub>*</sub>/M<sub>☉</sub>⟩ = 10.4. This is consistent with the color tracks in Figure 1 that indicate that only dusty galaxies are redder than the selection threshold at *z* < 3. That is, quiescent galaxies at these redshifts do not make the color cut. Similarly, since these galaxies are also relatively massive, they tend to be on the
brighter end of both the F356W and F150W mag, and consequently, the majority are HST-faint. Interestingly, many of these galaxies are located in the very high dust-obscuration region of the *UVJ* diagram with rest-frame colors exceeding $V-J > 2.0$ mag. Previous works have pointed out that this dusty corner of the *UVJ* diagram seems to exhibit a higher fraction of objects with late-type morphologies (or low Sérsic 1968 indices) and closer to edge-on inclinations that can increase the reddening along the line of sight due to dust within the disks, which lead to values of $A_V \gtrsim 2$ mag (e.g., Patel et al. 2012; Wang et al. 2018).

At $3 < z < 4$, there is roughly an even split between HST-dark and -faint galaxies. The average stellar mass is roughly the same, but the average extinction is much lower, $\langle A_V \rangle = 1.4$ mag, than at $z \sim 2.5$. Similarly, the average $V-J \sim 1.5$ mag is smaller and fewer galaxies are found in the dusty corner $V-J > 2$ mag. Notable exceptions are the submillimeter/radio galaxies with $V-J > 1.5$ mag. The lower average attenuation and $V-J$ color are likely driven by the inclusion of a larger number of quiescent galaxies ($\sim 17\%$), which enter the color selection at $z > 3$. We note also that at $z \sim 3.5$ some of the younger, perhaps recently quenched, quiescent galaxies with ages $\sim 1$ Gyr have bluer $U-V \sim 1.3$ mag colors, and some of them will not make the color threshold of the selection, F150W - F356W > 1.5 mag, even at $z = 3$ (orange template in Figure 1). To emphasize this point, in Figure 7 we show 6 out of the 17 galaxies in the recent Carnall et al. (2023) paper, which finds massive quiescent galaxies at $z > 3$ in the CEERS field. Three out of these six galaxies (green circles) are in our *UVJ* quiescent region but are not selected by our color criterion, as they are too blue in F150W - F356W. The other 11 galaxies with redder colors and higher redshifts are included in our sample and 9 of them are also classified as quiescent by our criteria.

At $4 < z < 7$, nearly all selected galaxies are HST-dark. The average extinction and color continue their downward trend relative to the lower-redshift bins with $\langle A_V \rangle = 1.2$ mag, $\langle V-J \rangle = 1.2$ mag, but there is also a noticeable decrease in the average stellar mass to $\langle \log M_*/M_{\odot} \rangle = 10$. Since stellar mass and dust are strongly correlated (e.g., Fang et al. 2018), a decrease in the average mass will also lower the attenuation. The lower average mass is caused by a decrease in the number of very massive (and dusty) galaxies (e.g., $\log M_*/M_{\odot} \gtrsim 10.8-11$) in our sample at redshift $z \gtrsim 4$ (see Figure 4). This is most likely caused by a combination of the steep decline in the number density of massive galaxies with redshift and our small surveyed area. For example, Alcalde Pampliega et al. (2019; see also Muzzin et al. 2013; Stefanon et al. 2015) find densities of $\log(N) \sim -4.75$ or $-5.00$ [Mpc$^{-3}$] at $3 < z < 6$, even after accounting for HST-dark galaxies. In the small area covered by CEERS, that density translates to less than $\sim 1$ galaxy (we find none), and even at $\log M_*/M_{\odot} > 10$ the number of galaxies is expected to be of the order of $\sim 20$ based on pre-JWST studies (e.g., extrapolating Muzzin et al. 2013, $\log(N) \sim -3.75$ [Mpc$^{-3}$] at $z \sim 3.5$). Consequently, more massive, dusty SFGs like those identified in submillimeter surveys (Wang et al. 2019) will likely be found in larger area surveys. In fact, two (of the five) submillimeter NOEMA galaxies at this redshift are among the most massive and higher-extinction galaxies. In particular, galaxy (a) was also recently discussed in Zavala et al. (2023) as

a source that might be misidentified with very high-z ($z \gtrsim 10$) candidates.

Aside from the average trends, the fraction of quiescent galaxies is a bit smaller compared to the previous redshift bin, $\sim 10\%$, but we note an increase in the overall scatter such that there are more galaxies with more extreme $V-J$ colors (e.g., $V-J > 2$ or $V-J < 0.25$). The majority of these galaxies are all at $z > 6$, in the redshift peak reported in previous sections. They also have peculiar SEDs, in terms of F356W - F410M color, and SEDs that are red at $\lambda > 2 \mu$m but get bluer at shorter wavelengths, similar to the SEDs of the galaxies presented in Labbe et al. (2022). At those redshifts, the $V-J$ color is not probed by the observed SED, and therefore, it is extrapolated from the best-fit SED template. Most of these galaxies are typically a combination of a dusty stellar population (hence the red $V-J$) plus a low mass, un-extinguished burst that causes the blue UV rest-frame color and a high, recent SFR (see SFH section). The true stellar properties of these galaxies are unclear and still debated, most remarkably whether the red color is caused predominantly by strong emission lines. Our 2D SED fits reveal strong (i.e., high equivalent width, EW) [O III] $\lambda4959,5007$ and H$\beta$ emission boosting the F356W flux (and sometimes the F140M emission), linked to a very young starburst on top of an older population. For this reason, we will identify these sources as extreme emission-like galaxies at $z \sim 6$ (XELG-z6, hereafter). They are a relatively numerous population; in fact, they dominate our sample at $z > 6$, and we might still be missing some, since their colors decrease with redshift, ranging from F150W - F356W = 3 mag down to our limit 1.5 mag (see postage stamps in Figure 8). In fact, the majority of galaxies from Labbe et al. (2022) at $z > 7$ are bluer and fainter at F356W, so we only have three in common with their sample. Given the selection of our sample, we are biased toward galaxies with emission lines entering the F356W passband, but if similar high-EW ELGs exist at lower redshifts, the lines could lie on other filters and affect the interpretation of JWST data. This is the case for one of the objects discussed in Zavala et al. (2023), for which $z = 4.6$ and $z = 16.7$ (Donnan et al. 2023) redshift solutions are obtained, depending strongly on the usage of templates with high-EW emission lines. In our 2D study, that exact galaxy, with ID nircam2-2159 in our catalog, is assigned with a redshift $z = 4.59$, with all its 74 pixels preferring the low-redshift solution over the $z \sim 16$ value. Apart from that, and interestingly, all XELGs-z6 share identical morphologies/appearances (see Section 6.3).

### 6.2.3. Selecting Subsets of Galaxies Based on Star Formation Activity

Based on Figures 5 to 7, in the following sections, we will divide our sample into three different galaxy types according to their star formation activity. We identify 25 galaxies as quiescent/dormant (QG) for which the SFMS plot provides $sSFR < 0.1$ Gyr$^{-1}$, 17 of them inside the quiescent wedge in the *UVJ* diagram. For the identification of QGs, we forced that the galaxies were below the quoted sSFR level for all the SFR timescales mentioned in Section 6.2.1 (ranging from 5 to 100 Myr), and the energy-budget sSFR. All other galaxies are considered as SFGs, except the special case of $z > 6$ sources, which we already introduced, XELG-z6. For QGs and SFGs, in some plots, we will also consider two redshift regimes, $z < 4$ and $z > 4$ for SFGs, and $z < 3.5$ and $z > 3.5$ for QGs, which divide the samples in roughly equal numbers (and help with the interpretation of the results).
brighter end of both the F356W and F150W mag, and consequently, the majority are HST-faint. Interestingly, many of these galaxies are located in the very high dust-obscuration region of the *UVJ* diagram with rest-frame colors exceeding $V-J > 2.0$ mag. Previous works have pointed out that this dusty corner of the *UVJ* diagram seems to exhibit a higher fraction of objects with late-type morphologies (or low Sérsic 1968 indices) and closer to edge-on inclinations that can increase the reddening along the line of sight due to dust within the disks, which lead to values of $A_V \gtrsim 2$ mag (e.g., Patel et al. 2012; Wang et al. 2018).

At $3 < z < 4$, there is roughly an even split between HST-dark and -faint galaxies. The average stellar mass is roughly the same, but the average extinction is much lower, $\langle A_V \rangle = 1.4$ mag, than at $z \sim 2.5$. Similarly, the average $V-J \sim 1.5$ mag is smaller and fewer galaxies are found in the dusty corner $V-J > 2$ mag. Notable exceptions are the submillimeter/radio galaxies with $V-J > 1.5$ mag. The lower average attenuation and $V-J$ color are likely driven by the inclusion of a larger number of quiescent galaxies ($\sim 17\%$), which enter the color selection at $z > 3$. We note also that at $z \sim 3.5$ some of the younger, perhaps recently quenched, quiescent galaxies with ages $\sim 1$ Gyr have bluer $U-V \sim 1.3$ mag colors, and some of them will not make the color threshold of the selection, F150W–F356W $> 1.5$ mag, even at $z = 3$ (orange template in Figure 1). To emphasize this point, in Figure 7 we show 6 out of the 17 galaxies in the recent Carnall et al. (2023) paper, which finds massive quiescent galaxies at $z > 3$ in the CEERS field. Three out of these 6 galaxies (green circles) are in our *UVJ* quiescent region but are not selected by our color criterion, as they are too blue in F150W–F356W. The other 11 galaxies with redder colors and higher redshifts are included in our sample and 9 of them are also classified as quiescent by our criteria.

At $4 < z < 7$, nearly all selected galaxies are HST-dark. The average extinction and color continue their downward trend relative to the lower-redshift bins with $\langle A_V \rangle = 1.2$ mag, $\langle V-J \rangle = 1.2$ mag, but there is also a noticeable decrease in the average stellar mass to $\langle \log M_*/M_{\odot} \rangle = 10$. Since stellar mass and dust are strongly correlated (e.g., Fang et al. 2018), a decrease in the average mass will also lower the attenuation. The lower average mass is caused by a decrease in the number of very massive (and dusty) galaxies (e.g., $\log M_*/M_{\odot} \gtrsim 10.8–11$) in our sample at redshift $z \gtrsim 4$ (see Figure 4). This is most likely caused by a combination of the steep decline in the number density of massive galaxies with redshift and our small surveyed area. For example, Alcalde Pampliega et al. (2019; see also Muzzin et al. 2013; Stefanon et al. 2015) find densities of $\log(N) \sim -4.75$ or $-5.00$ [Mpc$^{-3}$] at $3 < z < 6$, even after accounting for HST-dark galaxies. In the small area covered by CEERS, that density translates to less than $\sim 1$ galaxy (we find none), and even at $\log M_*/M_{\odot} > 10$ the number of galaxies is expected to be of the order of $\sim 20$ based on pre-JWST studies (e.g., extrapolating Muzzin et al. 2013, $\log(N) \sim -3.75$ [Mpc$^{-3}$] at $z \sim 3.5$). Consequently, more massive, dusty SFGs like those identified in submillimeter surveys (Wang et al. 2019) will likely be found in larger area surveys. In fact, two (of the five) submillimeter NOEMA galaxies at this redshift are among the most massive and higher-extinction galaxies. In particular, galaxy (a) was also recently discussed in Zavala et al. (2023) as

a source that might be misidentified with very high-$z$ ($z \gtrsim 10$) candidates.

Aside from the average trends, the fraction of quiescent galaxies is a bit smaller compared to the previous redshift bin, $\sim 10\%$, but we note an increase in the overall scatter such that there are more galaxies with more extreme $V-J$ colors (e.g., $V-J > 2$ or $V-J < 0.25$). The majority of these galaxies are all at $z > 6$, in the redshift peak reported in previous sections. They also have peculiar SEDs, in terms of F356W–F410M color, and SEDs that are red at $\lambda > 2 \mu$m but get bluer at shorter wavelengths, similar to the SEDs of the galaxies presented in Labbe et al. (2022). At those redshifts, the $V-J$ color is not probed by the observed SED, and therefore, it is extrapolated from the best-fit SED template. Most of these galaxies are typically a combination of a dusty stellar population (hence the red $V-J$) plus a low mass, un-extinguished burst that causes the blue UV rest-frame color and a high, recent SFR (see SFH section). The true stellar properties of these galaxies are unclear and still debated, most remarkably whether the red color is caused predominantly by strong emission lines. Our 2D SED fits reveal strong (i.e., high equivalent width, EW) [O III] $\lambda4959,5007$ and H$\beta$ emission boosting the F356W flux (and sometimes the F140M emission), linked to a very young starburst on top of an older population. For this reason, we will identify these sources as extreme emission-like galaxies at $z \sim 6$ (XELG-z6, hereafter). They are a relatively numerous population; in fact, they dominate our sample at $z > 6$, and we might still be missing some, since their colors decrease with redshift, ranging from F150W–F356W $= 3$ mag down to our limit 1.5 mag (see postage stamps in Figure 8). In fact, the majority of galaxies from Labbe et al. (2022) at $z > 7$ are bluer and fainter at F356W, so we only have three in common with their sample. Given the selection of our sample, we are biased toward galaxies with emission lines entering the F356W passband, but if similar high-EW ELGs exist at lower redshifts, the lines could lie on other filters and affect the interpretation of JWST data. This is the case for one of the objects discussed in Zavala et al. (2023), for which $z = 4.6$ and $z = 16.7$ (Donnan et al. 2023) redshift solutions are obtained, depending strongly on the usage of templates with high-EW emission lines. In our 2D study, that exact galaxy, with ID nircam2-2159 in our catalog, is assigned with a redshift $z = 4.59$, with all its 74 pixels preferring the low-redshift solution over the $z \sim 16$ value. Apart from that, and interestingly, all XELGs-z6 share identical morphologies/appearances (see Section 6.3).

### 6.2.3. Selecting Subsets of Galaxies Based on Star Formation Activity

Based on Figures 5 to 7, in the following sections, we will divide our sample into three different galaxy types according to their star formation activity. We identify 25 galaxies as quiescent/dormant (QG) for which the SFMS plot provides sSFR $< 0.1$ Gyr$^{-1}$, 17 of them inside the quiescent wedge in the *UVJ* diagram. For the identification of QGs, we forced that the galaxies were below the quoted sSFR level for all the SFR timescales mentioned in Section 6.2.1 (ranging from 5 to 100 Myr), and the energy-budget sSFR. All other galaxies are considered as SFGs, except the special case of $z > 6$ sources, which we already introduced, XELG-z6. For QGs and SFGs, in some plots, we will also consider two redshift regimes, $z < 4$ and $z > 4$ for SFGs, and $z < 3.5$ and $z > 3.5$ for QGs, which divide the samples in roughly equal numbers (and help with the interpretation of the results).
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Morphological properties of our sample of near-IR-faint, mid-IR-bright galaxies. On the left panel, we show the mass–size plot for HST-dark galaxies (hexagons) and HST-faint sources (circles), compared to the CEERS sample (gray-shaded density map). Sizes are represented by the effective radius of the Sérsic profile fit performed with the GALFIT software (Peng et al. 2002) on F356W images. MIPS24 emitters are marked with a cyan circle, submillimeter and radio galaxies with a blue star, quiescent galaxies by sSFR with a red square, and XELG-z6 with green triangles. We depict the scaling relationships for star-forming and quiescent galaxies at $z = 2.75$ (blue and red, respectively) from van der Wel et al. (2014). In the middle panels, histograms for the Sérsic index (top) and effective radius (bottom) are provided. We show estimations based on F356W images (light) as explained before, and on mass maps analyzed with the statmorph software (Rodriguez-Gomez et al. 2019). The middle-bottom panel also shows the distributions of Petrosian (1976) radii in mass maps. On the right panel, the mass–size relationship is again shown, this time based on mass map analyses performed with statmorph. Orange points show Petrosian radii. The blue and red lines show the ($r_{e,\text{mass}}$-based) mass–size relations for SFGs and quiescent galaxies from Suess et al. (2019). The cutouts on the bottom show the $5'' \times 5''$ images of XELG-z6 galaxies and other similar objects (not identified explicitly as such) at redshifts ranging from $5 < z < 7$, also providing IDs and redshifts. From bottom to top, the images are color composites of three HST bands and three JWST bands, short and long wavelengths, respectively.</figcaption></figure>

## 6.3. Structural Properties

### 6.3.1. Visual Appearances

Figures 5 and 7 show the cutouts of several interesting galaxies in our sample (HST-dark and -faint) at different redshifts. There are three cutouts per galaxy showing the color composite images based on three HST bands and three JWST bands, at short and long wavelengths. The purpose of these images is to not only illustrate the appearances of the galaxies as seen by JWST but also the substantial differences with respect to the HST images for galaxies that are faint and even dropouts in F160W.

At $2 < z < 3$, some of the very dusty galaxies at $V - J \gtrsim 2$ in Figure 7 do indeed look like nearly edge-on disks with signs of dust lanes visible in the NIRCam short-wavelength RGB images, while a smoother morphology is observed in the long-wavelength bands. Moreover, the comparison between the JWST and HST images for these galaxies reveals very striking differences, with the latter missing, in many cases, a significant fraction of the outskirts of the galaxies and, consequently, underestimating their true sizes (see the next section on parametric morphologies).

At higher redshift, $3 < z < 4$, the short-wavelength JWST bands start to probe the rest-frame UV of these (red) dusty and quiescent galaxies, and consequently, they can miss part of the structure in the outer parts of the galaxies. On the other hand, the higher spatial resolution of the short-wavelength bands

appears to identify one or more clumps in some galaxies that are not clearly seen in the smoother appearance of the long-wavelength JWST bands. The rest-frame UV dimming is even more extreme for the shallower HST bands, which, for some galaxies, can only detect a handful of pixels even in those previously identified in the CANDELS catalog.

Lastly, at $4 < z < 7$, galaxies become much fainter in the JWST short-wavelength bands, and in many cases, there is hardly anything to see on the images, with the exception of some quiescent galaxies that are still bright in F200W. In the long-wavelength bands, the galaxies are clearly detected but also faint and, in most cases, there are fewer distinct structural features (i.e., disks, or clear differences between the inner/outer parts of the galaxy) relative to what is seen at lower redshifts.

### 6.3.2. Morphological Parametric Measurements

Figure 8 shows the mass–size relation for the color-selected sample and the bulk of the CEERS sample (density map). The effective radii ($r_{e}$), measured along the major axis, and Sérsic (1968) indices ($n$) were determined from the F150W, F200W, and F356W images using GALFIT v3.0.5 (Peng et al. 2002). The code was run on the background-subtracted images, with the sky background held fixed at zero during fitting. Image cutouts that were fed to GALFIT have dimensions 2.5 times the Kron radius. The ERR array, which includes background sky, Poisson, and read noise, was used as the input noise map.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Morphological properties of our sample of near-IR-faint, mid-IR-bright galaxies. On the left panel, we show the mass–size plot for HST-dark galaxies (hexagons) and HST-faint sources (circles), compared to the CEERS sample (gray-shaded density map). Sizes are represented by the effective radius of the Sérsic profile fit performed with the GALFIT software (Peng et al. 2002) on F356W images. MIPS24 emitters are marked with a cyan circle, submillimeter and radio galaxies with a blue star, quiescent galaxies by sSFR with a red square, and XELG-z6 with green triangles. We depict the scaling relationships for star-forming and quiescent galaxies at $z = 2.75$ (blue and red, respectively) from van der Wel et al. (2014). In the middle panels, histograms for the Sérsic index (top) and effective radius (bottom) are provided. We show estimations based on F356W images (light) as explained before, and on mass maps analyzed with the statmorph software (Rodriguez-Gomez et al. 2019). The middle-bottom panel also shows the distributions of Petrosian (1976) radii in mass maps. On the right panel, the mass–size relationship is again shown, this time based on mass map analyses performed with statmorph. Orange points show Petrosian radii. The blue and red lines show the ($r_{e,\text{mass}}$-based) mass–size relations for SFGs and quiescent galaxies from Suess et al. (2019). The cutouts on the bottom show the $5'' \times 5''$ images of XELG-z6 galaxies and other similar objects (not identified explicitly as such) at redshifts ranging from $5 < z < 7$, also providing IDs and redshifts. From bottom to top, the images are color composites of three HST bands and three JWST bands, short and long wavelengths, respectively.</figcaption></figure>

## 6.3. Structural Properties

### 6.3.1. Visual Appearances

Figures 5 and 7 show the cutouts of several interesting galaxies in our sample (HST-dark and -faint) at different redshifts. There are three cutouts per galaxy showing the color composite images based on three HST bands and three JWST bands, at short and long wavelengths. The purpose of these images is to not only illustrate the appearances of the galaxies as seen by JWST but also the substantial differences with respect to the HST images for galaxies that are faint and even dropouts in F160W.

At $2 < z < 3$, some of the very dusty galaxies at $V - J \gtrsim 2$ in Figure 7 do indeed look like nearly edge-on disks with signs of dust lanes visible in the NIRCam short-wavelength RGB images, while a smoother morphology is observed in the long-wavelength bands. Moreover, the comparison between the JWST and HST images for these galaxies reveals very striking differences, with the latter missing, in many cases, a significant fraction of the outskirts of the galaxies and, consequently, underestimating their true sizes (see the next section on parametric morphologies).

At higher redshift, $3 < z < 4$, the short-wavelength JWST bands start to probe the rest-frame UV of these (red) dusty and quiescent galaxies, and consequently, they can miss part of the structure in the outer parts of the galaxies. On the other hand, the higher spatial resolution of the short-wavelength bands

appears to identify one or more clumps in some galaxies that are not clearly seen in the smoother appearance of the long-wavelength JWST bands. The rest-frame UV dimming is even more extreme for the shallower HST bands, which, for some galaxies, can only detect a handful of pixels even in those previously identified in the CANDELS catalog.

Lastly, at $4 < z < 7$, galaxies become much fainter in the JWST short-wavelength bands, and in many cases, there is hardly anything to see on the images, with the exception of some quiescent galaxies that are still bright in F200W. In the long-wavelength bands, the galaxies are clearly detected but also faint and, in most cases, there are fewer distinct structural features (i.e., disks, or clear differences between the inner/outer parts of the galaxy) relative to what is seen at lower redshifts.

### 6.3.2. Morphological Parametric Measurements

Figure 8 shows the mass–size relation for the color-selected sample and the bulk of the CEERS sample (density map). The effective radii ($r_{e}$), measured along the major axis, and Sérsic (1968) indices ($n$) were determined from the F150W, F200W, and F356W images using GALFIT v3.0.5 (Peng et al. 2002). The code was run on the background-subtracted images, with the sky background held fixed at zero during fitting. Image cutouts that were fed to GALFIT have dimensions 2.5 times the Kron radius. The ERR array, which includes background sky, Poisson, and read noise, was used as the input noise map.
Empirical PSFs were constructed using stars in all CEERS pointings. All galaxies in the image cutout within 3 mag of the primary source were fit simultaneously. All other sources were masked out during the fitting. Initial values for position, magnitude, effective radius, axis ratio, and position angle were set from the input photometry catalog. In addition, the following constraints were applied while fitting to keep values for neighboring sources within reasonable bounds: Sérsic index $0.2 \le n \le 8.0$, effective radius $0.3 \le r_e \le 400$ pixels, axis ratio $0.01 \le q \le 1$, magnitude $\pm 3$ mag from the initial value, and position $\pm 3$ pixels from the initial value. The radii in Figure 8 are based on the reddest band (F356W) since many of the HST-dark galaxies are quite faint in the bluer bands. The blue and red lines show the location of the mass-size scaling relations for star-forming and quiescent galaxies at $z=2.75$ based on HST/WFC3 data from van der Wel et al. (2014).

Overall, the bulk of the sample appears to follow the mass-size relation for SFGs, which is consistent with most of the sample being classified as star-forming by either the sSFR or rest-frame color criteria. More precisely, the color-selected sample has slightly lower average sizes than the star-forming relation at $z=2.75$, i.e., they follow a relation with a lower zero-point. We compare the F356W and F150W radii for the galaxies well detected in the bluest band (mostly HST-faint galaxies) and we find that the sizes are systematically smaller in F356W by $r_{e,F150W}/r_{e,F356W} = 1.7$. This agrees with the results of van der Wel et al. (2014), who also found smaller sizes at longer wavelengths when comparing HST/ACS and WFC3 measurements. If we account for this systematic effect by lowering the mass-size relation for SFGs at $z=2.75$ (blue line) by 0.23 dex (dashed gray line), it agrees well with the distribution of our sample. This is not surprising for the HST-faint galaxies with an average redshift of $\langle z \rangle = 3.2$, but it suggests that the HST-dark galaxies at higher redshifts, $\langle z \rangle = 4.2$, are not significantly smaller. That is, we do not find a strong size evolution for HST-dark galaxies and, in fact, many of them are quite large (see right panel of Figure 8 and discussion below), as recently pointed out by Nelson et al. (2022; blue circles in Figure 8). We also find that the axis-ratio ($q$) distribution peaks at $q \sim 0.6$ and ranges from $q = 0.2$ to 0.8, as expected for a sample of mainly disks galaxies, and we verify that the objects in Nelson et al. (2022) are indeed among those with the lowest axis-ratios.

As discussed in previous sections, the overall fraction of quiescent galaxies in the sample is relatively small (~15%). Among them, the most massive, $\log M_*/M_\odot > 10$, appear to follow relatively closely the mass-size relation of quiescent galaxies at $z=2.75$ whereas the least massive, $\log M_*/M_\odot < 10$, deviate from the relation toward larger masses. Nonetheless, most of these tend to have, on average, smaller sizes than SFGs of the same mass, as we would expect.

The most significant difference between the HST-faint and -dark samples is that the latter contains a population of galaxies with very small sizes $r_{e,F356W} < 0.25$ kpc and relatively small masses of $\log M_*/M_\odot \lesssim 9.5$. Roughly ~60% of these galaxies are XELG-z6, which reinforces the idea that they are indeed compact starbursts with strong emission lines coming from a small nuclear region. Similar galaxies in terms of compact morphology and rest-frame UV colors, but at $z=4-6$, are also found in our sample. Many of these galaxies have GALFIT fits with a Sérsic index of either $n > 3$ (spheroidal) or $n < 0.2$. The latter suggests that the size estimates are more uncertain, likely

because the small sizes and, in many cases fainter magnitudes, are close to the spatial-resolution limit of the F356W images. Nonetheless, the visual inspection of the residuals does not reveal any problems, and their visual appearances are indeed very small, featureless, and compact (see cutouts in Figure 8). Interestingly, three of these peculiar galaxies are quiescent by sSFR or UVJ colors and two more have low sSFRs and are adjacent to the quiescent region of the UVJ diagram (that is, post-starburst like or dormant, just starting a reignition).

The right panel of Figure 8 shows the mass-size relationship but this time obtained with statmorph (Rodriguez-Gomez et al. 2019) measurements performed on the mass maps. We provide several measurements of the size of our galaxies based on stellar-mass 2D distributions, namely, the Petrosian radius (which informs us of the total extent of the sources down to the surface brightness detection limit) and the effective radius (commonly used in this kind of mass-size plot). As mentioned in Section 4, most of our galaxies are quite extended, with a Petrosian radius extending from 1 to 10 kpc. The typical effective radius of our sample is 2 kpc, with HST-dark galaxies including a broader range of sizes and extending to smaller sizes than the HST-faint sample. The former follows a correlation that is steeper than what was obtained for $z \sim 2.5$ by Suess et al. (2019) for SFGs and very similar to that for quiescent systems.

Overall, the mass-size distribution of our sample in light and mass are quite similar. The most prominent differences are a larger scatter in the mass-based estimations and the larger sizes found in the mass maps for the very compact objects mentioned in the previous paragraph. Those galaxies have Sérsic indices below $n=1$. To aid in the discussion of light-to-mass morphological estimations, the middle panel of Figure 8 shows the distribution of Sérsic index values, which agrees with the results discussed above. First, we remark that the median (quartiles) ratio between effective and Petrosian radii is $1.8^{2.5}_{1.2}$. Second, the majority of the color-selected galaxies are disk-like, both for the HST-faint and -dark subsamples. The HST-dark galaxies have a tail of small Sérsic index values, which are precisely the very small galaxies mentioned above. The average Sérsic index of the HST-dark is marginally higher and extends to higher values than for the HST-faint galaxies. This is expected because most of the quiescent galaxies are in that subsample.

### 6.4. Star Formation Histories

Figure 9 shows the average global SFHs derived for the galaxies in our sample divided by star formation activity, namely, QGs, SFGs, and XELGs-z6. Global SFHs were computed by summing up the measurements from all pixels in a galaxy (each one fitted with a delayed exponentially decaying parameterization). This SFH was scaled up to reproduce the total stellar mass of the object by multiplying by the average (using all photometric bands) flux aperture correction obtained from the comparison of the sum of pixels arriving at our 2D stellar population analysis and the integrated-light elliptical aperture. After that, we transformed all SFHs to the common time frame of the age of the universe. Then, we averaged the functions in 100 Myr time intervals using all galaxies lying at a redshift below that corresponding to the given universe age. In order to make the average SFH curves more meaningful (in terms of combining coetaneous galaxies with similar natures), we divided the SFG and QG
Empirical PSFs were constructed using stars in all CEERS pointings. All galaxies in the image cutout within 3 mag of the primary source were fit simultaneously. All other sources were masked out during the fitting. Initial values for position, magnitude, effective radius, axis ratio, and position angle were set from the input photometry catalog. In addition, the following constraints were applied while fitting to keep values for neighboring sources within reasonable bounds: Sérsic index $0.2 \le n \le 8.0$, effective radius $0.3 \le r_e \le 400$ pixels, axis ratio $0.01 \le q \le 1$, magnitude $\pm 3$ mag from the initial value, and position $\pm 3$ pixels from the initial value. The radii in Figure 8 are based on the reddest band (F356W) since many of the HST-dark galaxies are quite faint in the bluer bands. The blue and red lines show the location of the mass-size scaling relations for star-forming and quiescent galaxies at $z=2.75$ based on HST/WFC3 data from van der Wel et al. (2014).

Overall, the bulk of the sample appears to follow the mass-size relation for SFGs, which is consistent with most of the sample being classified as star-forming by either the sSFR or rest-frame color criteria. More precisely, the color-selected sample has slightly lower average sizes than the star-forming relation at $z=2.75$, i.e., they follow a relation with a lower zero-point. We compare the F356W and F150W radii for the galaxies well detected in the bluest band (mostly HST-faint galaxies) and we find that the sizes are systematically smaller in F356W by $r_{e,F150W}/r_{e,F356W} = 1.7$. This agrees with the results of van der Wel et al. (2014), who also found smaller sizes at longer wavelengths when comparing HST/ACS and WFC3 measurements. If we account for this systematic effect by lowering the mass-size relation for SFGs at $z=2.75$ (blue line) by 0.23 dex (dashed gray line), it agrees well with the distribution of our sample. This is not surprising for the HST-faint galaxies with an average redshift of $\langle z \rangle = 3.2$, but it suggests that the HST-dark galaxies at higher redshifts, $\langle z \rangle = 4.2$, are not significantly smaller. That is, we do not find a strong size evolution for HST-dark galaxies and, in fact, many of them are quite large (see right panel of Figure 4 and discussion below), as recently pointed out by Nelson et al. (2022; blue circles in Figure 8). We also find that the axis-ratio ($q$) distribution peaks at $q \sim 0.6$ and ranges from $q = 0.2$ to 0.8, as expected for a sample of mainly disks galaxies, and we verify that the objects in Nelson et al. (2022) are indeed among those with the lowest axis-ratios.

As discussed in previous sections, the overall fraction of quiescent galaxies in the sample is relatively small (~15%). Among them, the most massive, $\log M_*/M_{\odot} > 10$, appear to follow relatively closely the mass-size relation of quiescent galaxies at $z=2.75$ whereas the least massive, $\log M_*/M_{\odot} < 10$, deviate from the relation toward larger masses. Nonetheless, most of these tend to have, on average, smaller sizes than SFGs of the same mass, as we would expect.

The most significant difference between the HST-faint and -dark samples is that the latter contains a population of galaxies with very small sizes $r_{e,F356W} < 0.25$ kpc and relatively small masses of $\log M_*/M_{\odot} \lesssim 9.5$. Roughly ~60% of these galaxies are XELG-z6, which reinforces the idea that they are indeed compact starbursts with strong emission lines coming from a small nuclear region. Similar galaxies in terms of compact morphology and rest-frame UV colors, but at $z=4-6$, are also found in our sample. Many of these galaxies have GALFIT fits with a Sérsic index of either $n > 3$ (spheroidal) or $n < 0.2$. The latter suggests that the size estimates are more uncertain, likely

because the small sizes and, in many cases fainter magnitudes, are close to the spatial-resolution limit of the F356W images. Nonetheless, the visual inspection of the residuals does not reveal any problems, and their visual appearances are indeed very small, featureless, and compact (see cutouts in Figure 8). Interestingly, three of these peculiar galaxies are quiescent by sSFR or UVJ colors and two more have low sSFRs and are adjacent to the quiescent region of the UVJ diagram (that is, post-starburst like or dormant, just starting a reignition).

The right panel of Figure 8 shows the mass-size relationship but this time obtained with statmorph (Rodriguez-Gomez et al. 2019) measurements performed on the mass maps. We provide several measurements of the size of our galaxies based on stellar-mass 2D distributions, namely, the Petrosian radius (which informs us of the total extent of the sources down to the surface brightness detection limit) and the effective radius (commonly used in this kind of mass-size plot). As mentioned in Section 4, most of our galaxies are quite extended, with a Petrosian radius extending from 1 to 10 kpc. The typical effective radius of our sample is 2 kpc, with HST-dark galaxies including a broader range of sizes and extending to smaller sizes than the HST-faint sample. The former follows a correlation that is steeper than what was obtained for $z \sim 2.5$ by Suess et al. (2019) for SFGs and very similar to that for quiescent systems.

Overall, the mass-size distribution of our sample in light and mass are quite similar. The most prominent differences are a larger scatter in the mass-based estimations and the larger sizes found in the mass maps for the very compact objects mentioned in the previous paragraph. Those galaxies have Sérsic indices below $n=1$. To aid in the discussion of light-to-mass morphological estimations, the middle panel of Figure 8 shows the distribution of Sérsic index values, which agrees with the results discussed above. First, we remark that the median (quartiles) ratio between effective and Petrosian radii is $1.8_{1.2}^{2.5}$. Second, the majority of the color-selected galaxies are disk-like, both for the HST-faint and -dark subsamples. The HST-dark galaxies have a tail of small Sérsic index values, which are precisely the very small galaxies mentioned above. The average Sérsic index of the HST-dark is marginally higher and extends to higher values than for the HST-faint galaxies. This is expected because most of the quiescent galaxies are in that subsample.

### 6.4. Star Formation Histories

Figure 9 shows the average global SFHs derived for the galaxies in our sample divided by star formation activity, namely, QGs, SFGs, and XELGs-z6. Global SFHs were computed by summing up the measurements from all pixels in a galaxy (each one fitted with a delayed exponentially decaying parameterization). This SFH was scaled up to reproduce the total stellar mass of the object by multiplying by the average (using all photometric bands) flux aperture correction obtained from the comparison of the sum of pixels arriving at our 2D stellar population analysis and the integrated-light elliptical aperture. After that, we transformed all SFHs to the common time frame of the age of the universe. Then, we averaged the functions in 100 Myr time intervals using all galaxies lying at a redshift below that corresponding to the given universe age. In order to make the average SFH curves more meaningful (in terms of combining coetaneous galaxies with similar natures), we divided the SFG and QG