# Massive Optically Dark Galaxies Unveiled by JWST Challenge Galaxy Formation
Models

M. Xiao1 P. A. Oesch1,2 D. Elbaz3 L. Bing4 E. Nelson5 A. Weibel1 R. P. Naidu6
E. Daddi3 R. J. Bouwens7 J. Matthee8 S. Wuyts9 J. Chisholm10 G. Brammer2 M.
Dickinson11 B. Magnelli3 L. Leroy3 P. van Dokkum12 D. Schaerer1 T. Herard-
Demanche7 L. Barrufet1 R. M. Endsley10 Y. Fudamoto13,14 C. Gómez-Guijarro3 R.
Gottumukkala1 G. D. Illingworth15 I. Labbe16 D. Magee15 D. Marchesini17 M.
Maseda18 Y. Qin19,20 N. Reddy21 A. Shapley22 I. Shivaei23 M. Shuntov2 M.
Stefanon24,25 K. Whitaker26,2 J. S. B. Wyithe19,20

###### Abstract

Over the past decade, the existence of a substantial population of optically
invisible, massive galaxies at $z\gtrsim 3$ has been implied from mid-infrared
to millimeter observations[Walter2012, Wang2019, Xiao2023]. With the
unprecedented sensitivity of the JWST, such extremely massive galaxy
candidates have immediately been identified even at $z>7$, in much larger
numbers than expected[Labbé2023, Barrufet2023, Adams2023, Pérez-González2023].
These discoveries raised a hot debate. If confirmed, early, high-mass galaxies
challenge the current models of galaxy formation[Boylan-Kolchin2023,
Lovell2023]. However, the lack of spectroscopic confirmations leads to
uncertain stellar mass ($M_{\star}$) estimates, and the possible presence of
active galactic nuclei (AGN) adds further uncertainty[Kocevski2023,
Matthee2023, Labbe2023b]. Here, we present the first sample of 36 dust-
obscured galaxies with robust spectroscopic redshifts at $z_{\rm spec}=5-9$
from the JWST FRESCO survey. The three most extreme sources at $z\sim 5-6$
($\sim$1 billion years after the Big Bang) are so massive
(log$M_{\star}/M_{\odot}$ $\gtrsim 11.0$) that they would require, on average,
about 50% of the baryons in their halos to be converted into stars – two to
three times higher than even the most efficient galaxies at later times. The
extended emission of these galaxies suggests limited contribution by AGN. This
population of ultra-massive galaxies accounts for 20% of the total cosmic star
formation rate density[Madau2014] at $z\sim 5-6$, suggesting a substantial
proportion of extremely efficient star formation in the early Universe.

Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix,
Switzerland

Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen,
Jagtvej 128, København N, DK-2200, Denmark

Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-
Yvette, France

Aix Marseille Université, CNRS, LAM, Laboratoire d’Astrophysique de Marseille,
Marseille, France

Department for Astrophysical and Planetary Science, University of Colorado,
Boulder, CO 80309, USA

MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts
Ave., Cambridge, MA 02139, USA

Leiden Observatory, Leiden University, NL-2300 RA Leiden, Netherlands

Department of Physics, ETH Zürich, Wolfgang-Pauli-Strasse 27, Zürich, 8093,
Switzerland

Department of Physics, University of Bath, Claverton Down, Bath, BA2 7AY, UK

Department of Astronomy, The University of Texas at Austin, 2515 Speedway,
Stop C1400, Austin, TX 78712-1205, USA

NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry
Ave., Tucson, AZ 85719, USA

Astronomy Department, Yale University, 52 Hillhouse Ave, New Haven, CT 06511,
USA

Waseda Research Institute for Science and Engineering, Faculty of Science and
Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan

National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo,
Japan

Department of Astronomy and Astrophysics, University of California, Santa
Cruz, CA 95064, USA

Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, Melbourne, VIC 3122, Australia

Department of Physics and Astronomy, Tufts University, 574 Boston Avenue,
Medford, MA 02155, USA

Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St.,
Madison, WI 53706 USA

School of Physics, University of Melbourne, Parkville, VIC 3010, Australia

ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D),
Australia

Department of Physics and Astronomy, University of California, Riverside, 900
University Avenue, Riverside, CA 92521, USA

Department of Physics & Astronomy, University of California, Los Angeles, 430
Portola Plaza, Los Angeles, CA 90095, USA

Steward Observatory, University of Arizona, Tucson, AZ 85721, USA

Departament d’Astronomia i Astrofìsica, Universitat de València, C. Dr.
Moliner 50, E-46100 Burjassot, València, Spain

Unidad Asociada CSIC “Grupo de Astrofísica Extragaláctica y Cosmología”
(Instituto de Física de Cantabria - Universitat de València)

Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA
✉e-mail<EMAIL_ADDRESS>
Recent studies have revealed a significant population of galaxies at $z\gtrsim
3$ that are completely missed or only tentatively detected in the deepest
optical observations, but brighter at longer wavelengths, such as Spitzer and
ALMA[Walter2012, Franco2018, Wang2019, Alcalde2019, Fudamoto2021, Xiao2023]
and more recently JWST[Barrufet2023, Gómez-Guijarro2023, Akins2023, Pérez-
González2023b, McKinney2023, Barro2023], the so-called optically dark/faint
galaxies. These galaxies are considered candidates of massive dusty star-
forming galaxies at high redshift[Wang2019], which are essential to constrain
the true star formation history, the number density, and formation mechanisms
of massive galaxies in the early Universe. However, the major challenge in
understanding the nature of these galaxies is the lack of systematic
spectroscopic observations.

Here, we exploit the unprecedented imaging and spectroscopic data provided by
the JWST FRESCO NIRCam/grism survey[Oesch2023]. Unlike slit spectroscopy,
which can only target pre-selected sources, grism spectroscopy can effectively
provide a complete sample of emission line galaxies, with the FRESCO survey
reaching down to $\sim 2\times 10^{-18}{\rm ergs}^{-1}{\rm cm}^{-2}$
($5\sigma$ depth) over a continuous region of $2\times 62$ arcmin2 in the
GOODS-South and North fields[Giavalisco2004]. This unique advantage of grism
spectroscopy can help to explore missed galaxy populations, e.g., optically
dark galaxies. With the deep 4-5$\mu$m grism spectra in the F444W filter, we
probe H${\alpha}$+[NII] emission lines at $z=4.86-6.69$ and
[OIII]$\lambda\lambda 4960,5008$+H${\beta}$ at $z=6.69-9.08$. The FRESCO
survey also provides high-resolution NIRCam images (F182M, F210M, and F444W;
typical 5$\sigma$ depth of $\sim$28.2 mag), that complement ancillary data
from both HST[Whitaker2019] and now also JWST/JADES[Rieke2023].

Using the multi-wavelength imaging and grism spectroscopy, we select optically
dark/faint galaxies that are red (F182M $-$ F444W $>$ 1.5 mag) and have
clearly detected emission lines of at least $>8\sigma$ for robust redshift
identification. The color cut together with emission line detections are
designed to select star-forming galaxies with significant Balmer breaks, i.e.,
massive and dusty star-forming galaxies. We identify 36 galaxies at $z_{\rm
spec}>5$ in our parent sample (Extended Data Fig. 1, Methods). Thanks to the
JWST spectroscopy, we can now obtain more accurate stellar masses by
accounting for the true emission line contributions in photometric bands in
the SED fit at a fixed $z_{\rm spec}$. The $M_{\star}$ are then determined by
fitting the UV-to-NIR SED using multiple codes (Bagpipes[Carnall2018],
CIGALE[Boquien2019]) and with different star formation histories, which
produce consistent stellar mass values within the error range (see Methods for
details). The redshift and stellar mass estimates are listed in Extended Data
Table. 1. The parent sample of spectroscopically confirmed, optically
dark/faint galaxies spans a redshift range of $z_{\rm spec}=5-9$, with stellar
masses covering a wide distribution of log($M_{\rm\star}$/$M_{\odot}$)
$=8.3-11.4$. This confirms the high-redshift, dusty, and high-mass nature of
this elusive class of galaxies that has so far only been studied from
photometry alone[Wang2019, Xiao2023, Barrufet2023, Labbé2023, Akins2023,
Pérez-González2023b, McKinney2023].

Among the parent sample, three optically dark galaxies at $z_{\rm spec}\sim
5-6$ show the most extreme properties, i.e., extremely high masses with
log($M_{\star}/M_{\odot}$) $\gtrsim 11.0$ (Chabrier initial mass function,
IMF[Chabrier2003]; Figs. 1, 2). These sources are named S1, S2, and S3. One of
these has been known previously as a sub-millimeter galaxy GN10[Riechers2020].
Their extreme mass makes metals and dust inefficient to escape through stellar
feedback from their ultra-deep gravitational potential, which is why they have
been missed by optical observations until now. The detection of these
optically dark galaxies shows another aspect of the power of JWST: it probes
not only the very distant, but also the highly dusty and hidden Universe.
Based on the deep JWST observations, we find that the three galaxies have
extremely red properties (F182M $-$ F444W $>$ 3.5 mag; Fig. 1, Extended Data
Fig. 1), due to their significant dust attenuation ($A_{\rm V}>3$ mag). Self-
consistently, we find that they also have highly dust-obscured star formation
rates (SFRIR), which are obtained from far-infrared SED fits with CIGALE,
specifically $625\pm 31M_{\odot}$ yr-1 for S1, 1030${}_{-150}^{+190}$ M⊙ yr-1
for S2, and $1137\pm 31$ M⊙ yr-1 for S3 (see Methods). This indicates that
they are in the process of very efficient stellar mass build-up. We find no
signs of a significant contribution to the rest-frame optical light by an AGN
based on our investigation of the emission lines, source morphology, and
multi-wavelength data. Therefore, we conclude that the ultra-massive nature of
these three galaxies is reliable.

