# Probing the earliest phases in the formation of massive galaxies
with simulated HST$+$JWST imaging data from Illustris

Ángela García-Argumánez Departamento de Física de la Tierra y Astrofísica,
Facultad de CC Físicas, Universidad Complutense de Madrid, 28040, Madrid,
Spain  Instituto de Física de Partículas y del Cosmos IPARCOS, Fac. CC
Físicas, Universidad Complutense de Madrid, 28040 Madrid, Spain  Pablo G.
Pérez-González Centro de Astrobiología (CSIC-INTA), Ctra de Ajalvir km 4,
Torrejón de Ardoz, E-28850, Madrid, Spain  Armando Gil de Paz Departamento de
Física de la Tierra y Astrofísica, Facultad de CC Físicas, Universidad
Complutense de Madrid, 28040, Madrid, Spain  Instituto de Física de Partículas
y del Cosmos IPARCOS, Fac. CC Físicas, Universidad Complutense de Madrid,
28040 Madrid, Spain  Gregory F. Snyder Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD 21218, USA  Pablo Arrabal Haro NSF’s National
Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ
85719, USA  Micaela B. Bagley Department of Astronomy, The University of Texas
at Austin, Austin, TX, USA  Steven L. Finkelstein Department of Astronomy, The
University of Texas at Austin, Austin, TX, USA  Jeyhan S. Kartaltepe
Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy,
Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY
14623, USA  Anton Koekemoer Space Telescope Science Institute, 3700 San Martin
Drive, Baltimore, MD 21218, USA  Casey Papovich Department of Physics and
Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA  George
P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy,
Texas A&M University, College Station, TX, 77843-4242 USA  Nor Pirzkal Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Harry C. Ferguson Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  L. Y. Aaron Yung Astrophysics Science Division, NASA
Goddard Space Flight Center, 8800 Greenbelt Rd, Greenbelt, MD 20771, USA
Marianna Annunziatella Centro de Astrobiología (CSIC-INTA), Ctra de Ajalvir km
4, Torrejón de Ardoz, E-28850, Madrid, Spain  Nikko J. Cleri Department of
Physics and Astronomy, Texas A&M University, College Station, TX, 77843-4242
USA  George P. and Cynthia Woods Mitchell Institute for Fundamental Physics
and Astronomy, Texas A&M University, College Station, TX, 77843-4242 USA  M.
C. Cooper Department of Physics & Astronomy, University of California, Irvine,
4129 Reines Hall, Irvine, CA 92697, USA  Luca Costantin Centro de
Astrobiología (CSIC-INTA), Ctra de Ajalvir km 4, Torrejón de Ardoz, E-28850,
Madrid, Spain  Benne W. Holwerda University of Louisville, Department of
Physics and Astronomy, 102 Natural Science Building, Louisville, KY 40292, USA
Rosa María Mérida González Centro de Astrobiología (CSIC-INTA), Ctra de
Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain  Caitlin Rose
Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy,
Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY
14623, USA  Mauro Giavalisco Astronomy Department, University of
Massachusetts, Amherst, MA 01003, U.S.A.  Norman A. Grogin Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Dale D.
Kocevski Department of Physics and Astronomy, Colby College, Waterville, ME
04901, USA

###### Abstract

We use the Illustris-1 simulation to explore the capabilities of the Hubble
and James Webb Space Telescope data to analyze the stellar populations in
high-redshift galaxies, taking advantage of the combined depth, spatial
resolution, and wavelength coverage. For that purpose, we use simulated broad-
band ACS, WFC3 and NIRCam data and 2-dimensional stellar population synthesis
(2D-SPS) to derive the integrated star formation history (SFH) of massive
(M${}_{\ast}>10^{10}\,$M⊙) simulated galaxies at $1<z<4$ that evolve into a
local M${}_{\ast}>10^{11}\,$M⊙ galaxy. In particular, we explore the potential
of HST and JWST datasets reaching a depth similar to those of the CANDELS and
ongoing CEERS observations, respectively, and concentrate on determining the
capabilities of this dataset for characterizing the first episodes in the SFH
of local M${}_{\ast}>10^{11}\,$M⊙ galaxies by studying their progenitors at
$z>1$. The 2D-SPS method presented in this paper has been calibrated to
robustly recover the cosmic times when the first star formation episodes
occurred in massive galaxies, i.e., the first stages in their integrated SFHs.
In particular, we discuss the times when the first 1% to 50% of their total
stellar mass formed in the simulation. We demonstrate that we can recover
these ages with typical median systematic offset of less than 5% and scatter
around 20%-30%. According to our measurements on Illustris data, we are able
to recover that local M${}_{\ast}>10^{11}\,$M⊙ galaxies would have started
their formation by $z=16$, forming the first 5% of their stellar mass present
at $z\sim 1$ by $z=4.5$, 10% by $z=3.7$, and 25% by $z=2.7$.

Astronomy data analysis (1858) — 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)

††journal: AJ††software: synthesizer (Pérez-González et al., 2003, 2008),
Photutils (Bradley et al., 2019), SExtractor (Bertin & Arnouts, 1996), IRAF
(Tody, 1993), Astropy (Astropy Collaboration et al., 2013, 2018), Colossus
(Diemer, 2018), NumPy (van der Walt et al., 2011), Pandas (Reback et al.,
2020), SciPy (Virtanen et al., 2020), and Matplotlib (Hunter, 2007).

##  1 INTRODUCTION

The most massive galaxies in the nearby Universe are generally quiescent and
present relatively old stellar populations (Renzini 2006 and references
therein). This is consistent with the ‘downsizing’ scenario, first introduced
by Cowie et al. (1996), and which states that the most massive galaxies
started and completed their stellar mass assembly earlier than less massive
ones. Nevertheless, the mechanisms involved and, in particular, the role of
feedback in the downsizing scenario of galaxy formation remain unclear. We can
summarize this downsizing scenario and its implications in our understanding
of galaxy formation by the following key questions: (i) When did the first
massive galaxies appear in the lifetime of the Universe? (ii) When did the
first star formation stages happen in massive galaxy formation, and how fast
and bursty were they? (iii) How did stellar mass assemble in massive galaxies?
One way to address these questions is to look further and fainter in order to
search for the most likely progenitors of local massive galaxies at different
redshifts, and subsequently analyze their stellar populations, especially
their Star Formation Histories (SFHs).

In the last few years, significant work has been conducted to study the
stellar populations of the most massive galaxies at $z>3$ from the integrated
emission of photometrically-selected (and sometimes spectroscopically-
confirmed) samples (e.g. Glazebrook et al. 2017, Alcalde Pampliega et al.
2019, Forrest et al. 2020a, b, Marsan et al. 2022 and references therein).
However, there are also spectroscopic and photometric studies that show that
the physical properties of galaxies present systematic trends both radially,
mainly, but also azimuthally in the local Universe (see, e.g., Di Matteo et
al. 2013, Tacchella et al. 2015, Nelson et al. 2016, Wang et al. 2017, Ho et
al. 2018, Sánchez-Menguiano et al. 2020, Abdurro’uf et al. 2022, Chamorro-
Cazorla et al. 2022, Bellardini et al. 2022). For these kinds of studies, the
stellar population synthesis (SPS) modeling in two dimensions (2D) has become
an essential tool to derive these subgalactic-scale resolved properties by
analyzing their spatially-resolved emission. If applied to massive progenitors
at high redshift, this spatially-resolved analysis can provide further clues
on the role of different mechanisms (e.g., internal vs. external, secular vs.
fast) on the evolution of massive galaxies, since the different mechanisms
proposed act at different scales and have a different impact as a function of
galactocentric distance (Dekel & Burkert 2014, Zolotov et al. 2015, Tacchella
et al. 2015, 2016, 2018, 2019, Akhshik et al. 2022) and for different
morphological components (e.g., Méndez-Abreu et al. 2021, Johnston et al. 2022
for local galaxies, and Costantin et al. 2021, 2022 for higher-redshift
galaxies).

Still, we cannot ignore that the main drawback of the stellar population
synthesis is the intrinsic degeneracies associated with the study of the
emission from stars, even within individual pixels. Indeed, limitations in the
spectral resolution, wavelength coverage and/or signal-to-noise of the data,
result in strong degeneracies in the parameter space of physical properties
such as age, star formation timescale, metallicity, or attenuation by dust.
The correlations among these parameters are more or less difficult to
disentangle depending on the value of each specific parameter (e.g., young
ages, such as those expected in high-redshift galaxies, are less prone to some
of these degeneracies), as well as the mentioned observational errors and
wavelength coverage (see e.g. Gil de Paz & Madore 2002).

The latter problems are especially difficult to tackle as we go to higher
redshifts and/or study less massive galaxies, since galaxies are generally too
faint to reach their continuum level with high spectral resolution. This makes
the SPS analysis of progenitors of massive nearby galaxies at high redshift
extremely challenging, as we cannot use a large enough sample of them reaching
faint magnitudes to account for a significant fraction of their whole mass
and, thus, obtain statistically robust results. These large samples are still
harder to get when using spectroscopic (instead of photometric) data, which
are resource-demanding and therefore relatively scarce. In addition, due the
technological limitations of our observatories, it is hard to cover wide
spectral ranges to account for different populations and effects such as dust
extinction or other degeneracies.

In this regard, the James Webb Space Telescope (JWST; Gardner et al. 2006)
will play an important part in overcoming these issues, thanks to its
unprecedented sensitivity, high spatial resolution, wide wavelength coverage,
sensitive to stars of different ages, and variety of spectral resolutions
(from spectroscopic $R\sim 3000$ to photometric $R\sim 7$ and spectro-
photometric $R=100$). Our approach is to use spatially-resolved
multiwavelength broad-band data from JWST, combined with already existing
broad-band data from the Hubble Space Telescope (HST), to extract Spectral
Energy Distributions (SEDs) of progenitors of nearby massive galaxies at high
redshift and apply 2D-SPS modeling on them. The determination of stellar
population parameters in two dimensions, and not only for the galaxy as a
whole, will help us to analyze stellar mass distributions inside galaxies and
to recover realistic integrated galactic SFHs. The power of analyzing the
stellar populations in 2D resides in the fact that smaller regions of a given
galaxy should have its observational properties driven by a more simple SFH,
which can therefore be characterized with fewer parameters compared to that
required to characterize an entire galaxy (which easily could count with
several, quite different stellar populations). Analyzing the stellar content
of high-redshift systems in 2D will be fundamental to identify the
evolutionary stages a galaxy can undergo regarding its stellar content
evolution and to determine when its assembly began.

This work aims at determining the robustness of spatially-resolved SFHs for
massive galaxies at $1<z<4$, derived from the 2D-SPS analysis of simulated
HST$+$JWST broad-band photometry, to infer when these massive progenitors
began their stellar mass assembly. The reason for focusing on this $1<z<4$
redshift interval, for which JWST is expected to provide high-quality data, is
that it includes the so-called ‘cosmic noon’, i.e., the epoch of the Universe
where the cosmic star formation rate density history was maximum at $z\sim 2$
and where a considerable fraction of the local stellar mass was formed: half
of the present-day stellar mass was formed before $z=1.3$ (see Madau &
Dickinson 2014 and references therein). The results of this work will help us
identify the first progenitors of massive galaxies around cosmic noon, an
epoch of the Universe of great interest for many JWST studies in the upcoming
years.

In order to establish a methodology and test its performance, we use galaxy
images in different bands simulated by the Illustris Project. One of the most
appealing features of recent cosmological simulations such as Illustris-1
(Vogelsberger et al. 2014a, Nelson et al. 2015) is that some of them provide
synthetic stellar images generated for simulated galaxies at different
redshifts in common broad-band filters from HST or JWST (Torrey et al. 2015)
and with “realistic” morphologies (at least, relative to the morphologies
provided by typical semi-analytical models). These images, when compared to
the available information of the stellar particles for each galaxy in the
simulation, make them the perfect benchmark to test how successful our 2D
stellar population synthesis method is at recovering the early formation of
massive galaxies at high redshift when upcoming JWST broad-band photometry
data become available.

The present paper is structured as follows. In Section 2, we give a brief
overview of the Illustris simulation and explain how the sample of simulated
galaxies has been selected from the synthetic images. In Section 3, we present
the processing of the images and the photometric method conducted to build
SEDs in 2D. Section 4 describes the SPS modeling of the SEDs measured in 2D
and the derivation of the integrated SFHs from these 2D-SPS analysis. In
Section 5, we evaluate the success of our method in recovering the earliest
phases in the formation of each galaxy in the sample. Finally, in Sections 6
and 7, we present our results and outline our conclusions, respectively.

Throughout this paper, we adopt the same cosmology as specified in the
Illustris-1 simulation: $\Lambda$CDM cosmology with $\Omega_{m}$ = 0.2726,
$\Omega_{\Lambda}=0.7274$, $\Omega_{b}=0.0456$, and $h=0.704$. We assume a
Kroupa (2001) initial mass function (IMF). All magnitudes presented in this
work are calculated using the AB system (Oke & Gunn 1983).

##  2 DATASET

In this paper, we present a method to study the star formation history of
massive, spatially resolved high-redshift galaxies using HST and JWST
photometric data. We apply this method to a sample of galaxies simulated by
the Illustris Project. Our aim is to assess the utility of HST plus JWST
combined datasets for the analysis of the earliest evolutionary phases of
nearby massive galaxies. For that purpose, our approach consists in selecting
the progenitor galaxies at high redshift of nearby massive galaxies, and then
analyzing their pixel-by-pixel and resulting integrated-galactic SFHs. In this
Section, we describe the simulations, the sample, and the data.