As shown in Fig. 2, we obtain robust spectroscopic redshift measurements at
$z_{\rm spec}=5-9$ from the JWST/FRESCO NIRCam grism survey to characterize
the population of 36 optically dark/faint galaxies. We highlight the three
ultra-massive galaxies at $z_{\rm spec}\sim 5.5$ and further reveal the
extremity of their high stellar masses at their redshifts in comparison with
the most massive galaxy which can be found according to the halo mass function
and cosmic baryon fraction[Boylan-Kolchin2023, Lovell2023]. Under this
paradigm, we derive the most massive dark matter halo mass at different
redshifts in the corresponding FRESCO survey volume ($\sim 1.2\times 10^{6}$
Mpc3 at $z=5-9$) according to the halo mass function. The maximum stellar mass
is inferred from the maximum dark matter halo mass, based on
$M_{\star}^{\mathrm{max}}=\epsilon f_{\rm b}M_{\mathrm{halo}}^{\mathrm{max}}$,
with a cosmic baryon fraction $f_{\rm b}=\Omega_{\rm b}/\Omega_{\rm
m}=0.158$[Planck2020], and the maximum theoretical efficiency ($\epsilon$) of
converting baryons into stars. Here, we consider two possible cases of
$\epsilon$, as shown in Fig. 2: the highest efficiency from observation-based
phenomenological modelings, such as abundance matching and halo occupation
distribution models ($\epsilon_{\rm max,obs}=0.2$; dashed line)[Moster2013,
Wechsler2018, Shuntov2022] and the maximum efficiency logically allowed
($\epsilon=1$; solid line).

By comparing the masses of our galaxies with the predictions in Fig. 2, it is
clear that these sources require an extremely efficient conversion of baryons
to stars of about 0.5 on average – two to three times the highest efficiency
observed at lower redshift ($\epsilon_{\rm max,obs}\simeq 0.2$). These three
galaxies lie at $z_{\rm spec}\sim 5.5$, demonstrating that the existence of
ultra-massive galaxies that challenge our galaxy formation models is not
restricted to the most distant Universe[Labbé2023], but extends over a wide
period, which was previously hidden by heavy dust obscuration. This is also
demonstrated by two further galaxies in our sample at $z_{\rm spec}>7$ that
imply $\epsilon>0.2$, similar to the photometrically selected sources from
previous studies[Labbé2023]. However, their stellar masses are less robust,
since the determination of stellar masses is mostly constrained by the F444W
band. At $z\sim 7-9$, this band traces the emission of the rest-frame B-band,
which is sensitive to dust attenuation and stellar age. Therefore, we mainly
focus on the three ultra-massive galaxies at $z_{\rm spec}\sim 5.5$. At these
lower redshifts, the determination of the stellar mass and SFR is more robust
because F444W traces the rest-frame optical.

The high required baryon-to-stellar conversion efficiency in the three ultra-
massive galaxies makes them extremely rare in simple galaxy formation models.
Conservatively assuming a maximum observed efficiency of $\epsilon_{\rm
max,obs}=0.2$[Moster2013, Wechsler2018, Shuntov2022], the stellar masses of
S1, S2, and S3 correspond to dark matter halo masses of log($M_{\rm
halo}/M_{\odot}$) $=12.88^{+0.11}_{-0.13}$, $12.68^{+0.23}_{-0.17}$, and
$12.54^{+0.17}_{-0.18}$, respectively. The observed volume density is $\sim
3.0\times 10^{-6}$ Mpc-3 at $z\sim 5-6$ in the 124 arcmin2 of the FRESCO
survey fields. For the most extreme case S1, compared to the theoretical
cumulative halo number densities of $\sim 2.8\times 10^{-9}$ Mpc-3, the
probability of detecting such a source in a random field as large as the
FRESCO survey is only 0.0009 (Fig. 3). It means that if the galaxy
distribution in the Universe is homogeneous, and if $\epsilon=0.2$, we would
only have expected to detect one source such as S1 in a field 1,100 times the
area of FRESCO. For S2 and S3, the probabilities are 0.02 and 0.09,
respectively, which are also low enough to require 50 and 10 times the FRESCO
field to detect them, respectively. This shows that the star-formation
efficiencies in these galaxies must be significantly higher than normally
found at lower redshifts in standard galaxy formation models within cold dark
matter halos[White1978, Boylan-Kolchin2023].

Thanks to the combined power of ALMA and JWST, we now have the first insight
into the physical mechanism at play in massive galaxy build-up (right panel of
Fig. 3). In S1, we observe dust emission (half-light radius $R_{\rm
e,1.1mm}=2.88\pm 0.91$ kpc) significantly more extended than the stellar
component ($R_{\rm e,F444W}=0.76\pm 0.01$ kpc), giving access to the
circumgalactic region. We note that the dust distribution in dusty star-
forming galaxies is typically more compact relative to the stellar
distribution[Tadaki2020, Puglisi2021, Gómez-Guijarro2022], which is thought to
be caused by angular momentum loss driven by mergers or counter-rotating gas
accretion. Our findings for S1 seem to favor an opposite scenario, in which
the galaxy activity may be triggered by a strong event of corotating and
coplanar gas accretion through the cold streams[Danovich2012, Danovich2015,
Kretschmer2022, Dekel2020, Xiao2022]. In addition, we find a 6$\sigma$ off-
nuclear clump in the dust distribution map of S1 (see the right panel of Fig.
3). Surprisingly, no stellar emission associated with this clump is detected
in the deep JWST data. It suggests that S1 might be undergoing a merger with
an even more obscured galaxy that is still missed, even by JWST. Future ALMA
data of the [CII] emission line will be able to test this scenario.

Extremely massive galaxies contribute significantly to the total cosmic star-
formation rate density in the early Universe. We derive this in the redshift
bin for our H${\alpha}$ sample ($z=4.9-6.6$). Including only the three ultra-
massive galaxies ($\epsilon>0.2$), the SFR density reaches 2.4 $\times$ 10-3
$M_{\odot}$ yr-1Mpc-3 at $z$ $\sim$ 5.8. Compared to the total cosmic star-
formation history[Madau2014] (Fig. 4), at $z\sim 5.8$, the existence of these
galaxies implies a surprisingly high efficiency of baryon conversion into
stars: about 20% of the stars formed above the peak value of $\epsilon_{\rm
max,obs}=0.2$. This finding suggests a substantial proportion of extremely
efficient star formation in the early Universe.

With the secure spectroscopic redshift and stellar mass measurements, our
results provide strong evidence that the early Universe has to be two to three
times more efficient in forming massive galaxies than the average trend found
by previous studies at later times. Our discovery, together with the possible
excess of UV-luminous galaxies at $z>8$ revealed by JWST observations,
indicates that early galaxy formation models may need to be revised. A top-
heavy IMF[Zhang2018] or feedback-free starburst in the presence of extremely
efficient gas infall from cold filaments[Dekel2023] could possibly reproduce
the observed extreme properties of our sample. However, these scenarios still
remain to be investigated with more detailed observations. With higher spatial
resolution and (or) higher sensitivity, future ALMA/NOEMA and deep JWST
spectroscopic observations could help consolidate the massive nature of these
galaxies with dynamical mass measurements and test the different scenarios for
their formation with the analysis of the kinematics and chemical composition
of the interstellar medium.

Figure 1: Images and Spectra of the three optically dark, ultra-massive
galaxies from the JWST FRESCO survey. From left to right: 4′′ $\times$ 4′′
stamps obtained in JWST/NIRCam filters (1.82$\mu$m, 2.10$\mu$m, and
4.44$\mu$m), RGB images (F182M in blue, F210M in green, and F444W in red),
H$\alpha$ line map, and 1D spectra (covering H$\alpha$, [NII], and [SII]
emission lines) obtained from NIRCam/grism observations with the F444W filter.
The black arrow in the lower-right corner of H$\alpha$ line map shows the
dispersion direction of the F444W grism.  Figure 2: Stellar masses of the
three ultra-massive and optically dark galaxies (red-filled large circles)
compared to model expectations. The red-filled circles show our whole parent
sample of 36 spectroscopically-confirmed optically dark/faint galaxies at
$z_{\rm spec}>5$. Unlike some massive sources reported in the
literature[Labbé2023] that have only $z_{\rm phot}$ (grey points), our
galaxies from the JWST FRESCO survey have $z_{\rm spec}$, as well as line-
corrected photometry, and thus more robust $M_{\star}$. The grey-shaded region
indicates the stellar mass prohibited by the standard $\Lambda$CDM cosmology
[Planck2020], which is calculated from the maximum halo mass
($M_{\mathrm{halo}}^{\mathrm{max}}$) that can be observed in the FRESCO survey
volume, based on $M_{\star}^{\mathrm{max}}=\epsilon f_{\rm
b}M_{\mathrm{halo}}^{\mathrm{max}}$, with a cosmic baryon fraction $f_{\rm
b}=\Omega_{\rm b}/\Omega_{\rm m}=0.158$, and the maximum theoretical
efficiency of converting baryons into stars, $\epsilon=1$. The grey dashed
line indicates the maximum efficiency of $\epsilon_{\rm max,obs}\sim 0.2$ from
empirical modeling methods based on observational data at lower
redshift[Moster2013, Wechsler2018, Shuntov2022]. The three most massive
galaxies in our sample require significantly higher star-formation
efficiencies of $\epsilon\sim 0.5$ on average.  Figure 3: The ultra-massive
sources S1, S2 and S3 with log(M⋆/M⊙) $>11.0$ at $z=5-6$. Left: Cumulative
comoving number density of dark matter halos as a function of halo mass at
different redshifts. We show the halo mass of S1, S2, and S3 assuming a
maximum observed efficiency of $\epsilon_{\rm max,obs}=0.2$. The observed
number density derived from the FRESCO survey is shown as red points. The
theoretical number density of S1, S2, and S3 in the same halo mass from the
standard $\Lambda$CDM cosmology (black points) is $\sim$1,100 times, 50 times,
and 10 times lower than the observations, respectively. Right: The most
extreme source S1 with log(M⋆/M⊙) $=11.37^{+0.11}_{-0.13}$ at $z=5.58$
exhibits an extended dust disk. We show the ALMA 1.1 mm dust contours overlaid
on the JWST RGB image. The 1.1 mm image (4′′ $\times$ 4′′) is shown in the top
right corner, and the beam size (0447 $\times$ 0418) is displayed in the lower
left corner. The dust is clearly much more extended (half-light radius $R_{\rm
e}=2.88\pm 0.91$ kpc) compared to the stellar distribution from JWST. It also
shows a 6$\sigma$ clump at 1.1mm in the off-nuclear region, which,
surprisingly, is undetected even in the deep, long wavelength JWST data.
Figure 4: Contribution of ultra-massive optically dark galaxies to the cosmic
SFR density. The black line is the computation of the total cosmic SFR density
(SFRD), $\psi$, corrected for dust attenuation, as a function of redshift by
Madau & Dickinson (MD14)[Madau2014]. At $z\gtrsim 4$, the curve is primarily
based on UV-selected galaxies (Lyman-break galaxies, i.e., LBGs; grey open
triangles)[Bouwens2012a, Bouwens2012b]. Red stars are spectroscopically-
confirmed ultra-massive optically dark galaxies from the FRESCO survey with
$\epsilon>0.2$. The existence of these galaxies implies that up to 20% of all
star-formation at $z\gtrsim 5$ occurs with an extremely high baryon-to-star-
formation efficiency.

## References

Methods

##  .1 Cosmology.

Throughout this paper, we assume a Planck cosmology[Planck2020] and adopt a
Chabrier IMF[Chabrier2003] to estimate stellar masses (M∗) and star formation
rates (SFR). When necessary, data from the literature have been converted with
a conversion factor of M∗ (Salpeter IMF)[Salpeter1955] = 1.7 $\times$ M∗
(Chabrier IMF). All magnitudes are in the AB system[Oke1983], such that
$m_{\rm AB}=23.9-2.5$ $\times$ log(Sν [$\mu$Jy]).

##  .2 FRESCO NIRCam grism spectra and imaging data.

The NIRCam data used in this paper stem from the JWST FRESCO
survey111https://jwst-fresco.astro.unige.ch (GO-1895; PI: P. Oesch;
see[Oesch2023] for details). FRESCO obtained direct images in three filters
(F182M, F210M, and F444W) as well as NIRCam/grism spectroscopy in the F444W
filter over $\sim$60 arcmin2 in each GOODS field, North and South, through two
2$\times$4 NIRCam/grism mosaics. The grism spectra span a wavelength of 3.8 to
5.0 $\mu$m, with some parts of the mosaic having slightly reduced coverage.
With an exposure time of 7043 s, the grism data reach an average 5 $\sigma$
line sensitivity of 2$\times$10-18 erg s-1 cm-2 at a resolution of
R$\sim$1600\. The FRESCO data were acquired between November 2022 and February
2023.

We use the publicly available grizli
pipeline222https://github.com/gbrammer/grizli to reduce the slitless
NIRCam/grism data (see also Brammer et al, in prep). The raw data are obtained
from the MAST archive, and calibrated with the standard JWST pipeline, before
being aligned to a Gaia-matched reference frame.

In order to produce a line-only dataset we remove the continuum along the
dispersion direction following[Kashino2022]. We subtract a running median
filter with a central 12-pixel gap (corresponding to $\sim 20\AA$ rest-frame)
along each row of the grism images. To minimize the self-subtraction of bright
lines, the filtering is done in two passes, masking pixels with significant
flux after the first pass.

Based on the direct F444W image an individual kernel is created for each
source in order to perform optimal extractions of 1D spectra. Slightly
modified sensitivity curves and spectral traces are used taking the publicly
available v4 grism configuration files333https://github.com/npirzkal/GRISMCONF
as the starting point.

In parallel to the long-wavelength grism data, FRESCO obtained short-
wavelength imaging in F182M and F210M. At the end of each grism exposure,
direct images in F444W are obtained to ensure that the grism data can be
aligned. The 5 $\sigma$ depths of the direct images in the F182M, F210M, and
F444W filters are 28.3, 28.1, and 28.2 mag, respectively, as measured in
circular apertures of 032 diameter.

##  .3 Ancillary data and multi-wavelength catalogs.

In addition to the FRESCO NIRCam data, we also make use of public HST and JWST
imaging over the GOODS fields. Most importantly, this includes the HST/ACS and
WFC3/IR data from the original GOODS survey[Giavalisco1996] as well as
CANDELS[Koekemoer2011, Grogin2011]. For a full list of all HST programs
covering this field see the Hubble Legacy Field (HLF) release
page444https://archive.stsci.edu/prepds/hlf/ (see also[Whitaker2019]
and[Illingworth2016]). In GOODS-South, our data cover the Hubble Ultra Deep
Field (HUDF[Beckwith2006]) and we include the deep JADES NIRCam imaging that
was released in June 2023[Eisenstein2023, Rieke2023]. The images are all co-
aligned and drizzled to a common 40 mas/pixel frame.

Multi-wavelength catalogs are derived using SExtractor[Bertin1996] in dual
image mode, using an inverse-variance weighted stack of the NIRCam long-
wavelength wide filters as the detection image. Fluxes are measured in 016
radius circular apertures in images that are PSF-matched to the F444W band.
Total fluxes are derived from the Kron AUTO aperture provided by the
SExtractor in the F444W band, in addition to a correction based on the
encircled energy of the Kron aperture on the F444W PSF.

Photometric redshifts are estimated through the SED-fitting code
EAZY[Brammer2008], using the blue_sfhz_13 template set which imposes redshift-
dependent SFHs, disfavoring SFHs that start earlier than the age of the
Universe at a given redshift555https://github.com/gbrammer/eazy-
photoz/tree/master/templates/sfhz. We apply an error floor of 5% prior to
running eazy to account for possible remaining systematic uncertainties in the
photometric fluxes[Labbé2023] and to allow for more flexibility in the SED-
fitting. The methods used to derive the multi-wavelength catalogs with
photometric redshifts will be presented in more detail in Weibel et al., in
prep.

Far-infrared and millimeter data. In the GOODS-South field, the FRESCO survey
overlaps with about 80% of the GOODS-ALMA 1.1mm survey[Franco2018, Gómez-
Guijarro2022, Xiao2023] (PI: D. Elbaz). The GOODS-ALMA survey covers a
continuous area of 72.42 arcmin2 and has a combined dataset[Xiao2023] of high-
and low-resolution 1.1mm observations obtained from ALMA Cycle 3 and Cycle 5.
The combined map achieves an rms sensitivity of $\sigma$ $\simeq$ 68.4 $\mu$Jy
beam-1 with a spatial resolution of 0447 $\times$ 0418\. The GOODS-ALMA 2.0
source catalog is presented in ref.[Gómez-Guijarro2022]. In addition, the
field is also covered by the JCMT/SCUBA-2 at 850$\mu$m[Cowie2018]. The source
S1, from our sample in GOODS-S, also has been observed by ALMA band-7
observations as a follow-up to a JCMT/SCUBA-2 target[Cowie2018]. In the GOODS-
North field, the FRESCO coverage overlaps with the 450$\mu$m and 850$\mu$m
surveys of JCMT/SCUBA-2[Cowie2017, Barger2022], as well as the 1.2mm and 2mm
observations of the NIKA2 Cosmological Legacy Survey[Bing2023] and their NOEMA
follow-up observations (PI: L. Bing, private communication).

##  .4 Optically dark/faint galaxies at $z_{\rm spec}\gtrsim 5$.