###  2.1 Illustris Simulation

Figure 1: Left panel: Main Sequence plot for those galaxies at $1<z<4$
included in the mock survey presented in S+17 which are progenitors of very
massive galaxies (M${}_{\ast}>10^{11}$ M⊙) at $z=0$. All progenitors are
color-coded by redshift. In this paper, we concentrate on the analysis of the
most massive progenitors (M${}_{\ast}>10^{10}$ M⊙), which are plotted with
star symbols. The Main Sequence found for all the Illustris simulated galaxies
at different redshifts (Sparre et al. 2015) has also been plotted. Middle
panel: postage stamp images for some representative examples of our galaxies
(size $2.5\arcsec\times 2.5\arcsec$). These RGB images are created using
ACS/F1814W, NIRCam/F200W, and NIRCam/F277W asinh-scaled images as B, G and R
filters, respectively. The position of these galaxies in the left panel has
been highlighted with a white dot inside. Right panel: $g-r$ color vs. stellar
mass diagram for descendants at $z=0$ of all galaxies at $1<z<4$
(independently of their stellar mass) in the S+17 mock survey images.
Descendants of the final sample of 221 galaxies analyzed in this paper are
plotted as stars.

The Illustris Project is a set of hydrodynamical simulations of a (106.5 Mpc)3
periodic cosmological volume that trace the evolution of dark matter, gas,
stars, and supermassive black holes from $z=127$ to $z=0$ (Vogelsberger et al.
2014a, Nelson et al. 2015). Illustris comprises six different runs:
Illustris-(1,2,3) and Illustris-(1,2,3)-Dark, where the former include a
baryonic physical component in addition to the dark-matter content, and the
latter the dark-matter component alone. We make use of the Illustris-1 run,
the one with the highest-resolution box of these six runs in terms of the
number of resolution elements and their masses: it follows the evolution of
18203 dark-matter particles with a mass of 6.26 $\times$ 106 M⊙ each, and
18203 baryonic particles with an initial mass of 1.26 $\times$ $10^{6}$ M⊙.

This work is based on the simulated galaxies and broad-band images from the
Illustris-1 supplementary data catalogs published by Snyder et al. (2017,
hereafter S+17) and available via
MAST111https://archive.stsci.edu/prepds/illustris/. These catalogs correspond
to three synthetic deep survey square images in three different fields, which
are labeled as field A, B, and C, each $2.8\arcmin\times 2.8\arcmin$ in size.
These three images, also called “mock ultra deep fields”, are provided down to
the spatial resolution of HST and JWST, among other observatories, in a wide
range of broad-band filters, imitating the conditions of real galaxy surveys.
Publicly available catalogs associated to these images include galaxies in the
survey images whose rest-frame g-band apparent magnitude is g $<30.0$ mag
(19,347 galaxies).

###  2.2 Sample Selection

Starting from all the galaxies included in the S+17 mock images, we select
those at $1<z<4$ which will evolve into a very massive (M${}_{\ast}>10^{11}$
M⊙) descendant at $z=0$, as tracked forward in time via the Illustris merger
trees (Rodriguez-Gomez et al., 2015). This means that we select $1<z<4$
progenitors in the mock images of $z=0$ very massive galaxies in the whole
Illustris-1 simulation. Among all progenitor galaxies at $1<z<4$ in the
images, we restrict our analysis to massive galaxy progenitors with
M${}_{\ast}>10^{10}$ M⊙, for which very high quality JWST data will be
available in the near future. This implies that, for the same 1011 M⊙ galaxy
at $z=0$, all its massive progenitors at $1<z<4$ in the images are considered
in our sample (e.g. two massive galaxies at $1<z<4$ could merge to form a
M${}_{\ast}>10^{11}$ M⊙ galaxy at $z=0$).

The left panel of Fig. 1 shows the location in the star formation rate (SFR)
vs. stellar mass plane of all the progenitor galaxies at $1<z<4$ in the S+17
simulated images that have a M${}_{\ast}>10^{11}$ M⊙ at $z=0$. These 2,994
progenitor galaxies shown in the left panel of Fig. 1 represent 19% of all
galaxies located at $1<z<4$ in the images. This is equivalent to saying that
81% of galaxies at $1<z<4$ in the images do not end up as a very massive
galaxy at $z=0$. Out of that fraction of $1<z<4$ galaxies with a very massive
descendant at $z=0$, our selected sample is composed of the 221 progenitor
galaxies with M${}_{\ast}>10^{10}$ M⊙ (represented by stars in Fig. 1, left
panel). Among them, two of them appear simultaneously in two of the three mock
ultra-deep fields with different orientations, and will be considered as
independent galaxies in this work regarding the photometric analysis and the
derivation of the SFHs. In fact, the initial selected sample from the images
is composed of 248 galaxies, but only 221 (+ 2 repetitions) are kept until the
final analysis (the rest are discarded due to problems when attempting to
download the data from the simulation database or when measuring photometry
for them, as explained in Section 3.2). These 248 galaxies in the initial
sample of massive progenitors represent the 64% of all M${}_{\ast}>10^{10}$ M⊙
galaxies at $1<z<4$ in the images (388 galaxies), i.e., if we were to select
all massive galaxies at $1<z<4$ in the images, only 64% of them would actually
become a M${}_{\ast}>10^{11}$ M⊙ galaxy at $z=0$. Interestingly, this
percentage rises to 68% and 85% when the mass cutoff is set to $10^{10.1}$ and
$10^{10.5}$ M⊙, respectively. In Fig. 1, we also show the galaxy Main Sequence
for Illustris galaxies, as determined by Sparre et al. (2015).

The right panel of Fig. 1 shows the $g-r$ color-stellar mass diagram for
descendants at $z=0$ of all galaxies at $1<z<4$ in the mock images of S+17.
Only descendants with log(M∗/M⊙)$\gtrsim 8.68$ are shown, for which integrated
photometry data is available (5,498 galaxies). Since descendants at $z=0$ of
galaxies in S+17 mock images have been traced via the Illustris merger tree
and do not appear in these survey images, the integrated photometry magnitudes
used to build this diagram have been extracted from the Torrey et al. (2015)
synthetic individual galaxy images. But these mock galaxy images are only
available for $z=0$ descendants with a minimum stellar mass of
log(M∗/M⊙)$\gtrsim 8.68$ ($\sim 63\%$ of all the 8,768 $z=0$ descendants of
$1<z<4$ galaxies in the images).

Our final 221 galaxies selected from these images evolve to 132 (unique) very
massive galaxies at $z=0$, which means that some galaxies in our final sample
at $1<z<4$ have the same descendant at $z=0$. These descendants are called
$z=0$ descendants of our main sample, hereafter, and are shown as violet stars
in the right panel of Fig. 1. We note that 51% of all $z=0$
M${}_{\ast}>10^{11}$ M⊙ with a progenitor galaxy at $1<z<4$ in the images do
not come from any of the galaxies in our initial sample of 248 $1<z<4$,
M${}_{\ast}>10^{10}$ M⊙ systems. The main progenitor at $z>1$ in the images of
those nearby massive galaxies presents a typical stellar mass
$\log($M∗/M${}_{\odot})=9.12^{9.65}_{8.48}$ and a typical redshift of
$z=2.71^{3.42}_{2.01}$ (median values and first and third quartiles).

The 221 galaxies in the sample account for 28% of the total stellar mass
present in their 132 $z=0$ descendants. These descendants are predominantly
Red Sequence galaxies (RS; approximately 73%), but there are also a few
galaxies in the Green Valley (GV; 14%) and in the Blue Cloud (BC; 13%). In the
case of the M${}_{\ast}>10^{11}$ M⊙ descendants of progenitor galaxies at
$1<z<4$ which are not massive, their location in the diagram is similar, but
with a higher fraction of them in the BC: 63% in the RS, 16% in the GV, and
21% in the BC.

Figure 2: Histograms for our $1<z<4$ sample of 221 simulated galaxies:
redshift (top left panel), total stellar mass (top right), SFR (bottom left),
and stellar half-mass radius (bottom right). These properties have been
extracted from the Illustris-1 database. Median and quartiles are shown as
segments on the top.

The galaxies in our sample are projected into three $2.8\times 2.8$ arcmin2
area (S+17). The typical stellar mass, star-formation rate (SFR), redshift and
stellar half-mass radius (i.e., the radius enclosing half of the total stellar
mass of the galaxy; r${}_{\text{hm}}$) of galaxies in the sample are
$\langle\log($M${}_{\ast})\rangle=10.4^{10.7}_{10.1}$ M⊙,
$\langle$SFR$\rangle=25^{43}_{12}$ M⊙/yr, $\langle z\rangle=2.0^{2.5}_{1.6}$,
and $\langle$r${}_{\text{hm}}\rangle=4.7^{5.8}_{3.6}$ kpc (median and
quartiles), respectively. Fig. 2 shows the histograms of these properties for
our sample of 221 galaxies.

##  3 PHOTOMETRIC DATA FROM THE ILLUSTRIS SIMULATION

In this section, we describe the simulated deep survey images used in this
work and how these have been processed to obtain SEDs for all galaxies in our
sample. We measure integrated-light photometry for each galaxy, and we also
consider spatially resolved SEDs, all constructed with data in several HST and
JWST filters. Apart from simulated images, we also use redshifts and physical
properties extracted from the Illustris database. Throughout this paper, we
refer to $`$ground-truth’ or $`$reference’ properties as those galaxy
properties derived either from the catalogs associated to the S+17 mock images
or from the Illustris database, in contrast to the properties obtained from
SED fits to stellar population synthesis models for different regions in each
galaxy, what we call 2D-SPS-derived properties.

###  3.1 Original Galaxy Images

As mentioned in subsection 2.1, this work uses the Illustris-1 “mock ultra-
deep fields” described in S+17, which consist of three synthetic deep survey
square images, each 2.8 arcmin across. Briefly, these images have been created
using the light cone technique described in Kitzbichler & White (2007) through
the periodic Illustris-1-simulation volume until z $\sim$ 18\. The output of
the simulation at each redshift is processed with the spectral synthesis code
SUNRISE (Jonsson 2006; Jonsson et al. 2010) to create synthetic images in
arbitrary filters. This is done by assigning SEDs to the stellar particles
according to their mass, age, and metallicity, and by projecting these
quantities from the simulation space to pre-defined hypothetical cameras. For
the creation of these images, Starburst99 stellar population models (Leitherer
et al. 1999) were used, with a Kroupa (2001) initial mass function (IMF). The
rest of SUNRISE configuration parameters are stored in the image headers.

These initial deep survey images are available to be used in 34 broad-band
filters onboard the HST/ACS and WFC3, JWST/NIRCam and MIRI, and Roman/Wfi. For
this work, we consider 15 HST and JWST broad-band filters in the optical and
near-infrared (see Table 1), which will be available for several cosmological
fields covered by HST legacy projects such as the Cosmic Assembly Near-
infrared Deep Extragalactic Legacy Survey fields (CANDELS; Grogin et al. 2011;
Koekemoer et al. 2011), and JWST Guaranteed Time Observations (GTO), Cycle 1
Guest Observers (GO1) or Early Release Science (ERS) programs such as the MIRI
Deep Survey (Norgaard-Nielsen & Perez-Gonzalez 2017), the JWST Advanced Deep
Extragalactic Survey (JADES; Williams et al. 2018, Rieke et al. 2019), the
Public Release IMaging for Extragalactic Research (PRIMER; Dunlop et al. 2021)
or the Cosmic Evolution Early Release Science (CEERS; finkelstein et al.
2017). The filterset covers the optical through mid-infrared observed spectral
region, corresponding to rest-frame UV to near-infrared wavelengths at
$1<z<4$, an adequate range to study the emission from young and old stellar
populations (i.e., providing information about the SFH). The survey images
released by S+17 are noise-free synthetic images with different spatial
resolution depending on the instrument and telescope used (see Table 1) and
available with or without considering the point spread function (PSF)
degradation for each filter. For this work, we use as starting images those
without the PSF model applied, with pixel scales provided in Table 1.

Table 1: Photometric broad-band images used in this work

|  $\mathrm{\lambda}\mathrm{{}_{central}}$ | Width | Pixel scaleaaPixel size in original images. Final (matched) pixel size is 0.06$\arcsec$. | 5$\mathrm{\sigma}$ depthbbHST depths from the CANDELS/3D-HST catalogs (Skelton et al. 2014): median $5\sigma$ depth calculated from the errors of objects in the final catalogs (apertures of 0.7″) for all 5 CANDELS fields. JWST depths corresponding to the CEERS proposal: values represent the planned 5$\sigma$ point source depths per filter, assuming a total integration time of $2867\,$s for all NIRCam filters except for F115W ($5734\,$s).
---|---|---|---|---
Instrument .Filter |  |  |  |  |  |
| ($\mathrm{\mu}$m) | ($\mathrm{\mu}$m) | (”/pix) | (mag)
ACS | F435W | 0.4318 | 0.0993 | 0.03 | 27.3 |  |
@HST | F606W | 0.5919 | 0.2225 | 0.03 | 27.4 |  |
| F775W | 0.7693 | 0.1491 | 0.03 | 26.9 |  |
| F814W | 0.8057 | 0.2358 | 0.03 | 27.2 |  |
| F850LP | 0.9036 | 0.2092 | 0.03 | 26.5 |  |
sWFC3 | F105W | 1.0585 | 0.2653 | 0.06 | 26.1ccThis limiting magnitude is not included in the CANDELS/3D-HST catalogs from Skelton et al. 2014. We assume the same value as F125W, since the magnitudes for this band are slightly fainter/brighter than those of F105W in the CANDELS survey. |  |
@HST | F125W | 1.2471 | 0.2867 | 0.06 | 26.1 |  |
| F140W | 1.3924 | 0.3760 | 0.06 | 25.6 |  |
| F160W | 1.5396 | 0.2744 | 0.06 | 26.4 |  |
NIRCam | F115W | 1.1512 | 0.2426 | 0.032 | 29.2 |  |
@JWST | F150W | 1.5017 | 0.3309 | 0.032 | 28.9 |  |
| F200W | 1.9905 | 0.4654 | 0.032 | 29.0 |  |
| F277W | 2.7861 | 0.7117 | 0.065 | 29.2 |  |
| F356W | 3.5594 | 0.8163 | 0.065 | 29.0 |  |
| F444W | 4.4457 | 1.1197 | 0.065 | 28.6 |  |