The traditional approach used to systematically select optically dark/faint
galaxies is based on photometric data only (pre- and present-JWST era). They
are required to be faint or completely undetected at optical wavelengths, but
bright in the infrared (e.g., $H>$ 27 mag & [4.5] $<$ 24 mag for
$H$-dropouts[Wang2019]; $H>$ 26.5 mag & [4.5] $<$ 25 mag for OFGs[Xiao2023];
see also[Alcalde2019, Barrufet2023, Gómez-Guijarro2023]). These magnitude
and/or color cuts are designed to select red galaxies with strong Balmer or
4000$\AA$ breaks at $z\gtrsim 3$, so they should be either quiescent/passive
galaxies or dusty star-forming galaxies with significant dust attenuation. The
$H$-band magnitude cut helps to avoid significant contamination from passive
galaxies and low-$z$ sources[Xiao2023], but obviously, they cannot be entirely
ruled out.

Now, by taking advantage of the unprecedented imaging and spectroscopic data
provided by the JWST FRESCO survey, we no longer need to use the magnitude
cuts described above to select optically dark/faint galaxies. The FRESCO grism
spectra help us more securely identify high-$z$ star-forming sources through
strong nebular emission lines. The spectral coverage of the FRESCO survey with
F444W filter allows detections of either H${\alpha}$+[NII]+[SII] lines or
[OIII]$\lambda\lambda 4960,5008$+H${\beta}$ lines at $z=4.9-9.0$ continuously.

We select candidate galaxies that are red (F182M $-$ F444W $>$ 1.5 mag;
Extended Data Fig. 1) and have strong [OIII]$\lambda\lambda
4960,5008$+H${\beta}$ or H${\alpha}$+[NII]+[SII] emission line detections. The
color threshold we set is similar to the lower color limit of optical
dark/faint galaxy selection criteria in the literature[Xiao2023, Gómez-
Guijarro2023, Akins2023, Pérez-González2023b, Barro2023], which helps to
further compare our sample with those of the pre-JWST/spectra sample. We
require all the galaxies to have $>$5$\sigma$ detections at F444W. With some
galaxies not detected at F182M, we use the 2$\sigma$ limit as an upper limit
in the color selection. For the emission line detection, we require the
strongest line to be $>8\sigma$ and the second strongest line to be
$>3\sigma$. For some galaxies with only one line detected, we require that the
line has $>8\sigma$ detection and their $z_{\rm spec}$ must be within the
16-84th percentile uncertainties of the $z_{\rm phot}$ from the UV-to-NIR SED
fitting with EAZY code. Finally, we have 36 galaxies at $z_{\rm spec}\gtrsim
5$ in our sample. Among these, 22 out of 36 galaxies have at least two
emission line detections or $z_{\rm spec}$ measurements confirmed from the
literature (see Extended Data Table. 1).

The Extended Data Fig. 1 presents our sample selection results on the full
FRESCO field, as well as the comparison with previous selection methods and
model expectations. Compared to previous selections of $H$-dropouts from
Spitzer/IRAC[Wang2019], our method extends the identification of optically
dark/faint galaxies to fainter F444W/4.5$\mu$m fluxes and bluer colors.

In Extended data Fig. 1, we present a comparison between the color and
magnitude of our sample and the model-predicted SED of galaxies. The redshift
evolution track of galaxies across the color-magnitude diagram is computed
using BC03[Bruzual2003] stellar population synthesis model with constant star
formation history starting at $z=10$. The SEDs are further reddened using a
Calzetti dust extinction law with E(B-V) at 0.2, 0.4, 0.6, and 0.8 to simulate
the color and magnitude evolution of dusty star-forming galaxies. The
evolution tracks of galaxies with a total stellar mass of $10^{10}M_{\odot}$
and $10^{11}M_{\odot}$ are presented.

While the optically faint selection criteria are defined to target dust-
obscured sources, galaxies can also exhibit red F182M $-$ F444W colors due to
emission line contamination in the F444W filter. Thanks to the FRESCO emission
line measurements, we can now correct this and check which sources only appear
red due to line emission. Therefore, we subtract the measured line fluxes from
the broad-band photometry to derive pure continuum flux measurements (F444Wc).
This shows that a small number ($\sim$10 sources) at the faint end only
satisfies our selection criteria due to line boosting (see Fig. 1). The
brighter galaxies discussed in the main text, however, are only a little
affected by this.

The emission-line corrected colors further ensure a fair comparison between
our sample with the evolution tracks. This clearly reveals that our JWST
imaging+spectroscopic selection could (1) select optically dark/faint galaxies
with high purity and without contamination of passive galaxies, and (2) extend
the selection of dusty galaxies to lower $M_{\star}$, lower amount of (but
still significant) dust reddening and higher $z_{\rm spec}$.

Our stellar mass estimation suggests the FRESCO dusty star-forming galaxy
sample has log($M_{\rm\star}/M_{\odot}$) $=8.3-11.4$ and $A_{V}=0.1-3.7$ at
$z_{\rm spec}>5$ (see Extended Data Table 1 and the section on stellar mass
measurements for further details). For comparison, the traditional methods of
selecting dusty galaxies (grey triangular regions[Wang2019] in Extended Data
Fig. 1; see also[Xiao2023, Alcalde2019, Barrufet2023, Gómez-Guijarro2023])
miss the majority of high-$z$ and(or) faint sources in our sample, but still
exclusively select the three ultra-massive and most reddened galaxies with
log($M_{\star}/M_{\odot}$) $\gtrsim 11.0$ at $z_{\rm spec}\sim 5.5$. This
further suggests the power and importance of JWST data in providing a complete
picture of dusty galaxies in the early Universe.

Stellar masses. The physical properties of these 36 optically dark/faint
galaxies in our parent sample are estimated by fitting the UV-to-NIR SED from
the JWST+HST photometry using BAGPIPES[Carnall2018], with the $z_{\rm spec}$
fixed. We assume a constant star formation history (SFH), the Bruzual stellar
population models[Bruzual2003], and a Calzetti dust attenuation
law[Calzetti2000]. We adopt a broad metallicity grid from 0.1 to 2.0
$Z_{\odot}$, an $A_{\rm V}$ grid from 0 to 5 mag, and an age grid from 10 Myr
to 1.5 Gyr. For our SED fits, we use the emission line subtracted F444W fluxes
in order to remove one free parameter of the emission line model and only fit
the stellar continuum emission. We have also ensured that no other strong
emission lines are contaminating other photometric bands. From this, we obtain
the main physical properties of the 36 optically dark/faint galaxies:
$M_{\star}$, $A_{\rm V}$, and SFR. These properties of our sample are listed
in Extended Data Table 1. Their locations in the star-formation main sequence
(SFMS) are shown in Extended Data Fig. 2. We would like to emphasize that due
to the lack of far-infrared and millimeter data for most of the sources, the
physical properties here are only obtained from UV-to-NIR SED fits. In this
case, the $M_{\star}$ and SFRs of these dusty star-forming galaxies (and the
SFR densities in Fig. 4) might be underestimated, as they may contain hidden
dust regions that absorb all the UV photons, which cannot be reproduced with
dust extinction corrections[Xiao2023, Elbaz2018, Puglisi2017]. However, the
possible underestimation of $M_{\star}$, SFR, and therefore the SFR density,
would make our main results in this paper even more significant. A further
comparison of the SFRs obtained from the UV-to-NIR fit with those derived from
the IR SED for the three ultra-massive galaxies can be found in the section on
Infrared luminosity and obscured star formation. Examples of best-fit SEDs for
three ultra-massive galaxies are in Extended Data Figs. 3, 4, and 5.

We note that although we could accurately remove the contamination of the
emission line in the F444W broad-band photometry, the contribution of the
nebular continuum could still be important at shorter wavelengths. To check
the possible impact of the nebular continuum, we performed SED modeling on our
sample with the original FRESCO and HST broad-band photometry and
stellar+nebular model of BAGPIPES. Here we model the nebular emission with a
range of ionization parameters log$U$ from -4 to -2 and set the other
parameters to be the same as the line-free modeling described above. With
these setups, we obtain stellar mass consistent with that from the line-free
modeling within 1$\sigma$ uncertainty. This indicates that the contribution of
the nebular continuum to the broad-band photometry of our sample is negligible
and the massive nature of our sample is robust.

To test the robustness of our results, we also use different SFH models (i.e.,
constant SFH, delayed SFH, and a combination of delayed SFH and a recent
burst) and different codes (i.e., BAGPIPES[Carnall2018] and
CIGALE[Boquien2019]) in the SED fitting for 36 optically dark/faint galaxies.
The $M_{\star}$ and SFR values derived from different codes are generally
consistent within the errors. However, there are global differences in SFR at
different SFHs. Overall, the delayed SFH+burst model returns the highest SFR,
followed by the constant SFH, and the delayed SFH returns the lowest SFR. This
is reasonable because the SFR values here trace the most recent star formation
(i.e., the last 10 Myr in BAGPIPES; instantaneous star formation in CIGALE).
In addition, we note that the fitting results are strongly constrained by the
range of initial parameter values, especially in the delayed SFH and the
delayed SFH+burst models (i.e., a long timescale of decrease $\tau$ will yield
a constant-like SFH and a short $\tau$ will lead to a burst-like SFH).
Considering the optically dark nature of our sources, so that the photometric
data points at UV-to-NIR wavelength are limited, we adopt the constant SFH –
the simplest and most constrained model – in this work. In the Extended Data
Table 2, we show the consistent $M_{\star}$ values derived from BAGPIPES and
CIGALE for the three ultra-massive galaxies studied in this paper under
different SFH models.

Overall, with the JWST FRESCO imaging and grism spectroscopic survey, we
present a robust stellar mass and redshift determination for a sample of
optically dark/faint galaxies at $z_{\rm spec}=5-9$ (Fig. 1). These galaxies
are red (F182M $-$ F444W $>1.5$), with a wide distribution of
log($M_{\rm\star}$/$M_{\odot}$) $=8.3-11.4$ and at $z_{\rm spec}>5$ (see
Extended Data Table 1). Most intriguingly, our sample contains three
spectroscopically confirmed ultra-massive galaxies
(log$M_{\rm\star}$/$M_{\odot}$ $>11.0$) at $z_{\rm spec}\sim 5.5$ and two
ultra-massive galaxies (log$M_{\rm\star}$/$M_{\odot}$ $>10.4$) at $z_{\rm
spec}\sim 7.5$. They all present very efficient star formation with
$\epsilon>0.2$ (Fig. 2). To conclude, our findings of the existence of ultra-
massive galaxies over a wide period (not only limited to the first few hundred
million years), suggest that the Universe is forming galaxies much more
efficiently than we expected.

Reliability of the $z_{\rm spec}$ for S1. Among our three ultra-massive
galaxies at $z=5-6$, S1 is the only source with a single emission line
(14$\sigma$). Although the solution of the H$\alpha$ line with $z_{\rm
H\alpha}=5.58$ is in excellent agreement with the $z_{\rm
phot}=5.61_{-0.68}^{+1.44}$ from EAZY, we here also investigate the
possibility of other emission lines. If the line with 14$\sigma$ detection was
[OIII]$\lambda 5008$ ($z_{\rm[OIII]}=7.62$) or [SIII]$\lambda 9533$
($z_{\rm[SIII]}=3.53$), we should also detect [OIII]$\lambda 4960$ or
[SIII]$\lambda 9071$ over the wavelength coverage of F444W at about 5$\sigma$,
considering the theoretical flux ratios of these two pairs ([OIII]$\lambda
5008$/[OIII]$\lambda 4960=2.92$ and [SIII]$\lambda 9533$/[SIII]$\lambda
9071=2.58$[Landt2015]). Over the wavelength coverage of the F444W filter (3.8
to 5.0 $\mu$m), the detection of only one line, i.e., the strongest line, can
only be the Pa$\alpha$ ($z_{\rm Pa\alpha}=1.30$) line or Pa$\beta$ ($z_{\rm
Pa\beta}=2.37$) line. Therefore, we mainly focus on the possibilities of
H$\alpha$, Pa$\alpha$, and Pa$\beta$ in the following discussion.

We perform SED fitting with Bagpipes, with the same parameter settings as
mentioned above. In the fitting process, we fix the $z_{\rm spec}$
corresponding to the H$\alpha$, Pa$\alpha$, and Pa$\beta$ lines, respectively.
The fitting results are shown in the Extended Data Fig. 6. The chi-square
$\chi^{2}$ values for the H$\alpha$, Pa$\beta$, and Pa$\alpha$ lines are 1.9,
6.4, and 29.9, respectively. This shows that the best fit is obtained for the
case of the H$\alpha$ line.

We also test the different redshift solutions using Bagpipes with more degrees
of freedom in the dust attenuation law. The CF00 dust attenuation
law[Charlot2000] combined with a wide Power-law slope value range from 0.5 to
2 returns $\chi^{2}$ values of 4.6, 12.3, and 47.5 for the H$\alpha$,
Pa$\beta$, and Pa$\alpha$ lines, respectively. The Salim dust attenuation
law[Salim2018] combined with a slope value range from -1.0 to 0.5 returns
$\chi^{2}$ values of 3.0, 6.0, and 22.6 for the H$\alpha$, Pa$\beta$, and
Pa$\alpha$ lines, respectively.

In addition, since S1 has a strong detection at ALMA 1.1mm ($0.95\pm 0.12$
mJy)[Gómez-Guijarro2022], we can calculate its SFRIR under different redshift
solutions, and then compare the SFRIR with the SFR value derived from UV-to-
NIR SED fitting (Extended Data Fig. 6). We calculate the SFRIR based on the
infrared luminosity, which is obtained from the IR template
library[Schreiber2018] normalized to the 1.1mm flux, following the method in
the literature[Xiao2023] (see section on Infrared luminosity and obscured star
formation). For the H$\alpha$, Pa$\beta$, and Pa$\alpha$ line solutions, the
derived SFRIR values are $643\pm 81\,M_{\odot}$yr-1, $565\pm
71\,M_{\odot}$yr-1, and $377\pm 48\,M_{\odot}$yr-1, respectively, and the SFRs
are $646^{+476}_{-199}\,M_{\odot}$yr-1, $30^{+27}_{-8}\,M_{\odot}$yr-1, and
$1.9^{+0.4}_{-0.2}\,M_{\odot}$yr-1, respectively. The SFRIR and SFR values for
the H$\alpha$ solution are in excellent agreement. However, for the Pa$\beta$
and Pa$\alpha$ solutions, their SFRIR are about 20 and 200 times higher than
their SFRs, respectively, resulting in a large discrepancy that cannot be
explained by the uncertainties caused by the differences in the SED fitting
models.

Furthermore, if we assume the observed line to be Pa$\beta$ (12821.6Å), it is
expected to be less affected by dust attenuation as it falls into the rest-
frame near-infrared. Therefore, we can make a direct comparison between SFRPaβ
and SFRIR. We calculate the SFRPaβ following the equation of SFRPaβ [M⊙yr-1] =
3.84 $\times$ 10-41 $\times$ $L$(Pa$\beta)$ [erg s-1][Reddy2023], where the
$L$(Pa$\beta)$ is the luminosity of Pa$\beta$. The derived SFRPaβ is $48\pm
3\,M_{\odot}$yr-1, which is about 12 times lower than the SFRIR.

Overall, all of the above tests provide poorer results to the Pa$\alpha$ and
Pa$\beta$ lines compared to the H$\alpha$ line. Therefore, we rule out these
two lines in S1. On the other hand, the solution of $z_{\rm spec}=5.58$ from
the H$\alpha$ line is consistent with (1) the $z_{\rm
phot}=5.61_{-0.68}^{+1.44}$ from EAZY; (2) the case of the excellent agreement
between the SFRIR and SFR values; (3) the evolutionary tracks of theoretical
dusty star-forming galaxy templates, suggesting the presence of extremely
dust-obscured massive galaxies at $z\gtrsim 5$ (see Extended Data Fig. 1); and
(4) the literature findings that optically dark/faint galaxies are generally
located at $z>3$[Wang2019, Zhou2020, Jin2022, Barrufet2023, Xiao2023].
Consequently, we consider the $z_{\rm spec}=5.58$ of H$\alpha$ line in S1
reliable.

Possible AGN contamination. To confirm the reliability of the ultra-massive
nature of the three galaxies S1, S2, and S3, we further examined whether their
line-subtracted broadband photometry could potentially be contaminated by AGN.
The contamination of AGN emission could lead to an overestimation of stellar
mass values in SED modeling that only accounts for stellar emission in models.
Type II (narrow emission lines only) AGN can produce strong, narrow line
emission, whose flux contribution has been securely subtracted from the
broadband photometry before we perform the Bagpipes SED fitting, following the
analysis procedure of typical star-forming galaxies. However, type I (broad
emission line) AGNs also show continuum emission of the accretion disk, which
cannot be removed with our method. Therefore, the most important thing is to
check if there is a type I AGN in our three ultra-massive galaxies.

In our parent sample of 36 optically dark/faint galaxies, seven galaxies have
distinct broad H$\alpha$ emission lines (full-width half maximum $v_{\rm
FWHM,H\alpha,broad}>1000$ km s-1) with a red point source morphology, in
agreement with the recently discovered “little red dot (LRD)”
population[Matthee2023, Labbé2023]. These seven galaxies are presented in
ref[Matthee2023] and marked with “LRD” in the Extended Data Table 1. Among our
three ultra-massive galaxies, none of them are in the sample of these seven
galaxies, suggesting the lack of the typical telltale signs of type I AGN in
their spectra (see Fig. 1). Specifically, we note that S3 shows a relatively
broader apparent line width compared to S1 and S2. However, the intrinsic
$v_{\rm FWHM,H\alpha}$ of S3 is only found to be about 75 km s-1. Its broad
appearance in the 1D spectra is mainly due to the coincidence of alignment
between the dispersion direction of the F444W grism and the elongated disk
morphology of S3. Furthermore, if S3 has a strong type I AGN, we should see
the H$\alpha$ line to be significantly wider than the [NII] line, but is not
actually the case. On the other hand, we also investigate the structural
properties of our three sources, to check whether they have a red point source
morphology as type I AGN. The low values of best fit Sérsic index and extended
morphology from Galfit[Peng2002] suggest that they are not dominated by a
point source (see the following section on Stellar structures of S1, S2, and
S3 for detailed discussion). All of the evidence above suggests that the
stellar masses of these galaxies are not significantly affected by an AGN and
that their ultra-massive nature is reliable.