Notes.

###  3.2 Photometric Measurements

We process the original idealized (without PSF and noise) images for HST and
JWST filters in order to imitate CANDELS data (Grogin et al. 2011; Koekemoer
et al. 2011) and the recently-begun CEERS observations (finkelstein et al.
2017), respectively. The reason for choosing the CEERS survey as a test case
is that, although depths for subsequent JWST surveys are expected to be
better, CEERS is expected to provide one of the first publicly available
datasets for deep-field JWST observations. Cycle 1 GO / Treasury programs such
as PRIMER are expected to achieve similar depths.

First, we match the pixel size of the images to that of HST/WFC3 (0.06
arcsec). Second, we convolve the image for each filter with the PSF model for
the WFC3/F160W filter (FWHM $\thickapprox$ 0.19 arcsec). Finally, we add
Gaussian sky noise to the HST and JWST images to achieve the same depth as in
the CANDELS and CEERS observations, respectively. We degrade HST images adding
noise in order to match the median background rms noise measured around
CANDELS massive galaxies (M${}_{\ast}>10^{10}$ M⊙) at $1<z<4$. For this
calculation, we use actual regions covered by CANDELS/Wide, CANDELS/Deep
(Grogin et al. 2011, Koekemoer et al. 2011), and the Hubble Ultra Deep field
(HUDF; Beckwith et al. 2006, Oesch et al. 2010a, Ellis et al. 2013, Koekemoer
et al. 2013, Illingworth et al. 2013). For JWST images, the sky noise is
estimated using the official JWST exposure time calculator
Pandeia222https://jwst.etc.stsci.edu/. Pandeia Version: 1.7 (Pontoppidan et
al. 2016), by measuring the predicted signal-to-noise ratio (SNR) as a
function of the pixels surface-brightness in CEERS observations for the
selected JWST broad-band filters. The initial spatial resolution of the images
for each filter (before registering) and final depths (after the noise
addition) are given in Table 1. The limiting magnitude has been estimated
(without aperture correction) from the sky rms value adopting a circular
aperture with a fixed radius of 0.2$\arcsec$ at a 5$\mathrm{\sigma}$ level.

We use these processed (registered, PSF-matched, and sky noise-added) deep
survey images to measure photometry. To do that, we first generate a
segmentation map by running SExtractor (Bertin & Arnouts 1996) on the
WFC3/F160W image. Typical parameters were used. The results did not vary much
when changing the relevant parameters given that we are dealing with
relatively bright galaxies. In general, there is good agreement between
sources detected in the segmentation map and galaxy positions given by the
Illustris catalogs based on the original images. Then, this segmentation map
is combined with the information of the sample in the Illustris catalogs to
create circular galaxy apertures that enclose the integrated emission from all
galaxies in the sample (see 3.2.1 below). We choose to use galaxy positions
given by the Illustris catalogs instead of the ones provided by SExtractor in
order to make the future comparison between our results derived from
2D-photometry and from the simulated particles in the galaxy more fair
(otherwise, the centers for the galaxy photometric apertures would be shifted
from the galaxy centers given by the simulated particles).

We build galaxy SEDs by measuring multi-wavelength photometry in two different
ways: inside these circular galaxy apertures, which enclose the entire flux of
the galaxy (hereinafter “integrated photometry”), and for small parts of the
galaxy defined after creating a grid inside this circular aperture with cell
size equal to 3$\times$3 pixels, roughly the area of the FWHM of the PSF-
homogenized dataset (“2D photometry”, hereafter). In both cases, the flux
uncertainty is estimated as $\sqrt{N}\times\sigma$, where $N$ is the number of
pixels within the aperture and $\sigma$ the rms of the sky for that filter.

####  3.2.1 Integrated Photometry

In order to define the integrated-photometry apertures, we use both the
segmentation map and information extracted from the Illustris database, in
particular, the galaxy centers and stellar half-mass radii. The apertures are
first defined in the segmentation map around galaxy centers with an initial
radius equal to twice the stellar half-mass radius ($r_{\text{hm}}$). We
decided to use the $r_{\text{hm}}$ values from the Illustris database because
many of the galaxy properties provided by these simulations refer to this
radius (and to $2\times r_{\text{hm}}$) and, typically, papers compare the
results based on observations to what Illustris provides for $2\times
r_{\text{hm}}$ (see, e.g., Vogelsberger et al. 2014b, Sales et al. 2015, Cook
et al. 2016, Elias et al. 2018, Valentino et al. 2020) . With the following
procedure, we test how those apertures compare to what is directly measured in
the simulated images. The initial radius is first reduced to minimize the
contamination from other sources until 80% of the pixels within the aperture
belong to the considered galaxy, i.e., we allow up to 20% of pixels from
neighbors. If the number of sky pixels (defined as those not belonging to any
source) within this new aperture is greater than 10%, we further reduce the
radius to decrease that number below 10%. We refer to the final photometric
aperture radius as $r_{\text{phot}}$ hereafter. We discard galaxies (from the
initial selected sample of 248 galaxies) whose final aperture radius is
$r<0.75\times r_{\text{hm}}$. This showed to be effective in removing sources
overdeblended by SExtractor and galaxies whose initial aperture (with
$r=2\times r_{\text{hm}}$) is partially beyond the edges of the simulated
images. We also discard galaxies which present surface brightness values
fainter than 25 mag/arcsec2 in WFC3/F160W in their brightest pixel within the
final aperture.

Typically, the fraction of pixels within the final integrated photometric
apertures belonging to other nearby galaxies is very low, less than 1% for 75%
of the sample of 221 galaxies. The maximum percentage of pixels from neighbor
galaxies within the integrated apertures is $\sim$20%, but only 27 out of 221
galaxies have more than 10% of pixels belonging to other sources.

In Fig. 3 (top panel) we show the histogram for the final radii of the
apertures as a fraction of $2\times r_{\text{hm}}$ for the 221 galaxies in the
final sample. This means we do not have apertures in this histogram at
$r/(2\times r_{\text{hm}})<0.375$, due to the threshold imposed on the final
apertures to be kept: $r>0.75\times r_{\text{hm}}$. We find the median and
quartiles for $r_{\text{phot}}/(2\times r_{\text{hm}})$ are
$0.83^{1.00}_{0.66}$. It can be noticed that a considerable number of galaxies
presents radii very similar or equal to the (starting) aperture radius of
$2\times r_{\text{hm}}$: $\sim 40$% and $\sim 37$% of galaxies have
$r_{\text{phot}}/(2\times r_{\text{hm}})\geq 0.95$ and 0.999, respectively. On
the contrary, only $\sim 8$% of galaxies have aperture radii smaller than
$r_{\text{hm}}$. But assuming $2\times r_{\text{hm}}$ as the best aperture to
compare our results with simulations has demonstrated not to work for a
significant number of galaxies.

Figure 3: Top panel: histogram for the radii of the photometric apertures of
the 221 galaxies in the final sample as a fraction of $2\times r_{\text{hm}}$
radius, the typical radius used by Illustris to compare with observations. The
cumulative fraction of galaxies in shown as a red line. Bottom panel:
histogram for the fraction of stellar mass enclosed by the photometric
apertures (blue filled histogram) and for a radius of $2\times r_{\text{hm}}$
(orange hatched) with respect to the total stellar mass in the galaxy. In both
cases, we neglect neighboring galaxies. These total masses have been extracted
from the simulated particles belonging to each galaxy in the Illustris
database. At the top of both panels, we show the median and quartiles for each
histogram.

The stellar masses enclosed by our final apertures, a fraction of the total
stellar mass of each galaxy, can be seen in the bottom panel of Fig. 3. As a
comparison, we also show the histogram for the fraction of stellar mass
enclosed by a radius of $2\times r_{\text{hm}}$. For both radii, the masses
have been extracted from the Illustris database, i.e., they are ground-truth
masses calculated by adding up the simulated stellar particles belonging to
each galaxy which are closer to the galaxy center than the radii considered.
The typical percentages of the stellar masses enclosed by these two apertures
are (median and quartiles) $65_{57}^{72}$% for the final photometric apertures
and $75_{71}^{80}$% for a radius of $2\times r_{\text{hm}}$. These percentages
refer to total stellar masses provided by the Illustris database or using all
particles (in both cases, referring to the whole dark matter halo, i.e., they
include regions with very low mass surface densities, whose emission is well
below our observational limits). This procedure aims at reproducing that to be
applied to actual HST+JWST observations, so we should be able to evaluate its
impact on reproducing ground-truth properties.

To measure the integrated flux for each galaxy, we only consider pixels whose
center lies within the final circular aperture and which do not belong to
other sources (according to the segmentation map). The fluxes of these pixels
within the aperture (considered galaxy and sky pixels) are added together to
build the integrated SED. In Table 2 we show the typical integrated magnitudes
(median and quartiles) of our 221 galaxies in several filters and in different
redshift bins.

Table 2: Integrated magnitudes of our $1<z<4$ sample

|  |  $m_{\mathrm{AB}}$ (mag)
---|---|---
Instrument | Filter | $1<z<2$ | $2<z<3$ | $3<z<4$
HST/ACSsS | F435W | $23.0^{23.6}_{22.1}$ | $23.5^{24.6}_{22.9}$ | $25.5^{26.4}_{24.9}$
HST/WFC3 | F160W | $21.7^{22.4}_{21.0}$ | $22.5^{23.1}_{21.9}$ | $23.3^{24.1}_{23.2}$
JWST/NIRCam | F115W | $22.0^{22.6}_{21.1}$ | $23.1^{23.9}_{22.5}$ | $23.8^{24.3}_{23.3}$
| F200W | $21.7^{22.4}_{20.9}$ | $22.5^{23.1}_{21.9}$ | $23.0^{23.7}_{22.8}$
| F356W | $21.4^{22.2}_{20.7}$ | $22.3^{22.8}_{21.7}$ | $23.0^{23.7}_{22.8}$

Note. Median values, first and third quartiles of the integrated magnitude
measured in different redshift bins and for several bands

####  3.2.2 2D Photometry

Separately, we also define a grid inside the circular aperture with cell size
equal to 3$\times$3 pixels ($0.18^{\prime\prime}\times 0.18^{\prime\prime}$)
and we measure SEDs for each of these spatial resolution elements. We only
consider cells inside the aperture which have at least one pixel that belongs
to the galaxy we are measuring. If there are any pixels in the cells belonging
to other galaxies (according to the segmentation map), the value of these
pixels is replaced by sky pixels by adopting the same Gaussian sky noise
distribution previously added to the images. Additionally, we only keep SEDs
for our analysis from cells with a SNR $>$ 3 in at least 5 bands and with
surface brightness brighter than 25 mag/arcsec2 in WFC3/F160W. Cells that do
not satisfy these conditions are discarded. The reason for imposing this
surface brightness limit of 25 mag/arcsec2 is to deal with the problem of the
larger-than-observed galaxy sizes in Illustris galaxies with
M${}_{\ast}\lesssim 10^{10.7}$ M⊙, which present larger half-light radii and
more extended discs than real observed galaxies (Snyder et al. 2015).

Fig. 4 shows an example of the integrated aperture and the grid for one of the
galaxies in the sample. We also include in this figure the segmentation map
around this galaxy and the initial integrated aperture of $r=2\times
r_{\text{hm}}$ (dotted line). We show the WFC3/F160W image for this galaxy
before and after replacing the values of the pixels of nearby galaxies by
random values drawn from the same Gaussian sky noise distribution previously
used to add the noise in this image. The initial grid covers all the region
within the final integrated aperture, but we only keep cells in cyan for the
stellar population analysis, as described in the previous paragraph.

Figure 4: Postage stamp images (size $2.5\arcsec\times 2.5\arcsec$) for one of
our galaxies: Illustris-1_066_0000006 at z = 2.192 (field A of S+17 images).
This galaxy code stands for galaxy id 6 and snapshot 66 ($z=2.21$) in the
Illustris-1 simulation. Top panel: Segmentation map, with colored regions
delimiting the galaxy (in orange), three nearby galaxies (green, yellow, and
red), and sky pixels (black). We also show the integrated photometric aperture
(dark cyan solid line) and the starting aperture of $r=2\times r_{\text{hm}}$
used to calculate the former (dotted line). Inside the integrated aperture, we
show the grid where the 2D photometry is measured: the cells that are finally
kept for the 2D-SPS analysis are shown in cyan, cells discarded for not
fulfilling the SNR and surface brightness criterion as pink hatched squares,
and cells discarded for not including any pixel from the considered galaxy as
white hatched squares. Bottom left panel: WFC3/F160W image (already
registered, PSF-matched, and sky noise-added) with the integrated aperture.
Bottom right: WFC3/F160W image where the values of the pixels from nearby
galaxies have been masked by replacing them with the same Gaussian sky noise
distribution previously added to the images. We also show the integrated
aperture and the grid.