Morphology of S1, S2, and S3 We place constraints on the structural properties
of our three sources in F444W and H$\alpha$ emission with Galfit[Peng2002].
Images are fit with a single Sersic profile convolved with the PSF allowing
the centroid, total brightness, half-light radius, Sersic index, axis ratio,
and position angle to vary. PSFs are derived from the WebbPSF
software[Perrin2014] and rotated to match the position angle of the
observations. Sources more than 3” away or more than 2.5 mag fainter are
masked; closer and brighter sources are fit simultaneously. All three sources
have a best fit Sérsic index $n<2.5$ and radii $R_{\rm e}>3$ pixels,
suggesting that they are not dominated by a point source (see Extended Data
Table 3 and Extended Data Fig. 7). In order to quantify the potential
contamination of their flux measurements by a Type I AGN, we compute the
fraction of the total flux that remains in the central resolution element of
the residual image. For all three sources, this is less than 2%, suggesting
that flux contamination from a potential AGN is not likely to be the driving
factor behind the large fiducial mass estimates.

Infrared luminosity and obscured star formation. S1, S2, and S3 are optically
dark galaxies that have star formation heavily obscured by dust, with $A_{\rm
V}>3$ mag. Although we have obtained their SFRs from UV-to-NIR SED fits
corrected for dust attenuation (Extended Data Table 1), we independently
calculated their obscured star formation here using infrared wavelength data
(SFRIR). Considering that the three sources have different far-infrared
datasets, we used different methods to obtain their total infrared luminosity
($L_{\rm IR}$; 8-1000$\mu$m rest-frame) and SFRIR.

For S1, we only have the 850$\mu$m flux from the JCMT/SCUBA-2 and ALMA band-7
observations[Cowie2018], and the 1.1mm flux from the GOODS-ALMA survey[Gómez-
Guijarro2022, Xiao2023]. For S3, we obtain 450$\mu$m and 850$\mu$m fluxes from
JCMT/SCUBA-2[Cowie2017, Barger2022], 1.2mm flux from NOEMA (PI: L. Bing,
private communication), and 1.2mm and 2mm fluxes from the NIKA2
survey[Bing2023]. These two sources are either not detected in Herschel or are
strongly affected by neighboring bright sources, so we do not include Herschel
values in our analysis here.

We first perform the FIR SED fitting to S3 with CIGALE using fixed $z_{\rm
spec}$. We use the Draine dust emission templates[Draine2014] to derive
$L_{\rm IR}$, with the same parameter settings as in ref.[Xiao2023]. The SFRIR
is then calculated based on the $L_{\rm IR}$, following the method of
ref.[Kennicutt2012]. The derived values for S3 are $L_{\rm IR}=(7.6\pm
0.38)\times 10^{12}\,L_{\odot}$ and SFR${}_{\rm IR}=1137\pm
57\,M_{\odot}$yr-1. In addition, for S1, due to its limited dataset, we use
two different approaches for the SED fitting analysis. The first method
assumes that it has the same SED shape as S3. We perform the SED fitting using
CIGALE as described above, but with fixed parameter settings obtained from the
best-fit model of S3. The second method is to re-normalize the IR
templates[Schreiber2018] to the ALMA 1.1mm flux, according to the method of
ref.[Xiao2023]. The two methods yielded very consistent values, i.e., $L_{\rm
IR}=(4.2\pm 0.21)\times 10^{12}\,L_{\odot}$ and SFR${}_{\rm IR}=625\pm
31\,M_{\odot}$yr-1 from the first methods and $L_{\rm IR}=(4.3\pm 0.55)\times
10^{12}\,L_{\odot}$ and SFR${}_{\rm IR}=643\pm 81\,M_{\odot}$yr-1 from the
second methods. In the main body of the paper, we show the values obtained by
the first method.

For S2, the so-called GN10, we obtain its $L_{\rm
IR}=1.03_{-0.15}^{+0.19}\times 10^{13}\,L_{\odot}$ and SFR${}_{\rm
IR}=1030_{-150}^{+190}\,M_{\odot}$yr-1 from ref.[Riechers2020]. These values
are also derived from CIGALE SED fits, based on extensive data from Herschel,
JCMT/SCUBA, SMA, and NOEMA observations.

The SFRIR of S1 is similar to the SFR value derived from the UV-to-NIR SED fit
(Extended Data Table 1). However, the SFRIR of S2 and S3 are about three times
higher than their SFRSED, implying that there may be hidden dust regions in
these galaxies that absorb all the UV photons, which cannot be reproduced with
a dust extinction correction[Xiao2023, Elbaz2018, Puglisi2017], or that our
assumed star-formation histories are not appropriate. Irrespective of this,
the high SFRIR values indicate that our galaxies are in a very efficient mass
assembly and accumulation process, which explains why they are so massive.

The maximal halo mass ($M_{\mathrm{halo}}^{\mathrm{max}}$). In this paper, we
define the $M_{\mathrm{halo}}^{\mathrm{max}}$ as the mass above which the
cumulative dark matter halo mass function predicts one halo to be detected in
the FRESCO survey volume ($V$). Under the most updated Planck
cosmology[Planck2020], we consider the full FRESCO survey volume as the
comoving cube enclosed by the 124 arcmin2 survey area between $z=4.9$ and
$z=9.0$, where either H$\alpha$+[NII] or H$\beta$+[OIII] fall into the
spectral coverage. This corresponds to a total survey volume of $\sim
1.2\times 10^{6}$ comoving Mpc3. As for the halo mass function, we compute it
following [Tinker2008] using the Python package hmf[Murray2013] under the same
Planck cosmology[Planck2020]. With the halo mass function, we further derive
the cumulative number density of dark matter halos above a certain mass ($n$),
as a function of redshift and halo mass. According to the definition, the
maximum halo mass at a given redshift (z) thus satisfies
n(M${}_{\mathrm{halo}}^{\mathrm{max}}$, $z$) $\times$ $V$ = 1, from which we
derive its redshift evolution and further constrain the maximal stellar mass
of galaxies expected by very simplistic galaxy assembly models in Fig. 2.

## References

This work is based on observations made with the NASA/ESA/CSA James Webb Space
Telescope. The data were obtained from the Mikulski Archive for Space
Telescopes at the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., under NASA
contract NAS 5-03127 for JWST. These observations are associated with program
# 1895.

Support for this work was provided by NASA through grant JWST-GO-01895 awarded
by the Space Telescope Science Institute, which is operated by the Association
of Universities for Research in Astronomy, Inc., under NASA contract NAS
5-26555.

This work has received funding from the Swiss State Secretariat for Education,
Research and Innovation (SERI) under contract number MB22.00072, as well as
from the Swiss National Science Foundation (SNSF) through project grant
200020_207349. The Cosmic Dawn Center (DAWN) is funded by the Danish National
Research Foundation under grant No. 140.

RPN acknowledges funding from JWST programs GO-1933 and GO-2279. Support for
this work was provided by NASA through the NASA Hubble Fellowship grant HST-
HF2-51515.001-A awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy,
Incorporated, under NASA contract NAS5-26555.

YF acknowledges support from NAOJ ALMA Scientific Research Grant number
2020-16B.

YQ acknowledges support from the Australian Research Council Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through
project number CE170100013.

MS acknowledges support from the CIDEGENT/2021/059 grant, from project
PID2019-109592GB-I00/AEI/10.13039/501100011033 from the Spanish Ministerio de
Ciencia e Innovación - Agencia Estatal de Investigación. MST also acknowledges
the financial support from the MCIN with funding from the European Union
NextGenerationEU and Generalitat Valenciana in the call Programa de Planes
Complementarios de I+D+i (PRTR 2022) Project (VAL-JPAS), reference
ASFAE/2022/025.

Cloud-based data processing and file storage for this work is provided by the
AWS Cloud Credits for Research program.

The authors declare no competing interests.

M.X. performed the majority of the present analysis and wrote the majority of
the text, with significant help from P.A.O., D.E., and L.B. E.N. performed
Galfit modeling of the sources and emission line maps. A.W. produced the
multi-wavelength catalogs used in this work. G.B. reduced the grism data. All
authors contributed to the manuscript and helped with the analysis and
interpretation.

All the raw data are publicly available through the Mikulski Archive for Space
Telescopes666https://archive.stsci.edu/ (MAST), under program ID 1895. The
FRESCO data are being released on MAST as a High-Level Science Product via
https://doi.org/10.17909/gdyc-7g80. Images are already available. The spectra
are being calibrated and will be discussed in an upcoming data paper (Brammer
et al., in prep.). For updates, please check the survey webpage: https://jwst-
fresco.astro.unige.ch/ or the MAST page
https://archive.stsci.edu/hlsp/fresco/.

Extended Data