The typical number of cells per galaxy decreases as we move to higher
redshifts, as expected given the smaller size of high-redshift galaxies (e.g.,
Bouwens et al. 2004, Oesch et al. 2010b, Ono et al. 2013). Our galaxies
present (median and quartiles) 98${}^{136}_{\phantom{1}69}$ cells at $1<z<2$,
52${}^{86}_{37}$ at $2<z<3$, and 24${}^{29}_{21}$ at $3<z<4$. We also
calculate the median SNR of the galaxy cells included in the analysis in
different redshift bins, by calculating the SNR as the ratio between the
measured flux and its uncertainty in these cells. Median and quartile values
for ACS/F435W (WFC3/F160W) are 31${}^{48}_{22}$ (30${}^{52}_{22}$), at
$1<z<2$, 30${}^{40}_{21}$ (27${}^{38}_{22}$) at $2<z<3$, and 17${}^{23}_{10}$
(24${}^{31}_{21}$) at $3<z<4$. For F115W, F200W and F356W NIRCam bands, the
median SNR of the galaxy cells are: 62${}^{93}_{51}$,
72${}^{110}_{\phantom{1}61}$ and 92${}^{132}_{\phantom{1}77}$, respectively,
at $1<z<2$, 43${}^{54}_{34}$, 66${}^{88}_{60}$ and
88${}^{110}_{\phantom{1}77}$ at $2<z<3$, and 53${}^{69}_{48}$,
81${}^{91}_{74}$ and 91${}^{108}_{\phantom{1}80}$ at $3<z<4$.

###  3.3 Redshifts