Extended Data Table 1: Physical properties of the whole sample of optically dark/faint galaxies in FRESCO-S&N fields at $z_{\rm spec}\gtrsim 5$. (1) IDs; (2) Source ID in FRESCO; (3)(4) Right ascension and declination (J2000) of sources from JWST; (5) Spectroscopic redshift identified from H${\alpha}$ ($z_{\mathrm{{spec}}}=$ 4.86 - 6.69) or [OIII] ($z_{\mathrm{{spec}}}=$ 6.69 - 9.08) emission lines. For sources with only one emission line detected at $>8\sigma$ (flagged with a “*” exponent), their $z_{\rm spec}$ must be within the 16-84th percentile uncertainties of the $z_{\rm phot}$; (6) signal-to-noise ratio (S/N) of H${\alpha}$ or [OIII] line; (7)(8)(9) $M_{\rm\star}$, SFR, and $A_{\rm V}$ derived from UV-to-NIR SED fitting at a fixed $z_{\rm spec}$ assuming a constant SFH and a Chabrier IMF[Chabrier2003]; (10) AGN candidates, showing broad H$\alpha$ emission lines ($v_{\rm FWHM,H\alpha,broad}>1000$ km s-1) with a “little red dot (LRD)” morphology[Matthee2023]; (11) Source IDs in other work.  S.No. | ID | RA | DEC | $z_{\mathrm{{spec}}}$ |  S/N${}_{\rm H\alpha\,or\,[OIII]}$ |  log($M_{\rm\star}$) | log(SFR) | $A_{\rm V}$ | AGN | other ID
---|---|---|---|---|---|---|---|---|---|---
|  | (deg) | (deg) |  |  |  log($M_{\odot}$) |  log($M_{\odot}$yr-1) | (mag) |  |
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
| S-15496 | 03:32:24.66 | -27:44:40.08 |  5.049∗ | 17.2 |  9.72${}_{-0.12}^{+0.12}$ |  1.05${}_{-0.16}^{+0.22}$ |  1.7${}_{-0.2}^{+0.3}$ |  |
| S-11579 | 03:32:20.82 | -27:47:14.76 |  5.330∗ | 20.9 |  9.44${}_{-0.15}^{+0.11}$ |  0.82${}_{-0.14}^{+0.19}$ |  1.4${}_{-0.2}^{+0.3}$ |  |
| S-4166 | 03:32:35.44 | -27:50:31.38 |  5.402∗ | 19.6 |  9.14${}_{-0.09}^{+0.08}$ |  0.42${}_{-0.08}^{+0.11}$ |  0.2${}_{-0.1}^{+0.1}$ |  |
| S-10944 | 03:32:27.46 | -27:47:31.63 | 5.445 | 9.14 |  8.64${}_{-0.15}^{+0.15}$ |  0.41${}_{-0.15}^{+0.14}$ |  0.3${}_{-0.2}^{+0.2}$ |  |  JADES-116930[Bunker2023]
| S-11150 | 03:32:33.26 | -27:47:24.91 | 5.483 | 30.4 |  9.96${}_{-0.06}^{+0.05}$ |  1.17${}_{-0.05}^{+0.06}$ |  1.0${}_{-0.1}^{+0.1}$ |  |  JADES-204851[Bunker2023]
S1 | S-18258 | 03:32:28.91 | -27:44:31.53 |  5.579∗ | 14.3 |  11.37${}_{-0.13}^{+0.11}$ |  2.81${}_{-0.16}^{+0.24}$ |  3.2${}_{-0.2}^{+0.3}$ |  |  OFG28[Xiao2023], A2GS33[Gómez-Guijarro2022], ID68[Cowie2018], ID20[Yamaguchi2019]
| S-5661 | 03:32:20.87 | -27:49:54.85 |  5.774∗ | 10.6 |  9.56${}_{-0.12}^{+0.12}$ |  1.00${}_{-0.14}^{+0.22}$ |  1.4${}_{-0.2}^{+0.3}$ |  |
| S-12686 | 03:32:34.63 | -27:46:47.48 |  6.090∗ | 15.5 |  8.29${}_{-0.13}^{+0.15}$ |  0.33${}_{-0.14}^{+0.11}$ |  1.0${}_{-0.2}^{+0.2}$ |  |
| S-19392 | 03:32:38.73 | -27:44:15.60 | 6.815 | 20.9 |  9.97${}_{-0.07}^{+0.06}$ |  1.30${}_{-0.07}^{+0.08}$ |  1.1${}_{-0.1}^{+0.1}$ |  |  JADES-219000[Bunker2023]
| S-6209 | 03:32:32.12 | -27:49:41.72 | 7.664 | 8.3 |  10.48${}_{-0.17}^{+0.18}$ |  2.05${}_{-0.20}^{+0.22}$ |  3.7${}_{-0.4}^{+0.4}$ |  |  JADES-90354[Bunker2023]
| S-8010 | 03:32:20.99 | -27:48:53.71 | 7.846 | 26.8 |  9.30${}_{-0.14}^{+0.12}$ |  0.88${}_{-0.12}^{+0.17}$ |  0.3${}_{-0.1}^{+0.2}$ |  |  JADES-GS-53.08745-27.81492[Hainline2023]
| N-1002 | 12:36:35.50 | +62:11:04.29 |  5.050∗ | 10.9 |  9.37${}_{-0.17}^{+0.13}$ |  0.82${}_{-0.17}^{+0.27}$ |  1.4${}_{-0.3}^{+0.4}$ |  |
| N-15498 | 12:37:08.53 | +62:16:50.82 | 5.086 | 17.0 |  10.12${}_{-0.13}^{+0.11}$ |  1.37${}_{-0.14}^{+0.15}$ |  2.0${}_{-0.2}^{+0.2}$ | LRD |
| N-6079 | 12:36:55.73 | +62:13:36.33 |  5.093∗ | 15.1 |  8.98${}_{-0.11}^{+0.08}$ |  0.27${}_{-0.10}^{+0.13}$ |  0.6${}_{-0.1}^{+0.1}$ |  |
| N-14409 | 12:36:17.30 | +62:16:24.35 | 5.146 | 30.6 |  9.09${}_{-0.09}^{+0.09}$ |  0.33${}_{-0.09}^{+0.08}$ |  0.2${}_{-0.1}^{+0.1}$ | LRD |
| N-1257 | 12:36:24.68 | +62:11:17.01 | 5.167 | 21.8 |  10.04${}_{-0.19}^{+0.14}$ |  1.82${}_{-0.16}^{+0.14}$ |  1.6${}_{-0.2}^{+0.2}$ |  |
S3 | N-2663 | 12:36:56.56 | +62:12:07.37 | 5.179 | 22.6 |  11.04${}_{-0.18}^{+0.17}$ |  2.48${}_{-0.20}^{+0.29}$ |  3.4${}_{-0.4}^{+0.6}$ |  |
| N-3188 | 12:36:51.97 | +62:12:26.04 | 5.187 | 15.7 |  10.58${}_{-0.17}^{+0.14}$ |  2.30${}_{-0.17}^{+0.18}$ |  1.8${}_{-0.2}^{+0.2}$ |  |  HDF850.1[Walter2012, Herard-Demanche2023]
| N-7162 | 12:37:16.90 | +62:14:00.90 | 5.189 | 21.7 |  10.32${}_{-0.14}^{+0.11}$ |  1.71${}_{-0.12}^{+0.15}$ |  1.6${}_{-0.1}^{+0.2}$ |  |
| N-16116 | 12:36:56.62 | +62:17:07.97 |  5.194∗ | 14.6 |  9.27${}_{-0.08}^{+0.09}$ |  0.51${}_{-0.08}^{+0.09}$ |  0.5${}_{-0.1}^{+0.1}$ |  |
| N-4014 | 12:37:12.03 | +62:12:43.36 | 5.221 | 34.1 |  9.76${}_{-0.10}^{+0.09}$ |  1.05${}_{-0.11}^{+0.09}$ |  1.2${}_{-0.1}^{+0.1}$ | LRD |
| N-12839 | 12:37:22.63 | +62:15:48.11 | 5.240 | 40.2 |  9.94${}_{-0.09}^{+0.08}$ |  1.16${}_{-0.08}^{+0.08}$ |  1.0${}_{-0.1}^{+0.1}$ | LRD |
| N-13733 | 12:36:13.70 | +62:16:08.18 | 5.243 | 25.0 |  9.54${}_{-0.11}^{+0.09}$ |  0.77${}_{-0.11}^{+0.12}$ |  1.0${}_{-0.2}^{+0.2}$ | LRD |
S2 | N-7496 | 12:36:33.42 | +62:14:08.57 | 5.306 | 14.3 |  11.18${}_{-0.17}^{+0.23}$ |  2.58${}_{-0.18}^{+0.30}$ |  3.4${}_{-0.5}^{+0.8}$ |  |  GN10[Riechers2020]
| N-16813 | 12:36:43.03 | +62:17:33.12 | 5.359 | 48.5 |  9.38${}_{-0.11}^{+0.07}$ |  0.67${}_{-0.03}^{+0.05}$ |  0.1${}_{-0.0}^{+0.1}$ | LRD |
| N-9771 | 12:37:07.44 | +62:14:50.31 | 5.535 | 101.9 |  10.24${}_{-0.19}^{+0.20}$ |  2.05${}_{-0.15}^{+0.12}$ |  1.8${}_{-0.2}^{+0.2}$ | LRD |
| N-6924 | 12:37:02.72 | +62:13:55.07 |  5.535∗ | 20.3 |  9.38${}_{-0.08}^{+0.07}$ |  0.64${}_{-0.08}^{+0.08}$ |  0.3${}_{-0.1}^{+0.1}$ |  |
| N-14178 | 12:36:41.69 | +62:16:18.39 |  5.569∗ | 14.7 |  8.90${}_{-0.25}^{+0.30}$ |  0.35${}_{-0.29}^{+0.37}$ |  1.2${}_{-0.5}^{+0.9}$ |  |
| N-1459 | 12:36:41.95 | +62:11:24.47 |  6.024∗ | 15.7 |  9.46${}_{-0.09}^{+0.07}$ |  0.78${}_{-0.08}^{+0.09}$ |  0.4${}_{-0.1}^{+0.1}$ |  |
| N-12207 | 12:36:36.47 | +62:15:34.72 |  6.760∗ | 23.5 |  9.75${}_{-0.11}^{+0.09}$ |  1.11${}_{-0.09}^{+0.11}$ |  0.6${}_{-0.1}^{+0.1}$ |  |
| N-9094 | 12:36:04.62 | +62:14:36.71 | 7.038 | 73.7 |  10.48${}_{-0.12}^{+0.11}$ |  1.91${}_{-0.12}^{+0.15}$ |  1.3${}_{-0.1}^{+0.2}$ |  |
| N-489 | 12:36:47.52 | +62:10:37.27 | 7.131 | 31.4 |  9.16${}_{-0.19}^{+0.15}$ |  0.80${}_{-0.20}^{+0.24}$ |  0.7${}_{-0.3}^{+0.3}$ |  |
| N-8697 | 12:37:00.04 | +62:14:29.00 | 7.146 | 19.3 |  9.17${}_{-0.13}^{+0.12}$ |  0.69${}_{-0.15}^{+0.17}$ |  0.4${}_{-0.2}^{+0.2}$ |  |
| N-2756 | 12:36:20.04 | +62:12:09.29 | 7.190 | 18.4 |  10.00${}_{-0.10}^{+0.08}$ |  1.42${}_{-0.11}^{+0.11}$ |  0.8${}_{-0.1}^{+0.1}$ |  |
| N-11316 | 12:36:10.84 | +62:15:16.32 |  7.403∗ | 14.2 |  9.55${}_{-0.14}^{+0.11}$ |  1.09${}_{-0.16}^{+0.20}$ |  0.8${}_{-0.2}^{+0.3}$ |  |
| N-4380 | 12:36:41.09 | +62:12:53.12 | 8.613 | 16.4 |  9.54${}_{-0.15}^{+0.16}$ |  1.14${}_{-0.16}^{+0.21}$ |  1.0${}_{-0.2}^{+0.3}$ |  |  JADES-GN-189.17121+62.21476[Hainline2023]
Extended Data Table 2: Stellar masses of S1, S2, and S3. We use the Bagpipes and CIGALE SED fitting codes, assuming a constant SFH (CSFH), a delayed SFH (DSFH), and a combination of delayed SFH and a recent burst (DSFH+burst), respectively.  | Bagpipes | CIGALE
---|---|---
S.No. | CSFH | DSFH | DSFH+burst | CSFH | DSFH | DSFH+burst
S1 | 11.37${}_{-0.13}^{+0.11}$ | 11.34${}_{-0.13}^{+0.11}$ | 11.33${}_{-0.12}^{+0.11}$ | 11.16${}_{-0.22}^{+0.15}$ | 11.27${}_{-0.16}^{+0.11}$ | 11.42${}_{-0.32}^{+0.18}$
S2 | 11.18${}_{-0.17}^{+0.23}$ | 11.17${}_{-0.21}^{+0.21}$ | 11.14${}_{-0.20}^{+0.22}$ | 11.11${}_{-0.32}^{+0.18}$ | 11.40${}_{-0.10}^{+0.08}$ | 11.36${}_{-0.07}^{+0.06}$
S3 | 11.04${}_{-0.18}^{+0.17}$ | 11.01${}_{-0.15}^{+0.16}$ | 10.96${}_{-0.18}^{+0.18}$ | 10.99${}_{-0.06}^{+0.05}$ | 11.00${}_{-0.12}^{+0.09}$ | 11.04${}_{-0.16}^{+0.12}$
Extended Data Table 3: Structural parameters of S1, S2, and S3. Half light radii ($R_{\rm e}$), sérsic indices (n), and minor over major axis ratios (q) in the F444W direct image and the grism image. Though much of the structure observed in the grism image is due to the H$\alpha$ morphology, some may also be contributed by the kinematic structure and the [NII] line. The spatially extended structures and similarity of the morphologies between the grism and direct images for each object suggest that the light of these objects is not dominated by AGN.  | direct image | H$\alpha$
---|---|---
S.No. |  $R_{\rm e}$ (kpc) | n | q |  $R_{\rm e}$ (kpc) | n | q
S1 | 0.76$\pm$0.01 | 1.57$\pm$0.03 | 0.54$\pm$0.01 | 0.89$\pm$0.08 | 0.20$\pm$0.19 | 0.72$\pm$0.08
S2 | 2.57$\pm$0.25 | 1.18$\pm$0.14 | 0.64$\pm$0.04 | 1.92$\pm$0.46 | 1.55$\pm$0.43 | 0.77$\pm$0.13
S3 | 2.22$\pm$0.09 | 2.31$\pm$0.09 | 0.42$\pm$0.01 | 2.10$\pm$0.09 | 0.20$\pm$0.11 | 0.18$\pm$0.02