Galaxy redshifts are taken from the Illustris database. Specifically, we use
the $`$inferred redshift’ values in the catalogs that were used to generate
the simulated images in S+17. This redshift value corresponds to the inferred
cosmological redshift obtained considering the contribution of both the true
cosmological redshift and the galaxy peculiar velocity. A discussion about the
possible differences in the derived stellar population properties due to
uncertainties in the redshifts of the corresponding observed galaxies is
beyond the scope of this paper. Although photo-z uncertainties can be
substantial, it will be likely that future spectroscopic redshift samples from
both JWST and the Atacama Large Millimeter/submillimeter Array (ALMA) will
increase the fraction of galaxies with spectroscopic redshifts and, at the
same time, improve the photo-z accuracy at such high redshifts.

Fig. 5 shows the stellar mass vs. redshift plot for our galaxies. Two
different values of stellar masses are shown for each galaxy: considering all
particles in the galaxy or only particles inside twice the stellar half-mass
radius. The median values and quartiles for these mass measurements, in
log(M∗/M⊙), are $10.43_{10.17}^{10.73}$ and $10.30_{10.05}^{10.63}$,
respectively, and $2.03_{1.57}^{2.47}$ for the redshifts.

Figure 5: Stellar mass vs. redshift plot for galaxies in our sample, with
values extracted from the Illustris database. We show the stellar mass when
considering all the particles in the galaxy (total stellar masses) as light
blue circles and the stellar mass obtained when only particles inside twice
the stellar half-mass radius as small dark blue dots. Histograms of stellar
masses and redshifts are shown at the top and on the right, with the median
and quartiles marked.

###  3.4 Ground-truth Physical Properties of each Galaxy

We use the Illustris database to extract the information for all the simulated
particles belonging to a galaxy. In particular, we extract total stellar
masses and SFRs of the whole galaxy. These total stellar masses and SFRs
values are the ones shown in Fig. 1.

We build ground-truth galaxy SFHs from the information of the individual
particles that belong to each galaxy. The SFH can be computed by first loading
from the database the stellar particles in the galaxy and then making a
histogram of their formation ages in lookback time. We assume time bins for
the formation ages of 25 Myr up to a formation age of 1 Gyr, and 250 Myr
afterwards. Subsequently, the stellar mass formed per time interval is
calculated by summing the masses of the stellar particles formed in each time
bin. We will also use in the following sections an additional SFR estimation
for each galaxy obtained from the SFH and calculated by loading the gas
particles in the galaxy for the snapshot corresponding to the observed
redshift and adding together their instantaneous SFR values. Since these
ground-truth SFHs will be compared with those derived from the SED-fitting in
Section 5, in order to make the comparison more fair, we only consider stellar
and gas particles which are closer to the galaxy center than the radius of the
photometric aperture. Still, our SFHs will be naturally affected by the
limitation inherent to detection and photometry procedures.

##  4 ESTIMATION OF THE SFH FROM 2D SED FITTING

Our aim is to study the earliest formation phases of nearby massive galaxies
by analyzing in detail the formation history and location of the star
formation in massive galaxy progenitors at high redshift. For that purpose,
our approach consists in using broad-band data covering the (observed-frame)
optical-to-mid-infrared spectral range provided by HST and JWST. In this
paper, we assess the robustness of the results derived from this type of
realistic but naturally simplified analysis using simulated imaging data from
Illustris (see Section 3). As described in the previous section, we build SEDs
for each source for both the integrated emission as well as in a 2D grid. In
this section, we describe the derivation of integrated SFHs for each galaxy
based on SPS modeling for these SEDs. We then analyze how successfully we can
recover the ground-truth SFHs when only broad-band HST and JWST photometric
data is used.

###  4.1 Stellar Populations Synthesis Modeling

The integrated and grid SEDs are compared with stellar populations models of
Starburst99 (Leitherer et al. 1999), assuming a Kroupa (2001) IMF and a
Calzetti et al. (2000) attenuation law. To perform these fits, we use the
synthesizer code described in Pérez-González et al. (2003, 2008). This code
combines the emission of both stars and gas, compares the models with the
observed data, and returns the model that best fits the data by performing a
$\chi^{2}$ minimization.

We assume a double-burst delayed exponential SFH, with each burst described by
$SFR(t)\propto t\cdot e^{-t/\tau}$ for t $>0$ up to $t=t_{\mathrm{burst}}$,
being $\tau$ the star formation timescale and $t_{burst}$ the age of the
burst. The age of the burst must be understood as the time passed between the
age of the Universe when the galaxy started to form stars in that burst and
the time corresponding to the observed redshift. The SFH form was chosen after
some testing which indicated that 2 bursts following a delayed-exponential SFH
instead of other simpler parametrizations (e.g., one single exponential burst)
as a more adequate parametrization to successfully reproduce the ground-truth
SFHs given by Illustris.

The free parameters in our SED fits are the age ($t$), star-formation
timescale ($\tau$), extinction ($A_{V}$), and metallicity ($Z$) for the two
stellar populations, in addition to the burst strength ($b$). This burst
strength describes the fraction of the total stellar mass that has been
created by the most recent burst. We consider an old burst as the one
occurring first in the galaxy formation history and a young burst as the one
occurring closer to the age of the Universe corresponding to the observed
redshift. The stellar mass (M∗) is derived from the SED fits by normalizing
the best-fitting model (which provides mass-to-light ratios at all
wavelengths) to the observed photometry.

In order to improve the resemblance of the estimated SFHs to the ground-truth
SFHs built from simulated particles in the simulation, we found that the
allowed age ranges for the old and young stellar populations must depend on
the galaxy redshift. In other words, the age frontier between the two stellar
populations is an important parameter to set a priori based on the galaxy
redshift. We thus tested the dependence of our results on this age separation
limit (age${}_{\text{lim}}$ in Table 3), considering values as a function of
the age of the Universe for a given redshift. The most accurate results in our
analysis (considering the mass-fraction formation ages discussed in Section 5)
are obtained when we impose an age limit between the 2 populations equal to
40% of the age of the Universe for galaxy redshifts at $1<z<2$ and to 50% of
the age of the Universe for $2<z<4$ (i.e., 40% of the age of the Universe for
the young population and 60% for the old population when $1<z<2$, and 50% of
the age of the Universe for both populations when $2<z<4$).

Similarly, we tested different values for the star formation timescales of
both bursts and found the best results are obtained when setting 200
$<\tau_{\text{old}}<$ 1000 Myr, and $\tau_{\text{young}}=$ 10 Myr. Table 3
shows the ranges within which each free parameter is allowed to vary in the
SED-fittings of this work.

It is impossible to present the results of the full set of tests performed to
explore the parameter space here, so instead we focus on those that yield the
best match and will be the most useful for upcoming JWST observations.

Table 3: Free parameters and their ranges in the SED fitting for a double
delayed-exponential SFH

Parameters | Values/Range |  ssUnitsss | Step
---|---|---|---
$t_{\text{young}}$ |  suss.0.1 - age${}_{\text{lim}}$aafootnotemark:  | Gyr | discretebbfootnotemark:
$\tau_{\text{young}}$ | 10 | Myr | –
$A_{\text{V,young}}$ | 0 - 2 | mag | 0.1 mag
$t_{\text{old}}$ |  ssssiage${}_{\text{lim}}$ \- age${}_{\text{Univ,z}}$ccfootnotemark:  | Gyr | discretebbfootnotemark:
$\tau_{\text{old}}$ |  i.200 - 1000 | Myr | 0.1 dex
$A_{\text{V,old}}$ | 0 - 1 | mag | 0.1 mag
$Z_{\text{young}}$ & $Z_{\text{old}}$ | [0.2, 0.4, 1] . | Z/Z⊙ | discrete
$b$ | 0.01 - 1 s.. | - | 0.01

Notes. a afootnotetext: The age separation limit between the young and old
populations, is measured (backwards) from the redshift of observation of the
galaxy. This value depends on the redshift of the galaxy: it is set to 40% of
the age of the Universe for 1 $<$ z $<$ 2 and to 50% of the age of the
Universe for 2 $<$ z $<$ 4. b bfootnotetext: We use all the discrete values
for the ages given by the SB99 models within the allowed range. c
cfootnotetext: The maximum allowed value for the age of the old population is
the age of the Universe at the galaxy redshift.

Figure 6: Example of one integrated SED fit for one of our galaxies:
Illustris-1_062_0004079 at $z=2.843$. This galaxy code stands for galaxy id
4079 and snapshot 62 ($z=2.73$) in the Illustris-1 simulation. The best fit is
shown as a red line. We also show the models corresponding to the old (green)
and young (orange) stellar populations. The transmission curves of the filters
have been included at the bottom of the figure. The best-fit stellar
parameters for the old and young populations are also given. An RGB image
(size $1.5\arcsec\times 1.5\arcsec$) with the integrated aperture is shown as
an inset.

The metallicity is allowed to adopt three different values: Z/Z⊙ = 0.2, 0.4,
and 1. These values are expected for our sample according to the mass-
metallicity relationship at the considered redshifts, at which this
relationship shows lower metallicities than locally for a given mass (Erb et
al. 2006; Maiolino et al. 2008; Mannucci et al. 2009; Zahid et al. 2011).
Regarding the dust attenuation, we limit the attenuation range to values
within $0<A_{\text{V,old}}<1$ mag for the old population, and within
$0<A_{\text{V,young}}<2$ for the young population.

Fig. 6 shows an example of an integrated SED fit and the values of the model
that best fits our data. For this particular galaxy, both the old and young
stellar populations contribute at the 40%-60% level throughout the full
spectral range. This analysis is performed for the 221 galaxies: both for
their integrated measurements and their 2D emission (given by the cells in the
grid).

Figure 7: Stellar masses for the $1<z<4$ sample derived from the integrated
SED-fits vs. their ground-truth values calculated from the simulated stellar
particles belonging to each galaxy (inside a sphere with the same radius as
that of the integrated photometric aperture). The one-to-one relation is shown
with the dotted black line. The histograms for both distributions are shown in
different colors, with the median and the 68% intervals as horizontal
segments.

To estimate uncertainties in the derived stellar population properties, Monte
Carlo (MC) simulations are performed in synthesizer by allowing the
photometric data to normally vary within their photometric uncertainty
(without correlation), and then repeating the fit again for 300 resampled SEDs
(more details in, e.g., Domínguez Sánchez et al., 2016). This results in 300
sets of solutions for each SED which also provide us with information about
the typical degeneracies present in these kinds of studies (e.g., age-
metallicity or $\tau$-age degeneracy).

Fig. 7 shows the stellar masses for our $1<z<4$ sample derived from the
integrated SED-fits vs. their ground-truth values and how the one-to-one
relation is recovered, with a median offset of 0.04 dex and a scatter of 0.2
dex. Each of the SPS-derived galaxy masses has been calculated as the median
stellar mass provided by the integrated SED fits. The ground-truth masses have
been calculated from the database by considering only the simulated particles
of each galaxy closer to the galaxy center than the radius of the integrated
photometric aperture. In both cases, stellar masses correspond to the current
mass of stars at the redshift of observation (i.e., calculated after taking
into account the time-dependent mass-loss and fraction of remnants as a
function of time). The histograms for both distributions are also shown. We
notice our SPS-derived stellar masses slightly overestimate the ground-truth
masses, although the differences are small:
log(M∗/M${}_{\odot})=10.33_{\phantom{1}9.96}^{10.80}$ for the integrated SED-
fits vs. $10.21^{10.72}_{\phantom{1}9.92}$ for the stellar particles (median
and 68% confidence interval).

###  4.2 Estimating SFHs from HST+JWST Photometry

Figure 8: Schematic diagram of the methodology followed to obtain the median
SFH for a galaxy from its 2D-SPS analysis. First, we show an RGB image of the
galaxy and some of the HST+JWST broad-band images (in Table 1) processed to
imitate CANDELS and CEERS observations for HST and JWST filters, respectively.
We measure the 2D photometry inside a grid (in cyan) with cell size equal to
the spatial resolution element, in addition to the integrated photometry (blue
aperture). As an example, we show the SEDs measured for three regions of the
grid: the center of the galaxy (teal cell), an arm region (light green), and a
diffuse emission zone (pink). To estimate the uncertainties and degeneracies
in the derived stellar populations, each SED is fitted 300 times by performing
MC simulations. We add the individual SFHs inferred from all the grid regions
to obtain the SFH for the whole galaxy. The black thin SFHs in the bottom
subfigure show the SFHs for the whole galaxy created by accounting for the
uncertainties in the SFHs of each grid region, where the median galaxy SFH
from this 2D-SPS method is shown in yellow. The blue dotted SFH is the galaxy
SFH inferred from the integrated photometry. The ground-truth SFH of the
galaxy, given by the stellar particles belonging to the galaxy, is the green
dash-dotted SFH.

Integrated SFHs are built for each galaxy from the 2D SED-fits, and compared
with the ground-truth SFHs extracted from the Illustris database (see Section
3.4). We build 300 integrated SFHs for each galaxy using the 2D stellar
population fits, which includes a MC method to estimate uncertainties and
degeneracies. We first create one SFH for each resolution element in the
galaxy grid. Subsequently, one global SFH is calculated by adding together all
these SFHs in the grid. This procedure is then repeated 300 times, one for
each of the MC particles presented in the previous section. By combining and
adding the 300 SED-fits for each cell in the grid, we obtain 300 independent
estimations of the integrated SFH for a given galaxy.

The resulting SFHs are smoothed using a 100 Myr square kernel. Then, we
normalize these 300 galaxy SFHs with the median stellar mass of the galaxy
derived from the 300 integrated SED-fits. Finally, we calculate the median SFH
of the galaxy from these 300 normalized SFHs.

In Fig. 8, we summarize the methodology followed for a galaxy to obtain this
2D-SPS median galaxy SFH from its 2D photometry measurements. We start from
the HST+JWST data (RGB image and some individual bands shown on the figure),
in our example, a galaxy presenting a red center (a protobulge?) and what
seems like a blue disk with some spiral structure. We measure integrated and
2D photometry in a grid. We show SEDs and their stellar population modeling
results for some representative grid regions in the center, blue arm on the
bottom left and a diffuse emission zone on the top right. We note the
difference in color shown in the SEDs between the center and the disk, and how
the integrated SED resembles more that of the galaxy center rather than the
disk regions. In this regard, considering 2D photometry in our analysis
facilitates the estimation of SFHs since small parts of a galaxy are expected
to have more simple forms than the entire galaxy, whose integrated photometry
is dominated by the regions with highest intensities, i.e., they overshine
fainter regions located in the outskirts for this example galaxy. The
individual SFHs for each grid region are added to obtain an integrated SFH. In
this example, we demonstrate how the derived SFHs with our method nicely
reproduce the ground-truth SFH, significantly improving what can be obtained
by analyzing the integrated SED. In the following section, we compare our
derived SFHs from our 2D analysis of broad-band imaging data with the
reference SFHs built using the Illustris-1 database in a statistical way.

##  5 VALIDATION OF THE METHOD

In this section, we evaluate the levels of success of our method in recovering
the SFHs of massive galaxies at high redshift. In particular, in this paper we
are mainly interested in the earliest phases of massive galaxy formation.
Therefore, we compare our results based on the analysis of broad-band HST and
JWST images with the ground-truth provided by the Illustris database for full
galaxies. We will analyze spatially resolved stellar population properties in
a future paper, here we concentrate on the results about the SFHs.

###  5.1 Characterization of the Earliest Phases in the Formation of Massive
Galaxies

In order to characterize when massive galaxies start their formation, we
calculate the formation times and redshifts when galaxies have formed a given
(small) fraction of their total stellar mass. In particular, we discuss here
the formation times when the first 5%, 10% and 25% of the stars in each galaxy
and the entire sample were formed. We define $t_{\text{k}}$ as the formation
time (measured from the galaxy redshift) at which a galaxy formed the k$\%$ of
its total stellar mass present at the galaxy redshift. These mass-fraction
formation times, $t_{\text{k}}$, can be directly calculated from the SFHs
derived with our 2D-SED fitting method and are similar to other parameters
used in recent literature (see, e.g., Ji & Giavalisco, 2022). We remark that
$t_{\text{k}}$ depend on both $t_{\text{0}}$ and $\tau$, so they might somehow
alleviate possible degeneracies between ages and timescales. We also calculate
mass-weighted ages, $t_{\text{mw}}$, from the SFHs. To test their robustness,
these $t_{\text{5}}$, $t_{\text{10}}$, $t_{\text{25}}$, and $t_{\text{mw}}$
extracted from the SFHs derived with our 2D-SPS method are compared with their
ground-truth values, which are obtained from the SFHs of the simulation
stellar particles following the same procedure and enclosed by the photometric
aperture.

Fig. 9 shows $t_{\text{5}}$, $t_{\text{10}}$, $t_{\text{25}}$, and
$t_{\text{mw}}$ calculated from the SFHs derived with 2D stellar population
modeling of broad-band data versus their reference values. For each galaxy, we
show 300 values of these quantities (vertically spread out), which correspond
to the 300 MC particles or 2D-SPS galaxy SFH as described in Section 4.2). The
black dotted line shows the one-to-one relation between values derived from
our 2D-SPS method (output) and those from the simulated particles (ground-
truth). As explained in Section 4.1, some a priori parameters of our method
(e.g. Table 3) have been optimized to reproduce this one-to-one relation for
these four mass-fraction formation times. In general, we find our ages are
consistent with this relation at all redshifts. As an inset in each panel, we
show the histogram of relative differences between median and ground-truth
values for all galaxies. We find our $t_{\text{5}}$, $t_{\text{10}}$,
$t_{\text{25}}$, and $t_{\text{mw}}$ values are consistent with the ground-
truth values with a median offset of 0.019, 0.0076, -0.0006, and -0.023 dex,
and a scatter (68% interval) of 0.21, 0.24, 0.28, and 0.26 dex, respectively.
In general, we do not see any systematic effect as a function of redshift.

Regarding the precision of our method, we find the scatter values are
relatively small and similar for all $t_{\text{k}}$, but we observe the
scatter tends to increase from $t_{\text{5}}$ to $t_{\text{25}}$. Taking all
this into account, we conclude that our 2D-SPS method successfully recovers
$t_{\text{5}}$, $t_{\text{10}}$, $t_{\text{25}}$ and $t_{\text{mw}}$ for
$1<z<4$ massive galaxies with a $\sim$25% uncertainty and small $<$3%
systematic effects.

Figure 9: Comparison between the mass-fraction formation times (measured
backwards from the redshift of observation of the galaxy) of our 2D-SPS method
and their ground-truth values: $t_{\text{5}}$ (upper-left), $t_{\text{10}}$
(upper-right), $t_{\text{25}}$ (lower-left) and $t_{\text{mw}}$ (lower-right)
calculated from our 2D-SPS-derived SFHs versus their reference values
calculated from the SFHs built from the simulated particles in galaxies. Each
galaxy is represented by 300 vertically spread points, which correspond to the
mass-fraction formation times of the 300 2D-SPS SFHs, built from the 300 MC
particles in each resolution element in the grid. The median of these values
for each galaxy is shown as bigger circles and error bars represent the
standard deviation of these values (68% interval). All points are color-coded
by galaxy redshift. As an inset, we include the offsets of the galaxies in the
sample.

##  6 EXPECTATIONS FOR THE DERIVATION OF THE SFH OF $z>1$ MASSIVE GALAXIES
WITH HST+JWST DATA

The aim of this section is to discuss when the early stages of stellar mass
assembly took place in very massive galaxies. To do this, we apply stellar
population synthesis in 2D on our galaxy sample at $1<z<4$ and we compare the
statistical results on the SFHs for our whole sample with the ground-truth
values inferred from the simulation stellar particles. This comparison gives
us information about the limitations and observational biases we will
encounter when using this method on galaxy samples constructed with real data.

###  6.1 When Did Massive Galaxies Begin to Form?

In Section 5, we showed that our 2D-SPS method successfully recovers
$t_{\text{5}}$, $t_{\text{10}}$, and $t_{\text{25}}$ for our galaxy sample.
These quantities represent the lookback times at which galaxies reach 5%, 10%,
25% of the stellar mass formed, respectively, and can be used to estimate the
beginning of star formation in galaxies. To address the question of when the
population of very massive galaxies began to form, our approach consists in
calculating $t_{\text{5}}$, $t_{\text{10}}$, and $t_{\text{25}}$ from a median
SFH of the sample built from the 2D-SED fits. At this point, we should remind
the reader that our sample consists of galaxies at $1<z<4$ with
M${}_{\ast}>10^{10}$ M⊙, all having very massive descendants
(M${}_{\ast}>10^{11}$ M⊙) at $z=0$. Hence, we expect the median SFH from our
2D-SPS modeling to resemble the median SFH of very massive galaxies at $z=0$
over the redshift interval of our sample (where both SFHs overlap). We say
resemble because our sample limitations might have implications on the derived
SFHs, i.e., we have to consider the systematic effects introduced by not
taking into account lower mass galaxies and by the relative number of galaxies
selected at different redshifts in our full interval.

Since the final goal of our 2D-SPS method is to be applied on real galaxy
samples at high redshift as soon as JWST data are available, we need to check
first whether our results from the 2D-SED fits regarding the first episodes of
stellar assembly in very massive galaxies are compatible with the ground-truth
of the simulation. To address this issue, we build two additional samples of
galaxies at $z=0$ in the simulation to be compared with our main $1<z<4$
galaxy sample. For these three samples of galaxies, we build the typical
(median) ground-truth SFH from the simulated particles in their galaxies,
which will be compared with the typical SFH from the 2D-SPS on the $1<z<4$
studied sample. The three samples of galaxies in the Illustris-1 simulation
considered in this section are:

  * •

Our main galaxy sample at $1<z<4$ used for the 2D-SPS analysis. It consists of
the 221 (out of the 248) M${}_{\ast}>10^{10}$ M⊙ progenitors at $1<z<4$ in the
S+17 images with a very massive descendant at $z=0$, and which have not been
discarded during the analysis procedure (see Section 2.2). This is the
“HST$+$JWST high-redshift sample” or “2D-SPS sample”, hereafter.

  * •

The descendants at $z=0$ of our main high-redshift galaxy sample, all of them
with M${}_{\ast}>10^{11}$ M⊙ (132 descendants). This sample is a subset of the
whole population of very massive galaxies at $z=0$.

  * •

The whole population of very massive galaxies (M${}_{\ast}>10^{11}$ M⊙) at
$z=0$. This sample is composed by 856 galaxies.

We remind the reader that the S+17 images only contain a limited number of
galaxies from the full Illustris-1 simulation. Thus, the descendants at $z=0$
of the galaxies we have studied through their HST+JWST simulated imaging data
are also a limited subset of all the galaxies at $z=0$ in the simulation. Our
aim is to see if, by analyzing the SFHs for high-redshift,
M${}_{\ast}>10^{10}$ M⊙ progenitor galaxies in the images of $z=0$ very
massive galaxies, we can learn when the whole population of very massive
galaxies began to form.

If we consider the masses of all the simulated stellar particles in galaxies,
our main $1<z<4$ sample of 221 galaxies accounts for only 28% of the total
stellar mass present in their 132 $z=0$ descendants. If we also take into
account the total stellar masses of less massive progenitors at $1<z<4$ in the
images, none of them considered by the mass cutoff, this number raises to 32%.
This means that the remaining 68% of the stellar mass in nearby massive
galaxies must be explained by either more recent in-situ star formation events
or by subsequent mergers at $z<1$. We refer the reader to Section 2.2 for more
information about the main $1<z<4$ galaxy sample and their $z=0$ descendants.

In Fig. 10 we compare the histograms for the ground-truth stellar masses of
galaxies in the three samples, calculated from the simulated particles in the
database. On different panels we show the histograms for the stellar masses
corresponding to all particles in galaxies and to only particles inside
2$\times$r${}_{\text{hm}}$. In both cases, we find that the distribution of
masses for the population of very massive galaxies at $z=0$ and for the
specific descendants are very similar: the median and quartile values are
log(M∗/M⊙)= $11.3_{11.1}^{11.7}$ vs. $11.2_{11.1}^{11.6}$ when all particles
are considered, and $11.1_{11.0}^{11.5}$ vs. $11.1_{11.0}^{11.4}$ for
particles inside 2$\times$r${}_{\text{hm}}$. As discussed below in this
section, this difference in the median masses is a consequence of the mass-
cutoff imposed for the selection of the $1<z<4$ sample and will have an impact
on the mass-fraction formation times of both samples. In summary, we conclude
that our limited sample would provide representative results for the full
population of massive galaxies at $z=0$.

Figure 10: Ground-truth stellar mass histograms for the three subsamples
extracted from the Illustris-1 database: when considering all particles in
galaxies (top panel) or only particles inside twice the stellar half-mass
radius (bottom panel). We show our main 2D-SPS sample at $1<z<4$ in red, their
descendants at $z=0$ in orange, and the whole population of very massive
galaxies at $z=0$ in dark gray. The median and the 68% intervals are shown at
the top of each panel.

In subsection 6.1.1, we calculate the typical ground-truth SFH for these three
subsets of galaxies using only the simulated particles in them, while in
subsection 6.1.2 we compare these results with the median SFH of our main
high-redshift sample of M${}_{\ast}>10^{10}$ M⊙ precursors built using our
2D-SPS method. Finally, in subsection 6.1.3 we show the variation of the onset
of star formation of massive galaxies.

####  6.1.1 Ground-truth SFHs: very Massive $z=0$ Galaxies in Illustris and
the $1<z<4$ Sample

Fig. 11 shows the sample-averaged (median) SFHs built from the Illustris
simulated particles for the three different subsets of galaxies considered:
the $1<z<4$ sample (blue line), the descendants at $z=0$ of these high-
redshift galaxies (orange lines), and all massive galaxies at $z=0$ in
Illustris-1 (green lines). To study whether there is any effect of the galaxy
aperture on the SFHs, the SFHs of the two subsamples of $z=0$ galaxies have
been built by selecting simulated particles in two different ways: either by
considering all the stellar particles belonging to each galaxy (dotted lines),
or only the stellar particles inside a sphere with radius 2$\times
r_{\text{hm}}$ (measured at $z=0$; solid lines). In the case of the HST$+$JWST
high-redshift sample, the radius of each sphere has been set to match the
radius of the photometric aperture used for the 2D-SED fits, rphot, and only
particles within this radius have been considered. These typical SFHs have
been smoothed using a square 250 Myr kernel. Shaded areas correspond to the
uncertainty in the median SFH calculated as the 95% confidence interval:
$t_{N-1}\times\sigma/\sqrt{N}$, where $t$ is the $t$ value from the
$t$-distribution for 95% confidence, $\sigma$ has been estimated from the
dispersion of the different SFHs, and $N$ is the number of galaxies at each
age. We only show median SFHs of the galaxies at $z=0$ down to $z=1$, the
lower limit of our HST$+$JWST sample. Additionally, we point out that the
typical SFH of the HST$+$JWST $1<z<4$ sample has been calculated from
different numbers of galaxies at different ages (what we would be able to
perform when we have real data), in contrast to the typical SFHs of the other
two $z=0$ subsamples, which count with the same number of galaxies in the
entire redshift range.

Regarding the different levels of the SFHs in Fig. 11, these are a consequence
of the fact that no normalization has been applied to them. Therefore, these
differences in the SFR values can be explained by the construction method of
the SFH followed in each case: SFHs which have been built using all particles
present generally higher SFRs than those built using only a limited number of
particles. Both for the population of $z=0$ $10^{11}$ M⊙ galaxies and for the
$z=0$ descendants of the $1<z<4$ sample, the stellar mass inferred from the
typical SFHs of particles within $2\times$r${}_{\text{hm}}$ (i.e., the area
enclosed by the SFH) is $\sim 66$% of the stellar mass inferred from the
typical SFH built from all the particles in galaxies. In addition, the average
SFH of the $1<z<4$ sample is the one that has the lowest values of SFR. This
is expected, since the radii of the photometric apertures used to build the
sphere that contains the simulated particles in the case of the 2D-SPS sample
are usually smaller than 2$\times$r${}_{\text{hm}}$ measured at the observed
redshift, as explained in Section 3.2, and should be smaller than
2$\times$r${}_{\text{hm}}$ measured at $z=0$. In fact, we remind the reader
that the typical photometric aperture radius is nearly 20% smaller than
$2\times$r${}_{\text{hm}}$ and enclose around 65% of the total stellar mass in
the galaxies (see Fig. 3).

Figure 11: Ground-truth median SFHs calculated from the simulated particles
for the subsamples of galaxies: all M${}_{\ast}>10^{11}$ M⊙ galaxies at $z=0$
(green lines), our HST$+$JWST $1<z<4$ sample (blue), and the descendants at
$z=0$ of this high-redshift sample (orange). The results from the 2D-SPS
analysis are not included in this figure. The median SFHs for the $z=0$
subsets have been calculated in two ways: considering all the particles in
galaxies (dotted lines) or only the particles inside a sphere with radius
2$\times$r${}_{\text{hm}}$ to the galaxy center (solid). The median SFH of the
high-redshift sample has been built using only particles whose distance to the
galaxy center is lower than the radius of the photometric aperture used for
each galaxy. Shaded areas represent the uncertainty of the median. As vertical
lines, we show the $t_{\text{5}}$, $t_{\text{10}}$, $t_{\text{25}}$, and
$t_{\text{mw}}$ calculated for each SFH shown. As a comparison, we show the
median SFH calculated in Iyer et al. (2020) for all galaxies in the Illustris
simulation with $10^{10.75}<$ M${}_{\ast}<10^{11.25}$ M⊙.

In Fig. 11, we also show with vertical lines the mass-fraction formation times
$t_{\text{5}}$, $t_{\text{10}}$, $t_{\text{25}}$, and $t_{\text{mw}}$
calculated for each sample-averaged SFH. As mentioned before, we only take
into account for these calculations the stellar mass formed in each SFH at
ages of the Universe below that of the minimum redshift of the 2D-SPS sample
($z\sim 1$), as at lower redshifts we lack information on the photometric
properties of any of the galaxies in this sample-averaged SFH. This means we
only consider the at $z>1$ part of the median SFHs for the $z=0$ galaxy
subsets and the whole SFH for the high-redshift sample. We find the mass-
fraction formation times for the SFH of the HST$+$JWST high-redshift sample
are systematically earlier (expressed in terms of the age of the Universe; or
correspond to older lookback times), than the ones for the descendants of this
sample at $z=0$, and these, in turn, are younger than those from the whole
population of massive galaxies at $z=0$. These $200-400$ Myr shifts ($\sim
10$% in relative terms) could be explained by a progenitor bias, which appears
as a consequence of the mass cutoff of M${}_{\ast}>10^{10}$ M⊙ imposed on our
2D-SPS sample selection. We first focus on the younger (or earlier) formation
epochs of the descendants with respect to those of all the very massive
galaxies at $z=0$. The mass cutoff makes the sample of $z=0$ descendants to be
biased towards larger masses and more massive galaxies tend to present older
stellar population ages. In contrast, when taking into account all
M${}_{\ast}>10^{11}$ M⊙ galaxies at $z=0$, i.e., a complete sample of local
massive galaxies, we obtain slightly later formation ages. In fact, as it can
be seen in Fig. 10, the median and quartile stellar masses are $0.1-0.2$ dex
larger for the descendants of our sample of M${}_{\ast}>10^{10}$ M⊙ galaxies
at $1<z<4$ compared to the complete sample of M${}_{\ast}>10^{11}$ M⊙ $z=0$
galaxies. In addition, our mass cut at $1<z<4$ implies losing 17% of massive
M${}_{\ast}>10^{11}$ M⊙ galaxies at $z=0$.

Indeed, this progenitor bias also explains why the time-averaged SFR for $z=0$
descendants is higher than that for the population of very massive galaxies.
This applies to the galaxy apertures considered in Fig. 11 (for which no mass-
normalization has been applied). The differences increase for ages
corresponding to smaller fraction of the total stellar mass. Furthermore, both
for all massive galaxies at $z=0$ and for $z=0$ descendants of the high-
redshift sample, mass-fraction formation times are older (smaller values in
age of the Universe) when all stellar particles are included in the SFH
computation than when considering only those inside
$2\times$r${}_{\text{hm}}$. This may be explained by an outward migration of
stars, which would result in larger apertures adding a larger fraction of
older stars in the SFH computation and, consequently, earlier formation times.

Regarding the shift to earlier formation times of the 2D-SPS sample with
respect to the $z=0$ descendants, possibly also due to the progenitor bias, we
propose a similar explanation as the one given above. Our 2D-SPS sample is
only a biased subset (cut in mass and with a given redshift distribution) of
all the progenitors that evolve to a galaxy from the whole sample of $z=0$
descendants. The descendants, as mentioned above, also have other progenitors
that do not fulfill the mass cutoff at $z>1$ and which may have undergone a
merger with a galaxy from the 2D-SPS sample at a lower redshift. These minor
progenitors at $1<z<4$ would not be included in the SFH of the 2D-SPS sample,
and they would probably host younger stellar populations, which would explain
the shift towards more recent formation times for their $z=0$ descendants.

####  6.1.2 Recovering the SFH of Massive Galaxies with 2D-SPS Modeling of
HST+JWST Data

We now calculate the median SFH of the $1<z<4$ galaxy sample obtained from the
2D-SED fits (not from the Illustris database, as done in Section 6.1.1). For
this calculation, we take the SFH for each galaxy in this sample constructed
from the 2D-SPS analysis (see Section 4.2) and combine them to build the
typical SFH for the whole sample. We remark that this implies that the SFR at
a given age of the Universe involves the combination of a different number of
galaxies, namely, all that lie at lower redshifts compared to the redshift
corresponding to that age of the Universe.

Fig. 12 shows this median SFH of the $1<z<4$ sample derived from the 2D-SED
fits (red solid line). We compare this SFH, which could be obtained following
the same procedure on real galaxy samples when JWST (plus HST) data is
available, with the ground-truth. First, we compare with the averaged SFH for
the same 2D-SPS sample but built from the simulated particles inside the
radius of the photometric aperture (blue). This comparison allows to evaluate
the accuracy of our method. We also compare with the SFH of $z=0$ descendants
using the stellar particles inside 2$\times$r${}_{\text{hm}}$ (orange). This
comparison allows us to understand the bias of a sample based on imaging data
linked to selection effects. The two latter comparison SFHs are the same shown
in Fig. 11, but normalized to recover the same stellar mass as the median
2D-SPS-derived SFH ($10^{10.57}$ M⊙ at $z=1$). Note that, without
normalization, the corresponding stellar masses at $z=1$ would be
$\log($M∗/M${}_{\odot})=10.67$ for the median ground-truth SFH of this high-
redshift sample and 10.91 for the median SFH of their descendants at $z=0$
(built from particles inside 2$\times$r${}_{\text{hm}}$). Shaded areas
correspond to the uncertainty in the median SFH calculated as the 95%
confidence interval.

Figure 12: Comparison between the SFH of our main high-redshift sample derived
from the 2D-SED fits (in red) and from the simulated particles (in blue). We
also include the SFH of the descendants of the high-redshift sample at $z=0$
built with particles inside 2·r${}_{\text{hm}}$ (in orange). Shaded areas
represent the uncertainty of the median. The SFHs have been normalized to
recover the same stellar mass as the median 2D-SPS-derived SFH over the same
redshift intervals ($z\gtrsim 1$). The vertical lines depict the
$t_{\text{5}}$, $t_{\text{10}}$, $t_{\text{25}}$, and $t_{\text{mw}}$ mass-
fraction formation times for each SFH. The top panel shows the evolution of
the stellar mass assembly calculated from the median SFHs of the high-redshift
sample (2D-SPS and ground-truth, same color code). Gray diamonds show the
Illustris SFMS level at $z=1$, 2, and 4 (Sparre et al. 2015). On the top
panel, we show the evolution of the integrated stellar mass at each redshift
for the HST$+$JWST high-redshift sample (2D-SPS and particles, same color
code).

We find that the general shapes of the three SFHs shown are not very different
from each other, except for the SFR increase near $z=1$ (in the last 0.5-1.0
Gyr) in the SFH built from the 2D-SPS analysis. This SFH of 2D-SED fits
follows the general trend of the SFH built from the simulated particles in the
sample, but for this raise. This discrepancy might be explained by the fact
that the number of galaxies considered in the calculation of these SFHs drops
as the redshift decreases (as clearly seen in the increase of uncertainties).
Thus, there are fewer galaxies near $z=1$ from which to calculate the typical
SFH for the high-redshift sample. In addition to this, as explained in Section
4.1, the individual galaxy SFHs built from the 2D-SPS are the combination of
the SFH of one young and one old stellar population. This causes the galaxy
SFHs built from the 2D-SED fits to normally have a peak in SFR near the
redshift of the galaxy, which is not necessary present in the galaxy SFH built
from the simulated particles (see, for example, the 2D-SPS median galaxy SFH
shown in Fig. 8). When the number of galaxies considered to build the SFH of
the sample is high, these peaks in the individual galaxy SFHs are not
reflected in the median SFH. But they do lead to a rise in the SFR value as
the number of galaxies drops.

In Fig. 12, we include the SFR values expected from the Illustris SFMS at
$z=1$, 2, and 4 (Sparre et al. 2015) for the stellar masses obtained by
integrating our median 2D-SPS SFH down to each of those redshifts. We find
that our median SFR values (in red in Fig. 12) are consistent within the
errors with those expected from the Sparre et al. SFMS. On the top panel, we
show the evolution of stellar mass assembly for the HST$+$JWST high-redshift
sample, calculated both from the median SFH of the 2D-SPS analysis and that of
the stellar particles.

The mass-fraction formation times for each SFH in the figure are shown as
vertical lines. We remark the good agreement between the $t_{5}$ and $t_{10}$
formation times derived for the high-redshift sample using the 2D-SED fits and
the simulated particles (differences of -0.009 and -0.002 dex in age of the
Universe, respectively, with the ground-truth values being larger).
Additionally, the differences in the other two formation times, although
higher, are also low: $+$0.02 dex for both $t_{25}$ and $t_{mw}$. Regarding
the comparison of these formation times for the high-redshift sample (2D-SPS
and particles) with the $z=0$ descendants, the agreement is also high, with
maximum differences in age of the Universe of -0.06, -0.04, -0.04 and -0.02
dex for $t_{5}$, $t_{10}$, $t_{25}$ and $t_{mw}$, respectively. As mentioned
before, the older values of the $t_{5}$, $t_{10}$ and $t_{25}$ formation times
for the high-redshift sample with respect to the $z=0$ descendants are due to
the progenitor bias.

Table 4: Mass-fraction formation times and redshifts for different samples

Galaxy subset | SFHs from | $t_{\text{5}}$ | $t_{\text{10}}$ | $t_{\text{25}}$ | $t_{\text{mw}}$ | $z_{\text{5}}$ | $z_{\text{10}}$ | $z_{\text{25}}$ | $z_{\text{mw}}$
---|---|---|---|---|---|---|---|---|---
Main $1<z<4$ sample of | 2D-SED fits (r${}_{\text{phot}}$) |  .S$1.4_{1.1}^{1.7}$ |  .S1.7${}_{1.5}^{2.1}$ |  .S$2.5_{2.2}^{3.1}$ |  .S$3.6_{3.2}^{4.1}$ |  .S$4.5_{3.8}^{5.2}$ |  .S$3.7_{3.1}^{4.2}$ |  .S$2.7_{2.2}^{2.9}$ |  .S$1.9_{1.6}^{2.1}$
M${}_{\ast}>10^{10}$ M⊙ progenitors | Particles | in r${}_{\text{phot}}$ | 1.40 | 1.73 | 2.39 | 3.39 | 4.41 | 3.70 | 2.78 | 1.98
All galaxies at $z=0$ | Particles | all | 1.60 | 1.98 | 2.76 | 3.69 | 3.97 | 3.29 | 2.43 | 1.81
with M${}_{\ast}>10^{11}$ M⊙ | in 2$\times$r${}_{\text{hm}}$ | 1.66 | 2.05 | 2.83 | 3.74 | 3.85 | 3.19 | 2.37 | 1.78
Descendants at $z=0$ | Particles | all | 1.52 | 1.88 | 2.60 | 3.58 | 4.14 | 3.45 | 2.57 | 1.87
of $1<z<4$ sample | in 2$\times$r${}_{\text{hm}}$ | 1.56 | 1.91 | 2.62 | 3.59 | 4.04 | 3.41 | 2.56 | 1.86

Note. Values have been calculated from the different SFHs shown in figs. 11
and 12. The mass-fraction formation times are measured in Gyr from Big Bang

According to the 2D-SPS analysis, local M${}_{\ast}>10^{11}$ M⊙ galaxies would
have formed 5% of their stellar mass present at $z\sim 1$ by $z=4.5$, 10% of
their mass by $z=3.7$, and 25% by $z=2.7$. These redshifts correspond to ages
of the Universe of 1.4, 1.7, and 2.5 Gyr, respectively. This is equivalent to
saying that 5% of the stellar mass present in very massive local galaxies at
$z\sim 1$ was already assembled when the Universe was only $\sim 10$% of its
current age, 10% of the mass by 13% of its age, and 25% of the mass by 18% of
the cosmic time.

The first stellar formation episodes for these galaxies are located at
$z=16\pm 1$ (calculated as the median redshift where the SFR starts to be
larger than 0). Then it would rise up to $z\sim 3$. After that, it would
remain approximately constant down to $z\sim 2$, and a slightly decreasing
trend would follow. This can be clearly seen in the ground-truth SFH derived
from the simulated particles in the database and agrees with the median SFHs
calculated in Iyer et al. (2020) for $10^{10.75}<$ M${}_{\ast}<10^{11.25}$ M⊙
galaxies at $z=0$ in the Illustris simulation (considering all the particles
in galaxies). Nevertheless, we hardly reproduce the most recent part of the
SFH from the 2D-SPS analysis, where the number of galaxies used for our
averaged SFH also drops. We note that our 2D-SED fitting method was calibrated
in order to successfully reproduce the first instants of the stellar mass
assembly in massive galaxies, i.e., $t_{\text{5}}$, $t_{\text{10}}$,
$t_{\text{25}}$, along with $t_{\text{mw}}$. Thus, even though we successfully
determine these mass-fraction formation times for the typical SFH of massive
galaxies, our method fails to reproduce this SFH at higher ages of the
Universe. However, the recovery of the whole SFH of massive galaxies and,
consequently, their typical SFH throughout the whole redshift range, is beyond
the scope of this work.

We do not find any starburst epoch in the typical SFH from the 2D-SPS analysis
on the $1<z<4$ sample, and the other SFHs from the simulated particles, albeit
34% (13%) of galaxies in our $1<z<4$ sample reach 50 (100) M${}_{\odot}/$yr at
some time of their 2D-SPS-derived galaxy SFH because of the young population
assumptions in the SED-fits. These starburst episodes are also present in the
ground-truth SFHs of this $1<z<4$ sample with similar numbers: 31% (14%) of
galaxies reach 50 (100) M${}_{\odot}/$yr at some time of their ground-truth
galaxy SFH. Starburst events, which are usually short-lived, $\sim 100$ Myr
(Tacconi et al. 2008, Wuyts et al. 2011, Espino-Briones et al. 2022), have
been proven to occur during the evolution of some galaxies, as it can be
inferred from the high SFR values of Main Sequence outliers or those of the
population of very luminous high-redshift submillimeter galaxies, discovered
by Smail et al. (1997). The reason why these brief intense star formation
episodes are not present in our median SFHs is that they averaged out when we
consider a whole population of galaxies, in the same way that they do not
appear when calculating the cosmic star formation history (Madau & Dickinson,
2014). In addition to this, the Illustris simulation is known to have has a
paucity of strong starburst galaxies, i.e., a fewer fraction of galaxies that
lie significantly above the SFMS when compared to observations (Sparre et al.
2015), which appears to be a consequence of the insufficient resolution of
Illustris to resolve the sub-kiloparsec starbursting regions (Sparre &
Springel 2016).

####  6.1.3 The Variety of the Start of SFHs in Massive Galaxies

Figure 13: Histograms of $t_{5}$ and $t_{1}$ calculated from the SFHs of
galaxies in the different samples considered. Our main sample at $1<z<4$ is
represented by the unfilled histograms: values calculated from the 2D-SPS-
derived galaxy SFHs in red and those calculated from the ground-truth SFHs
built from the database (DB), with a dashed blue line. For the 2D-SPS-derived
values, we only show the median $t_{5}$ or $t_{1}$ for each galaxy (calculated
out of the 300 SFHs per galaxy). The filled histograms represent the
distribution of both mass-fraction formation times calculated from the SFHs
built by considering all particles belonging to galaxies in the database:
descendants at $z=0$ (in orange) and to the whole population of very massive
galaxies at $z=0$ (grey). Hatched histograms show the distributions for these
two samples when only particles inside 2$\times$r${}_{\text{hm}}$ are taken
into account, with the same color code as before. Median and quartiles are
shown as segments on the top.

finally, we briefly concentrate on when the star formation started for the
individual galaxies in the subsamples. We assume the start of the star
formation for each galaxy can be given by its $t_{5}$ (calculated from their
individual galaxy SFH) or, even, by its $t_{1}$.

Fig. 13 shows the diversity in the values of $t_{5}$ (left) and $t_{1}$
(right) for the three subsamples of galaxies and the different methods used to
calculate their SFHs: the main sample at $1<z<4$ (2D-SPS method and ground-
truth from particles in r${}_{\text{phot}}$), their descendants at $z=0$
(ground-truth from all particles and from particles inside $2\times$rhm), and
the whole population of very massive galaxies at $z=0$ (all particles and
particles inside $2\times$rhm). The median values of both $t_{5}$ and $t_{1}$
for the main $1<z<4$ sample given by the 2D-SPS are consistent with the
ground-truth for this sample (with differences of 0.26 and 0.19 Gyr,
respectively), even though the histogram for our 2D-SPS method predicts that
some galaxies start to form a bit earlier ($\sim 0.2-0.4$ Gyr earlier) than
expected according to the database. Regarding the distributions of the two
subsamples at $z=0$, their median $t_{5}$ and $t_{1}$ are similar and also
consistent with each other.

Fig. 13 can be useful if we want to know where to look for the first star
formation episodes in massive galaxies. According to the ground-truth
distribution of $t_{5}$ and $t_{1}$ for the $1<z<4$ sample of massive
galaxies, 25% of these galaxies would have formed 1% (5%) of its stellar mass
by $z\sim 8$ (6), and 50% of the galaxies by $z\sim 7$ (5).

In Table 4, we summarize the mass-fraction formation times (and corresponding
redshifts) for all the SFHs shown in figs. 11 and 12. We have calculated an
estimation of the uncertainties for the 2D-SPS-derived values by taking into
account the (300) individual SFHs for each galaxy and by calculating the mass-
fraction formation times for all of them. We estimate the uncertainty for the
2D-SPS-derived formation times as the quartile values of these distributions.

##  7 SUMMARY AND CONCLUSIONS

We assess the power of broad-band HST+JWST imaging data for the analysis of
the earliest evolutionary phases of galaxies that evolve into local massive
objects (M${}_{\ast}>10^{11}$ M⊙) by analysing the SFHs of progenitors at high
redshift in deep cosmological surveys. For that purpose, we apply stellar
population synthesis in 2D to a sample of 221 Illustris-1 simulated
M${}_{\ast}>10^{10}$ M⊙ galaxies at $1<z<4$ that evolve into
M${}_{\ast}>10^{11}$ M⊙ galaxies at $z=0$. We use ACS, WFC3, and NIRCam data
in the optical and near-infrared from the synthetic Illustris-1 deep survey
images (S+17). We measure SEDs for galaxies in the sample from both integrated
and 2D multi-wavelength photometry on these images, previously processed to
mimic the depths of CANDELS imaging data from HST and ongoing CEERS
observations with JWST along with the source extraction and analysis
procedures to be applied to data from these surveys. From the SED modeling, we
derive the SFH of each source by combining the 2D information from the fits,
and compare it with the ground-truth SFH given by the simulated particles
belonging to each galaxy in the Illustris-1 database. In this work, we focus
on determining the capabilities of broad-band HST+JWST for determining the
first episodes in the stellar mass assembly, and our main findings are the
following:

(i) We evaluate the success of our 2D-SPS method in recovering the earliest
phases of the stellar mass assembly by comparing the mass-fraction formation
times $t_{5}$, $t_{10}$, and $t_{25}$ calculated for each galaxy from its
2D-SPS-derived galaxy SFH with those given by its ground-truth SFH. We find
our 2D-SPS method successfully recovers these quantities with a median
relative offset between our values and ground-truth of 4.3%, 1.7%, and 0.1%
for $t_{5}$, $t_{10}$, and $t_{25}$, respectively, and scatters of 0.21, 0.24,
and 0.28 dex. Additionally, no systematic effects are observed as a function
of galaxy redshift.

(ii) We build the median SFH of our $1<z<4$ sample from the 2D-SPS-derived
individual galaxy SFHs to infer the mass-fraction formation times of the
sample as a whole. Thus, local M${}_{\ast}>10^{11}$ M⊙ galaxies would have
assembled 5% of their stellar mass present at $z\sim 1$ by $z=4.5$, 10% of
their mass by $z=3.7$, and 25% by $z=2.67$, or equivalently, when the age of
the Universe was 1.38, 1.72, and 2.50 Gyr, respectively. These ages agree with
their ground-truth values derived from the typical SFH built from the
simulated particles of galaxies in the sample, with differences of 0.009,
0.002, and 0.02 dex for $t_{5}$, $t_{10}$, and $t_{25}$, respectively.
Nevertheless, our method fails to reproduce the shape of this SFH for at
$z\lesssim 1$, mainly because of the fewer number of galaxies used to compute
the typical SFH at this redshift range.

(iii) We compare the formation times derived from the 2D-SPS analysis with the
values obtained from the ground-truth median SFH at $z<1$ of local
M${}_{\ast}>10^{11}$ M⊙ galaxies and of the descendants of the high-redshift
sample. The mass-fraction formation epochs from the 2D-SPS are systematically
earlier than those inferred from the descendants of this sample. Besides, the
latter are higher than the ones inferred from the whole population of
M${}_{\ast}>10^{11}$ M⊙ galaxies. In both cases, these shifts in the formation
redshifts, of $\sim 200-400$ Myr, are due to the progenitor bias. This bias
arises as a consequence of the mass-cutoff imposed on the high-redshift
progenitor sample, aimed at reproducing the selection suffered in actual
HST+JWST survey observations, and which makes this sample to contain only the
most massive (and, thus, older) progenitors of local M${}_{\ast}>10^{11}$ M⊙
galaxies. Given the way our progenitor sample has been selected, their
descendants at $z=0$ also present a median stellar mass that is slightly
higher than the whole population of local M${}_{\ast}>10^{11}$ M⊙ galaxies,
which leads to mass-fraction formation times closer to the Big Bang.

(iv) Regarding the variety of the $t_{1}$, and $t_{5}$ values, the
distribution of these formation times shows that 25% of our $1<z<4$ sample of
progenitors formed 1% (5%) of their stellar mass present at the redshift of
observation by $z\sim 8$ (6), and 50% of the galaxies by $z\sim 7$ (5).

The results from this work should be taken into consideration when applying
2D-SPS analysis methods to massive high-redshift galaxies observed with
HST+JWST in broad-band imaging, even if the numerical values of the mass-
fraction formation times inferred from real data could vary with respect to
those of the simulated data, on average or as a function of galaxy properties,
such as mass, redshift or even environment.

##  8 Acknowledgments

This research has been funded by grants PGC2018-093499-B-I00 and
RTI2018-096188-B-I00 funded by MCIN/AEI/10.13039/501100011033. ÁGA
acknowledges the support of the Universidad Complutense de Madrid through the
predoctoral grant CT17/17-CT18/17. AY is supported by an appointment to the
NASA Postdoctoral Program (NPP) at NASA Goddard Space Flight Center,
administered by Oak Ridge Associated Universities under contract with NASA. LC
acknowledges financial support from Comunidad de Madrid under Atracción de
Talento grant 2018-T2/TIC-11612. RMM acknowledge support from Spanish
Ministerio de Ciencia, Innovación y Universidades through grant
PGC2018-093499-B-I00, MDM-2017-0737 Unidad de Excelencia “Maria de
Maeztu”-Centro de Astrobiología (INTA-CSIC) by the Spanish Ministry of Science
and Innovation/State Agency of Research MCIN/AEI/ 10.13039/501100011033 and by
“ERDF A way of making Europe”, and also from the Instituto Nacional de Técnica
Aeroespacial SHARDS-JWST project through the PRE-SHARDSJWST/2020 grant.

## References

  * Abdurro’uf et al. (2022) Abdurro’uf, Lin, Y.-T., Hirashita, H., et al. 2022, ApJ, 926, 81. doi:10.3847/1538-4357/ac439a
  * Akhshik et al. (2022) Akhshik, M., Whitaker, K. E., Leja, J., et al. 2022, arXiv:2203.04979
  * Alcalde Pampliega et al. (2019) Alcalde Pampliega, B., Pérez-González, P. G., Barro, G., et al. 2019, ApJ, 876, 135. doi:10.3847/1538-4357/ab14f2
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33. doi:10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123. doi:10.3847/1538-3881/aabc4f
  * Beckwith et al. (2006) Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729. doi:10.1086/507302
  * Bellardini et al. (2022) Bellardini, M. A., Wetzel, A., Loebman, S. R., et al. 2022, arXiv:2203.03653
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393. doi:10.1051/aas:1996164
  * Bouwens et al. (2004) Bouwens, R. J., Illingworth, G. D., Blakeslee, J. P., et al. 2004, ApJ, 611, L1. doi:10.1086/423786
  * Bradley et al. (2019) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2019, Zenodo
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000. doi:10.1046/j.1365-8711.2003.06897.x
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682. doi:10.1086/308692
  * Chamorro-Cazorla et al. (2022) Chamorro-Cazorla, M., Gil de Paz, A., Castillo-Morales, A., et al. 2022, A&A, 657, A95. doi:10.1051/0004-6361/202141930
  * Chabrier (2003) Chabrier, G. 2003, ApJ, 586, L133. doi:10.1086/374879
  * Cook et al. (2016) Cook, B. A., Conroy, C., Pillepich, A., et al. 2016, ApJ, 833, 158. doi:10.3847/1538-4357/833/2/158
  * Costantin et al. (2021) Costantin, L., Pérez-González, P. G., Méndez-Abreu, J., et al. 2021, ApJ, 913, 125. doi:10.3847/1538-4357/abef72
  * Costantin et al. (2022) Costantin, L., Pérez-González, P. G., Méndez-Abreu, J., et al. 2022, ApJ, 929, 121. doi:10.3847/1538-4357/ac5a57
  * Cowie et al. (1996) Cowie, L. L., Songaila, A., Hu, E. M., et al. 1996, AJ, 112, 839. doi:10.1086/118058
  * Dekel & Burkert (2014) Dekel, A. & Burkert, A. 2014, MNRAS, 438, 1870. doi:10.1093/mnras/stt2331
  * Diemer (2018) Diemer, B. 2018, ApJS, 239, 35. doi:10.3847/1538-4365/aaee8c
  * Di Matteo et al. (2013) Di Matteo, P., Haywood, M., Combes, F., et al. 2013, A&A, 553, A102. doi:10.1051/0004-6361/201220539
  * Domínguez Sánchez et al. (2016) Domínguez Sánchez, H., Pérez-González, P. G., Esquej, P., et al. 2016, MNRAS, 457, 3743. doi:10.1093/mnras/stw201
  * Dunlop et al. (2021) Dunlop, J. S., Abraham, R. G., Ashby, M. L. N., et al. 2021, JWST Proposal. Cycle 1, 1837
  * Elias et al. (2018) Elias, L. M., Sales, L. V., Creasey, P., et al. 2018, MNRAS, 479, 4004. doi:10.1093/mnras/sty1718
  * Ellis et al. (2013) Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, ApJ, 763, L7. doi:10.1088/2041-8205/763/1/L7
  * Erb et al. (2006) Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813. doi:10.1086/503623
  * Espino-Briones et al. (2022) Espino-Briones, N., Pérez-González, P. G., Zamorano, J., et al. 2022, MNRAS, 513, 1175. doi:10.1093/mnras/stac728
  * finkelstein et al. (2017) finkelstein, S. L., Dickinson, M., Ferguson, H. C., et al. 2017, JWST Proposal ID 1345. Cycle 0 Early Release Science, 1345
  * Forrest et al. (2020a) Forrest, B., Annunziatella, M., Wilson, G., et al. 2020, ApJ, 890, L1. doi:10.3847/2041-8213/ab5b9f
  * Forrest et al. (2020b) Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020, ApJ, 903, 47. doi:10.3847/1538-4357/abb819
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485. doi:10.1007/s11214-006-8315-7
  * Gil de Paz & Madore (2002) Gil de Paz, A. & Madore, B. F. 2002, AJ, 123, 1864. doi:10.1086/339480
  * Glazebrook et al. (2017) Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Nature, 544, 71. doi:10.1038/nature21680
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35. doi:10.1088/0067-0049/197/2/35
  * Ho et al. (2018) Ho, I.-T., Meidt, S. E., Kudritzki, R.-P., et al. 2018, A&A, 618, A64. doi:10.1051/0004-6361/201833262
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90. doi:10.1109/MCSE.2007.55
  * Illingworth et al. (2013) Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6. doi:10.1088/0067-0049/209/1/6
  * Iyer et al. (2020) Iyer, K. G., Tacchella, S., Genel, S., et al. 2020, MNRAS, 498, 430. doi:10.1093/mnras/staa2150
  * Ji & Giavalisco (2022) Ji, Z. & Giavalisco, M. 2022, arXiv:2204.02414
  * Johnston et al. (2022) Johnston, E. J., Häußler, B., & Jegatheesan, K. 2022, MNRAS. doi:10.1093/mnras/stac1725
  * Jonsson (2006) Jonsson, P. 2006, MNRAS, 372, 2. doi:10.1111/j.1365-2966.2006.10884.x
  * Jonsson et al. (2010) Jonsson, P., Groves, B. A., & Cox, T. J. 2010, MNRAS, 403, 17. doi:10.1111/j.1365-2966.2009.16087.x
  * Kitzbichler & White (2007) Kitzbichler, M. G. & White, S. D. M. 2007, MNRAS, 376, 2. doi:10.1111/j.1365-2966.2007.11458.x
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36. doi:10.1088/0067-0049/197/2/36
  * Koekemoer et al. (2013) Koekemoer, A. M., Ellis, R. S., McLure, R. J., et al. 2013, ApJS, 209, 3. doi:10.1088/0067-0049/209/1/3
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231. doi:10.1046/j.1365-8711.2001.04022.x
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3. doi:10.1086/313233
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, ARA&A, 52, 415. doi:10.1146/annurev-astro-081811-125615
  * Maiolino et al. (2008) Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&A, 488, 463. doi:10.1051/0004-6361:200809678
  * Mannucci et al. (2009) Mannucci, F., Cresci, G., Maiolino, R., et al. 2009, MNRAS, 398, 1915. doi:10.1111/j.1365-2966.2009.15185.x
  * Marsan et al. (2022) Marsan, Z. C., Muzzin, A., Marchesini, D., et al. 2022, ApJ, 924, 25. doi:10.3847/1538-4357/ac312a
  * Méndez-Abreu et al. (2021) Méndez-Abreu, J., de Lorenzo-Cáceres, A., & Sánchez, S. F. 2021, MNRAS, 504, 3058. doi:10.1093/mnras/stab1064
  * Nelson et al. (2015) Nelson, D., Pillepich, A., Genel, S., et al. 2015, Astronomy and Computing, 13, 12. doi:10.1016/j.ascom.2015.09.003
  * Nelson et al. (2016) Nelson, E. J., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2016, ApJ, 828, 27. doi:10.3847/0004-637X/828/1/27
  * Norgaard-Nielsen & Perez-Gonzalez (2017) Norgaard-Nielsen, H. U. & Perez-Gonzalez, P. G. 2017, JWST Proposal. Cycle 1, 1283
  * Oesch et al. (2010a) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2010, ApJ, 709, L16. doi:10.1088/2041-8205/709/1/L16
  * Oesch et al. (2010b) Oesch, P. A., Bouwens, R. J., Carollo, C. M., et al. 2010, ApJ, 709, L21. doi:10.1088/2041-8205/709/1/L21
  * Oke & Gunn (1983) Oke, J. B. & Gunn, J. E. 1983, ApJ, 266, 713. doi:10.1086/160817
  * Ono et al. (2013) Ono, Y., Ouchi, M., Curtis-Lake, E., et al. 2013, ApJ, 777, 155. doi:10.1088/0004-637X/777/2/155
  * Pérez-González et al. (2003) Pérez-González, P. G., Gil de Paz, A., Zamorano, J., et al. 2003, MNRAS, 338, 508. doi:10.1046/j.1365-8711.2003.06077.x
  * Pérez-González et al. (2008) Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234. doi:10.1086/523690
  * Pontoppidan et al. (2016) Pontoppidan, K. M., Pickering, T. E., Laidler, V. G., et al. 2016, Proc. SPIE, 9910, 991016. doi:10.1117/12.2231768
  * Reback et al. (2020) Reback, J., McKinney, W., Jbrockmendel, et al. 2020, Zenodo
  * Renzini (2006) Renzini, A. 2006, ARA&A, 44, 141. doi:10.1146/annurev.astro.44.051905.092450
  * Rieke et al. (2019) Rieke, M., Arribas, S., Bunker, A., et al. 2019, BAAS, 51, 45
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, MNRAS, 449, 49. doi:10.1093/mnras/stv264
  * Sales et al. (2015) Sales, L. V., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 447, L6. doi:10.1093/mnrasl/slu173
  * Sánchez-Menguiano et al. (2020) Sánchez-Menguiano, L., Sánchez, S. F., Pérez, I., et al. 2020, MNRAS, 492, 4149. doi:10.1093/mnras/staa088
  * Skelton et al. (2014) Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24. doi:10.1088/0067-0049/214/2/24
  * Smail et al. (1997) Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJ, 490, L5. doi:10.1086/311017
  * Snyder et al. (2015) Snyder, G. F., Torrey, P., Lotz, J. M., et al. 2015, MNRAS, 454, 1886. doi:10.1093/mnras/stv2078
  * Snyder et al. (2017) Snyder, G. F., Lotz, J. M., Rodriguez-Gomez, V., et al. 2017, MNRAS, 468, 207. doi:10.1093/mnras/stx487
  * Sparre et al. (2015) Sparre, M., Hayward, C. C., Springel, V., et al. 2015, MNRAS, 447, 3548. doi:10.1093/mnras/stu2713
  * Sparre & Springel (2016) Sparre, M. & Springel, V. 2016, MNRAS, 462, 2418. doi:10.1093/mnras/stw1793
  * Tacchella et al. (2015) Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Science, 348, 314. doi:10.1126/science.1261094
  * Tacchella et al. (2016) Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 458, 242. doi:10.1093/mnras/stw303
  * Tacchella et al. (2018) Tacchella, S., Carollo, C. M., Förster Schreiber, N. M., et al. 2018, ApJ, 859, 56. doi:10.3847/1538-4357/aabf8b
  * Tacchella et al. (2019) Tacchella, S., Diemer, B., Hernquist, L., et al. 2019, MNRAS, 487, 5416. doi:10.1093/mnras/stz1657
  * Tacconi et al. (2008) Tacconi, L. J., Genzel, R., Smail, I., et al. 2008, ApJ, 680, 246. doi:10.1086/587168
  * Tody (1993) Tody, D. 1993, Astronomical Data Analysis Software and Systems II, 52, 173
  * Torrey et al. (2015) Torrey, P., Snyder, G. F., Vogelsberger, M., et al. 2015, MNRAS, 447, 2753. doi:10.1093/mnras/stu2592
  * Valentino et al. (2020) Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93. doi:10.3847/1538-4357/ab64dc
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22. doi:10.1109/MCSE.2011.37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261. doi:10.1038/s41592-019-0686-2
  * Vogelsberger et al. (2014a) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Nature, 509, 177. doi:10.1038/nature13316
  * Vogelsberger et al. (2014b) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518. doi:10.1093/mnras/stu1536
  * Wang et al. (2017) Wang, W., Faber, S. M., Liu, F. S., et al. 2017, MNRAS, 469, 4063. doi:10.1093/mnras/stx1148
  * Williams et al. (2018) Williams, C. C., Curtis-Lake, E., Hainline, K. N., et al. 2018, ApJS, 236, 33. doi:10.3847/1538-4365/aabcbb
  * Wuyts et al. (2011) Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, ApJ, 742, 96. doi:10.1088/0004-637X/742/2/96
  * Zahid et al. (2011) Zahid, H. J., Kewley, L. J., & Bresolin, F. 2011, ApJ, 730, 137. doi:10.1088/0004-637X/730/2/137
  * Zolotov et al. (2015) Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327. doi:10.1093/mnras/stv740