List of Extended Data Figures 1 Color-magnitude diagram of optically
dark/faint galaxies (circles) color-coded by $z_{\rm spec}$. F444Wc (right) is
the corrected magnitude at 4.44$\mu$m (F444W; left) after subtracting the
emission line contribution. Our parent sample of optically dark/faint galaxies
at $z_{\rm spec}\gtrsim 5$ is selected with F182M $-$ F444W $>$ 1.5 mag (red
line). Galaxies with different $M_{\star}$ show different sizes in circles.
All the optically dark/faint galaxies have $>$5$\sigma$ detections at F444W.
For those undetected at F182M, we present their 2$\sigma$ limits with arrows.
The blue-shaded region describes the distribution of all F444W-detected
sources in the FRESCO survey. For comparison, the grey triangular region shows
traditional selection criteria of $H$-dropouts[Wang2019] before the JWST
spectral era. We here assume that the HST $H$-band and Spitzer/IRAC 4.5 $\mu$m
band used in the pre-JWST selection criteria are approximate to those of JWST
F182M and F444W, respectively. We overlap the evolutionary tracks of
theoretical dusty star-forming galaxy templates (blue, green, orange, and red
regions tracing different degrees of reddening) at $z=2-6$. The numbers on the
tracks indicate the redshift. The solid and dashed tracks correspond to
log($M_{\star}/M_{\odot}$) = 11 and 10, respectively.  List of Extended Data
Figures 2 Locations of optically dark/faint galaxies compared to the SFMS in
the SFR-$M_{\star}$ plane. The extrapolated Schreiber SFMS[Schreiber2015] at
$z=5.5$, 1$\sigma$ scatter (0.5 $<$ $\Delta$MS $<$ 2, i.e., $\sim$0.3 dex),
and $\pm 3\times\Delta$MS region (0.33 $<$ $\Delta$MS $<$ 3, i.e., $\sim 0.5$
dex) are highlighted with a grey line, a dark grey shaded area, and a light
grey shaded area, respectively. $\Delta$MS $>$ 3 is commonly used to separate
MS and SB galaxies. We also plot the Popesso SFMS[Popesso2023] at $z=5.5$ with
a dashed line.  List of Extended Data Figures 3 The most extreme source S1: an
optically dark, ultra-massive galaxy at $z_{\rm spec}=5.58$. Top-left: 4′′
$\times$ 4′′ stamp images of HST/WFC3 (F160W)[Whitaker2019], ZFOURGE ($K_{\rm
s}$)[Straatman2016], Spitzer/IRAC (3.6 $\mu$m & 4.5 $\mu$m)[Stefanon2021], and
ALMA band 6 (1.1 mm)[Gómez-Guijarro2022, Xiao2023]. Before the JWST era, this
source was only detected by the (sub)millimeter observations[Cowie2018,
Yamaguchi2019, Gómez-Guijarro2022, Xiao2023]. Bottom-left: Stamps from JWST
FRESCO survey with three filters. We also show SED fitting results from
Bagpipes with fixed $z_{\rm spec}$ and the likelihood distributions of $z_{\rm
phot}$ from EAZY, which perfectly match the $z_{\rm spec}$. Right: 2D (top)
and 1D (bottom) spectra from the FRESCO NIRCam/grism (F444W filter), with the
Gaussian fit of the H${\alpha}$ emission line and the other lines around it.
List of Extended Data Figures 4 Same as Extended Data Fig. 3 but for the S2.
List of Extended Data Figures 5 Same as Extended Data Fig. 3 but for the S3.
List of Extended Data Figures 6 Different redshift solution for the S1. Since
we detect only one emission line ($\sim 14\sigma$) in S1, which could be the
H$\alpha$, Pa$\alpha$, and Pa$\beta$ lines. We show SED fitting results from
Bagpipes with the fixed $z_{\rm spec}$ corresponding to the solution of
different emission lines under the same parameter settings, respectively. The
solution of the H$\alpha$ line has the least chi-square $\chi^{2}$ (see more
discussion in Methods).  List of Extended Data Figures 7 Structural parameter
fits. We show our inference of the structural parameters of each source for
both the F444W direct image and the H$\alpha$ map from the grism spectrum. The
extended morphologies of the sources in both the direct and H$\alpha$ images
suggest that the light from these objects is not dominated by AGN.