# Revealing Galaxy Candidates out to $z\sim 16$ with JWST Observations of the
Lensing Cluster SMACS0723

Hakim Atek1, Marko Shuntov1, Lukas J. Furtak2, Johan Richard3, Jean-Paul
Kneib4, Guillaume Mahler5, Adi Zitrin2, H.J. McCracken1, Stéphane Charlot1,
Jacopo Chevallard6, Iryna Chemerynska1
1Institut d’Astrophysique de Paris, CNRS, Sorbonne Université, 98bis Boulevard
Arago, 75014, Paris, France
2Physics Department, Ben-Gurion University of the Negev, P.O. Box 653, Be’er-
Sheva 84105, Israel
3Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230, Saint-Genis Laval, France
4Laboratoire d’Astrophysique, Ecole Polytechnique Fédérale de Lausanne,
Observatoire de Sauverny, CH-1290 Versoix, Switzerland
5Institute for Computational Cosmology, Durham University, South Road, Durham
DH1 3LE, UK
6Department of Physics, University of Oxford, Denys Wilkinson Building, Keble
Road, Oxford OX1 3RH, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

One of the main goals of the JWST is to study the first galaxies in the
Universe. We present a systematic photometric analysis of very distant
galaxies in the first JWST deep field towards the massive lensing cluster
SMACS0723. As a result, we report the discovery of two galaxy candidates at
$z\sim 16$, only $250$ million years after the Big Bang. We also identify two
candidates at $z\sim 12$ and 6 candidates at $z\sim 9-11$. Our search extended
out to $z\lesssim 21$ by combining color information across seven NIRCam and
NIRISS filters. By modelling the Spectral Energy Distributions (SEDs) with
EAZY and BEAGLE, we test the robustness of the photometric redshift estimates.
While their intrinsic (un-lensed) luminosity is typical of the characteristic
luminosity L∗ at $z>10$, our high-redshift galaxies typically show small sizes
and their morphologies are consistent with disks in some cases. The highest-
redshift candidates have extremely blue UV-continuum slopes $-3<\beta<-2.4$,
young ages $\sim 10-100$ Myr, and stellar masses around
$\log(M_{\star}/\mathrm{M}_{\odot})=8.8$ inferred from their SED modeling,
which indicate a rapid build-up of their stellar mass. Our search clearly
demonstrates the capabilities of JWST to uncover robust photometric candidates
up to very high redshifts, and peer into the formation epoch of the first
galaxies.

###### keywords:

galaxies: high-redshift – dark ages, reionization, first stars – galaxies:
dwarfs – galaxies: evolution – gravitational lensing: strong – cosmology:
observations

††pubyear: 2022††pagerange: Revealing Galaxy Candidates out to $z\sim 16$ with
JWST Observations of the Lensing Cluster SMACS0723–LABEL:lastpage

##  1 Introduction

Understanding the emergence of the first structures that hosted star formation
has been one of the most important goals in modern astrophysics. One of the
main routes to address this question is to identify the first generation of
galaxies that we believe formed at $z\gtrsim 15$. Great progress has been made
in the last decade to uncover the most distant galaxies out to $z\sim 11$,
essentially thanks to the Hubble Space Telescope (HST) and 8-10 m class
ground-based telescopes. The record-holder for the farthest galaxy known with
a spectroscopic confirmation is GN-z11 (Oesch et al., 2016; Jiang et al.,
2021), at a redshift of $z=10.96$. Remarkably, GN-z11 is a relatively bright
and massive galaxy ($M_{\mathrm{UV}}=-22.1$ AB), which challenges current
galaxy formation models. The existence of a galaxy with a billion solar masses
only $\sim 400$ Myr after the Big Bang indicates that galaxy formation was
well underway very early in the history of the Universe.

At slightly lower redshifts $z\sim 7-9$, large samples of nearly 2,000
galaxies have been photometrically identified, thanks mostly to the infrared
sensitivity of the Wide Field Camera Three (WFC3) aboard HST (e.g. Finkelstein
et al., 2015; Bouwens et al., 2021) or ground-based observations (e.g.,
Kauffmann et al., 2022). These large samples resulted in strong constraints on
the shape of the ultra-violet (UV) luminosity function of galaxies at those
epochs and its evolution with redshift. Moreover, Spitzer Space Telescope
observations allowed to probe the rest-frame optical emission of these high-
redshift galaxies, which enabled us to constrain their stellar mass (e.g. Song
et al., 2016; Bhatawdekar et al., 2019; Kikuchihara et al., 2020; Stefanon et
al., 2021), albeit with important uncertainties due to source confusion and
nebular emission contamination (Grazian et al., 2015; Furtak et al., 2021).

More recently, using deep Spitzer and ground-based observations, Harikane et
al. (2022) reported the detection of two galaxy candidates at $z\sim 12-13$.
Atacama Large Millimeter/submillimeter Array (ALMA) follow-up observations
have revealed the tentative ($2\sigma$) detection of an [O iii]$\lambda
88\mu\mathrm{m}$ emission line, putting HD1 at a redshift of $z=13.27$. Again,
these two candidates do not align with theoretical models, which predict a
significantly lower density of such bright galaxies at $z>10$. While wide
extra-galactic surveys have unveiled these surprisingly large numbers of
luminous galaxies at $z\sim 9-11$, deeper observations at longer wavelengths
are needed to, (i) push the redshift frontier to uncover farther away galaxies
and (ii) constrain the number density of faint galaxies at early times.

The advent of the JWST marks a new era in the detection and, more importantly,
the characterization of early star-forming galaxies. The most important
improvement is the significantly higher spatial resolution of JWST compared to
Spitzer, which enables accurate photometric measurements in the near and mid-
infrared range corresponding to the rest-frame optical for galaxies at $z>6$.
This is crucial for deriving physical properties of galaxies through spectral
energy distribution (SED) fitting, since it allows us to account for a larger
dynamic range in the age of the underlying stellar populations. Also, for the
first time, we will be able to obtain rest-frame optical spectra for these
early-star forming galaxies, giving us access to several optical emission
line, which are the gold standards for measuring galaxy properties such as
star formation rate (SFR), gas-phase metallicity, dust content, etc.

The first JWST observations consist of the Early Release Observations (ERO),
which were obtained to showcase the observatory capabilities, and the Early
Release Science (ERS) programs, aimed at testing the variety of instruments
observing modes and to help the community understand and exploit JWST data.
One of the primary targets of the ERO observations is the galaxy cluster SMACS
J0723.3-7327 (SMACS0723 hereafter; Ebeling et al., 2001; Repp & Ebeling, 2018)
which was recently observed with HST as part of the Reionization Lensing
Cluster Survey program (RELICS; Coe et al., 2019; Salmon et al., 2020). One of
the accepted ERS programs (1324, PI: Treu) also targets a lensing cluster,
Abell 2744 (A2744), which is one of the Hubble Frontier Fields (HFF; Lotz et
al., 2017) clusters. These observations follow in the steps of numerous HST
observing campaigns and efforts to use massive galaxy clusters as
gravitational telescopes. The magnification provided by the strong
gravitational lensing (SL) enables us to reach intrinsically fainter and
lower-mass galaxies and thus to constrain the very faint and low-mass ends of
the luminosity and stellar mass function (e.g. Atek et al., 2014a; Atek et
al., 2018; Livermore et al., 2017; Bouwens et al., 2017; Ishigaki et al.,
2018; Bhatawdekar et al., 2019; Kikuchihara et al., 2020; Furtak et al.,
2021).

In this work, we use the JWST ERO data of SMACS0723 to identify and
characterize the most distant galaxies in the Universe. We use imaging data in
7 broad-band filters to select $z\gtrsim 9$ galaxies based on the Lyman break
technique coupled with SED fitting to estimate their photometric redshifts and
physical properties. This paper is structured as follows. In section 2, we
present the observational data. In section 3, we describe how the galaxy
sample was selected. The SL model used in this work and our SED-fitting
analysis are shown in sections 4 and 5 respectively. Finally, we compare our
estimate with a compilation of literature results and discuss the implications
on galaxy formation models in section 6 before presenting our summary and
conclusions in section 7. Throughout the paper, magnitudes are in the AB
system (Oke & Gunn, 1983) and we adopt a cosmology with H${{}_{0}}=70$ km s-1
Mpc-1, $\Omega_{\Lambda}=0.7$, and $\Omega_{m}=0.3$.

##  2 Observations

Table 1: Limiting AB magnitudes (5$\sigma$) of the JWST/NIRCam and JWST/NIRISS imaging data used in this work. The depths were computed using random 0.3″circular apertures. SW filter | F090W | F115W | F150W | F200W
---|---|---|---|---
Depth | 28.0 | 27.8 | 28.8 | 28.9
LW filter | F277W | F356W | F444W |
Depth | 29.2 | 29.3 | 29.1 |
Figure 1: JWST/NIRCam and JWST/NIRISS through-put curves of the filter set of
the JWST ERO observations of SMACS0723 used in this work.

The first JWST science observations, targeting the lensing cluster SMACS0723,
were obtained as part of the Early Release Observations (ERO) program (ID
2736; PI: Pontoppidan), and have been released on Wednesday, July 13th. They
consist of deep multi-wavelength imaging with NIRCam (Near-Infrared Camera)
and MIRI (Mid-Infrared Instrument), NIRSpec (Near-Infrared Spectrograpgh)
multi-object Spectroscopy, and NIRISS (Near-Infrared Imager and Siltless
Spectrograpgh) wide-field slitless spectroscopy. We retrieved the data from
the Mikulski Archive for Space Telescopes (MAST). Both the short wavelength
(SW) and the long wavelength (LW) channels of NIRCam were used simultaneously
to obtain uniform spectral coverage from $\lambda\sim 1$ µm to $\lambda\sim 5$
µm. Tab. 1 summarizes the set of filters used in this work, and Fig. 1 shows
their transmission curves. Note that we complement the 6 available NIRCam
filters with F115W-band imaging data from the Near-Infrared Imager and
Slitless Spectrograph (NIRISS; Doyon et al., 2012). These however have a
smaller field-of-view and only cover the core of the cluster.

For each of the filters, we reduced the raw uncal files using the most recent
JWST calwebb pipeline v1.6.2 and the calibration context CRDS_CTX =
jwst_0970.pmap. In particular, the latest reference files contain updates from
in-flight calibration, whereas ground calibrations were used in the early
reductions of the ERO data. First, we applied the Detector1 stage to all
files. Second, we ran the images through the second reduction stage of the
pipeline Image2, which performs flat-field correction, WCS registration, and
photometric calibration. Due to significant levels of background residuals in
the calibrated images, and following recommendations from the ERS CEERS
program, we independently ran background subtraction on individual cal files.
In addition, the first NIRCam calibration observations have shown that in most
cases the detectors/filters were more sensitive than predicted (Rigby et al.,
2022). Therefore, we updated the photometric zeropoints of the cal files with
recent measurements of Boyer et al. (2022) and G.
Brammer111https://github.com/gbrammer/grizli/pull/107. Finally, cal files were
processed through the Image3 stage of the pipeline to create the final
drizzled mosaic in each filter. All images were registered to GAIA DR2
astrometry and the relative registration was refined using TweakReg. We
adopted a pixel scale of 0.03″ for both SW and LW channels.

In order to increase the detection efficiency in the extraction procedure
described in section 3.1, we created a deep LW image by stacking the three LW
filters using their respective weight maps. We also create a stacked SW image
of the non-detection filters blue-ward of the break.

One common limitation for the detection of faint sources behind massive
cluster fields is the intra-cluster light (ICL). This emission in the central
region of the cluster is emitted by tidally stripped stars and bright cluster
members (e.g. Atek et al., 2014a; Montes, 2022). In order to reduce the impact
of this diffuse light, we apply a median filter to the detection image. We
adopted a filter size of $\sim 2$″ $\times$ 2″, which significantly reduces
the diffuse ICL and BCG light while minimizing residual artifacts. Note that
the photometry measurement itself is still performed in the original image of
each filter.

Finally, we also use spectroscopic observations of SMACS0723 obtained with the
Near-Infrared Spectrograph (NIRSpec; Jakobsen et al., 2022) as part of the
same ERO program. We obtained level 3 data products from MAST which include 2D
and 1D spectra of selected targets in the cluster field.

##  3 Sample

###  3.1 Photometric Catalogs

Figure 2: Magnitude number counts from the SE++ extraction in the seven NIRCam
and NIRISS filters. Also over-plotted are the HST WFC3/IR F140W number counts
from RELICS (Coe et al., 2019; Salmon et al., 2020). The plot illustrates the
improvement in depth and the increase of the number of faint sources detected
with JWST.

We carry out source detection and photometry measurements using
SourceXtractor++ v0.17 (Bertin et al., 2020; Kümmel et al., 2020).
SourceXtractor++ (SE++) is the successor of the classic and widely used
SExtractor2, rewritten in C++ in order to optimize computational efficiency
and parallelization across multiple cores. The novelty and advantage of the
code is its support for two-level detection thresholding, flexible multi-
object and multi-frame model fitting as well as allowing for custom definition
of models. Additionally, it can operate on multiple images based on their WCS
information, bypassing the need to re-sample images on the same pixel scale.
SE++ has been extensively tested in its performance on photometry and
morphological parameter recovery from model fitting in the Euclid Morphology
Challenge, where it achieved the highest scores among all the tested codes
(Euclid Collaboration et al., 2022b, a).

For the detection of sources, we use the ICL-subtracted LW stack (cf. section
3.1) and its corresponding weight image as a weight map. The background map is
computed on a grid with cell size of 128 pixels and smoothed with a box-filter
of 5 pixels per side. Sources are detected via two-level thresholding which
ensures the detection of faint and low surface brightness sources while
minimizing false positives around bright sources. First, we detect all sources
with at least 6 contiguous pixels and an integrated SNR threshold of 1.75. A
second, higher threshold of 2.55 and 18 contiguous pixels is applied, whose
main impact is to clear false positives around bright sources. Finally, a
minimum area of 15 pixels is applied to clean spurious detections near bright
sources and random noise peaks. Detected sources that are close enough such
that they have connected pixels above the thresholds are grouped together and
later fitted simultaneously.

Measurements are then carried out on all 6 NIRCam bands simultaneously.
Additionally, we perform measurements on the NIRISS F115W-band separately
because of its smaller size and shallower depth (cf. Tab. 1), but using the
same detection image. Errors in the photometry are estimated from an RMS map,
obtained as the inverse of the square root of the weight image. The weight
image contains the contribution of the sky, dark, flat and read noise to each
pixel. As such, it considers the Poisson noise of the background, but not of
the sources. Regardless, given the fact that of interest for this paper are
the faintest sources, which are typically dominated by read and background
noise, this approach is sufficient. The photometry for each source is obtained
using AUTO apertures, as well as from model-fitting. The models and the priors
on their parameters can be user-defined, and in our case we use a Sérsic model
(Sérsic, 1963) with the radius, Sérsic index, angle, and axis ratio as free
parameters. We find that this model performs well for the faint and small
sources that we are interested in. The models are convolved with the point-
spread-function (PSF) obtained using WebbPSF v1.0.0. The PSF models are
generated on a pixel scale of $0.03$ ″/pix, on a stamp with a field-of-view of
$\sim 10$″ (301 pixels).

For the purpose of this paper, we use the model magnitudes from the
simultaneous fit of the NIRCam bands and the AUTO magnitudes for the separate
NIRISS measurement. The resulting magnitude number counts are shown in Fig. 2
for all bands used in this work. In order to check the photometric calibration
of these early JWST observations, we conducted a comparison with well-
calibrated HST images of the same cluster from RELICS. We measured fluxes in
the HST/WFC3 F140W and in the JWST/NIRCam F150W images, using the same
apertures and the HST mosaic as a detection image. The results are discussed
in appendix A and show a good agreement with a small residual offset. Fig. 2
also shows consistent number counts and magnitude depths across the JWST
filters. For comparison, we also show the number counts from the RELICS
F140W-band.

###  3.2 Dropout selection

Figure 3: Color-color selection windows using the NIRCam broadband filters for
high-redshift dropouts. Left: for $9<z<11$ Middle: for $11<z<15$; and Right:
for $16<z<21$. The solid lines, from blue to turquoise, represent the color
tracks of star-forming galaxy templates at high redshift generated with BEAGLE
(cf. text for details). The dashed lines, from yellow to red, are quiescent
galaxies, illustrating low-redshift galaxy contaminants, generated with GRASIL
(Silva et al., 1998). For each star-forming template we applied three values
of dust attenuation following the SMC dust law, $A_{V}$=[0,0.25,0.5], which
are illustrated with the brighter blue or red colors in increasing order to
illustrate the effect of dust attenuation on the color tracks. For the
quiescent galaxies, we apply attenuation values in the range $A_{V}$=[1,2,3].
The purple circles represent the colors of a library of M-class stars and
brown dwarfs (Chabrier et al., 2000; Allard et al., 2001), which also
represent potential contaminants. The high-redshift dropouts selected in this
work are represented as green diamonds.

Figure 4: Examples of NIRSpec spectra using the G395M grism for 2 dropout
galaxies at $z\sim 7-9$, displaying very strong H$\beta$ and [O iii] emission
lines.

We proceed to identify high-redshift galaxies using the Lyman-break, or
dropout, selection technique (e.g. Steidel et al., 1996; Atek et al., 2015;
Kawamata et al., 2015; Bouwens et al., 2015). In order to determine the best
color-color selection criteria, we ran a set of simulations using different
galaxy templates. We generated starburst galaxy templates using BEAGLE
(Chevallard & Charlot, 2016, cf. section 5) at increasing redshifts from $z=6$
to $z=28$, including the intergalactic medium (IGM) attenuation from Inoue et
al. (2014). In addition, we explore a range of dust attenuation from $A_{V}=0$
to $A_{V}=0.5$ using the SMC extinction law (Pei, 1992), which has been found
to well match the dust attenuation in high-redshift galaxies (e.g. Capak et
al., 2015; Reddy et al., 2015, 2018a). We then computed synthetic photometry
in the NIRCam broad-band filters and plotted the color tracks of these mock
galaxies in Fig. 3. In order to mitigate contamination from low-redshift
interlopers that mimic the IGM absorption by their intrinsically red colors,
we also explore the position of low-redshift quiescent galaxies, generated
with GRASIL (Silva et al., 1998) and available in the framework of the SWIRES
template library (Polletta et al., 2007), on this color-color diagram. These
again include several values of dust attenuation between $A_{V}=1$ and
$A_{V}=3$. Finally, we include the last significant source of contamination,
cold red stars and brown dwarfs, which can also mimic the Lyman break. We used
the stellar templates from Chabrier et al. (2000) and Allard et al. (2001).
With this information in hand, we define the best color selection window that
maximizes identification of the high-redshift sources while minimizing
contamination for each redshift range. The resulting selection criteria are
illustrated in Fig. 3.

According to these simulations, we select galaxy candidates in the redshift
ranges $z\sim 9-11$ and $z\sim 11-15$ respectively based on the following
criteria:

$\begin{array}[]{l}M_{115}-M_{150}>0.5\\\
M_{115}-M_{150}>1.0+1.4(M_{150}-M_{200})\\\ M_{150}-M_{200}<0.5\end{array}$

and

$\begin{array}[]{l}M_{150}-M_{200}>0.5\\\
M_{150}-M_{200}>1.6+1.8(M_{200}-M_{277})\\\ M_{200}-M_{277}<0.2\end{array}$

In addition to these color criteria, galaxy candidates must also satisfy a
detection significance above 5$\sigma$ in the deep IR stacked image, and
4$\sigma$ in the individual detection filters. We also reject any candidate
that shows a significant detection (at 2$\sigma$ level) in the stacked image
that combines all filters blue-wards of the Lyman break. If a source is not
detected in the dropout filter, we adopt the 1$\sigma$ lower limit measured in
the same aperture to calculate the Lyman break.

For $z\sim 16-21$ candidates, we adopt the following criteria:

$\begin{array}[]{l}M_{200}-M_{277}>0.5\\\
M_{200}-M_{277}>1.2+2.0(M_{277}-M_{356})\\\ M_{277}-M_{356}<0.4\end{array}$

Similarly to the $11<z<15$ selection, all candidates must satisfy the same
detection and non-detection criteria detailed above. All candidates are
visually inspected to check for spurious detection, artifacts, or bad PSF
residuals. We also carefully check the candidates against other potential
sources of contamination. For instance, very strong emission lines in
intermediate-redshift ($z\sim 3-5$) galaxies can contaminate the broadband
flux and create a photometric dropout (e.g. Atek et al., 2011; de Barros et
al., 2014). In particular, strong [O iii]$\lambda\lambda$4959,5007 + H$\beta$
lines are expected in high-redshift galaxies, (e.g. De Barros et al., 2019)
which can lead to red colors. To address this issue, we perform SED fitting
that takes nebular emission lines into account (cf. section 5). Note also that
most of the intermediate-redshift contaminants should be easily detected in
the bluer filters, given the depth of the F090W and F150W filters and the
observed hypothetical excess in the red filters. Finally, we discard all
sources that have point-source morphologies in order to sieve out any stars
that might scatter into our selection window. The colors of our selected
galaxies are shown in Fig. 3 as green diamonds. Note that there are a few
objects selected as high-redshift galaxies in the $z\sim 11-15$ bin which do
not completely satisfy all of the color-criteria. These objects were selected
from their photometric redshift estimates with EAZY (cf. section 5) because
they have unambiguous high-redshift solutions and do not present morphologies
consistent with stars.

In the present study, we have also analyzed the NIRSpec observations to
investigate the spectroscopic properties of the high-redshift galaxies
targeted by this program. Two examples at $z_{\mathrm{spec}}=7.66$ and
$z_{\mathrm{spec}}=8.5$ are shown in Fig. 4. These are the first rest-frame
optical spectra observed for high-redshift ($z>6$) galaxies (cf. also Schaerer
et al., 2022). They show strong H$\beta$ and [Oiii] emission lines. The rest-
frame equivalent width (EW) of these lines reaches up to $\sim 400$Å, which is
still a lower limit because the continuum is hardly detected. This means that
the flux contribution of emission lines to the broadband flux can reach up to
$\sim 30$ % in some galaxies.

Another consequence of the strong lensing is that some background sources will
have multiple images. Using the most recent available parametric lensing model
(cf. section 4), we have not identified any possible multi-image systems. The
magnification factors of the candidates range from $\mu\sim 1$ to $\sim 4$.

##  4 Lensing Model

Figure 5: The cumulative surface area behind SMACS0723 at $z=9$ as a function
of gravitational magnification $\mu$ (cf. section 4), expressed in magnitudes.
The solid line shows the area corresponding to the NIRCam field of view,
whereas the dashed line is for the NIRISS pointing. The un-lensed total survey
area is 7.23 arcmin2, and 2.1 arcmin2 for the NIRCam and NIRISS observations
respectively. For reference, the observed area is $\sim 10.8$ arcmin2 and
$\sim 4.9$ arcmin2 for NIRCam and NIRISS, respectively.

In order to compute the gravitational magnifications of our objects, we use
the most recent strong lensing mass model by Mahler et al. (2022). The model
was built using the lenstool (Kneib et al., 1996; Jullo et al., 2007; Jullo &
Kneib, 2009) software to optimize a parametric SL model for the mass
distribution of the cluster. Following the approach used by Richard et al.
(2014) for the HFF clusters (Lotz et al., 2017) and Fox et al. (2022) for the
RELICS clusters, this model uses a combination of double Pseudo Isothermal
Elliptical potentials (dPIE; Elíasdóttir et al., 2007) for the total mass
distribution at cluster- and galaxy-scales. The parameters of the galaxy-scale
dPIE potentials are matched to the identification of cluster members on the
red sequence from the RELICS HST catalogs (cf. Coe et al. 2019 for details),
fixing the center and elliptical parameters of the potential to its isophotes.
The mass parameters of the dPIE profile (velocity dispersion $\sigma$, core
radius rcore, cut radius rcut) are scaled according to the luminosity L∗ of a
reference galaxy using a constant mass-to-light ratio.

The SL constraints used in Mahler et al. (2022) are based on the released
JWST/NIRCam images (cf. section 2), and make use of spectroscopic redshifts
measured with public data from the Multi Unit Spectroscopic Explorer (MUSE;
Bacon et al., 2010) on ESO’s Very Large Telescope (VLT) by Golubchik et al.
(2022) and with JWST/NIRSpec data in Mahler et al. (2022). In summary, 48
secured images of 16 separate sources were identified, including 5
spectroscopic redshifts. The Mahler et al. (2022) SL model reproduces all
systems with a residual lens plane RMS of $\sim 0.7\arcsec$ between the
predicted and observed locations of all images.

The cumulated surface area of magnification according to the model is shown in
Fig. 5. We make use of this model to compute the magnification factors and
associated errors on each high redshift galaxy candidates. All magnifications
assume $z=15$ as there is only very little variation of magnification with
redshift at $z>7$. lenstool allows us to compute a statistical error on the
magnification by sampling the posterior probability distribution of each
parameter of the model.

##  5 Spectral Energy Distribution Fitting

We continue our selection of high-redshift galaxy candidates by computing
photometric redshifts by fitting galaxy spectral energy distributions (SEDs)
to the multi-wavelength photometry.

The photometric redshifts are obtained with EAZY (Brammer et al., 2008), a
Python-based SED fitting code. EAZY fits a non-negative linear combination of
a set of basis templates to the observed flux densities for each galaxy. The
fit is performed on the flux densities obtained from the model-fitting in SE++
(cf. section 3.1), converted to units of $\mu$Jy, and corrected for Milky Way
extinction internally in the code. Sources with valid flux density
measurements available in at least three bands are considered for fitting. We
use the template set "tweak_fsps_QSF_1_v3.param" that is derived from the
Flexible Stellar Population Synthesis models (FSPS) (Conroy et al., 2009;
Conroy & Gunn, 2010). The allowed redshift range is set to $0.01<z<22$, with a
flat prior on the luminosity. Sources of interest are then selected by
visually inspecting all sources with a best-fit $z_{\mathrm{phot}}>9$. In a
second selection step, candidates are required to show a reasonably good SED
fit.

In addition, we simultaneously infer physical parameters and photometric
redshifts of the selected high-redshift candidates by fitting galaxy SEDs to
the JWST photometry with the BayEsian Analysis of GaLaxy sEds tool (BEAGLE;
Chevallard & Charlot, 2016). This is done to refine our selection and obtain
accurate uncertainty estimates on the physical parameters of the galaxies.
Because the inclusion of nebular emission has been found to be crucial in
fitting broad-band photometry of galaxies (Schaerer & de Barros, 2009, 2010;
Atek et al., 2011; Atek et al., 2014b; McLure et al., 2011; Smit et al., 2014;
Duncan et al., 2014; Reddy et al., 2018b), BEAGLE uses galaxy templates
computed by Gutkin et al. (2016) which combine the stellar population
templates by Bruzual & Charlot (2003) with the photoionization code CLOUDY
(Ferland et al., 2013). It then accounts for IGM attenuation using the Inoue
et al. (2014) models and applies a dust attenuation law (cf. section 5.1 for
details) to the galaxy emission. The BEAGLE tool is fully Bayesian, i.e.
performs a Monte-Carlo Markov Chain (MCMC) analysis of the SED, which makes it
ideally suited to determine both photometric redshifts and several physical
parameters at the same time while robustly quantifying and combining their
uncertainties.

###  5.1 BEAGLE set-up

For our SED-fit with BEAGLE, we assume a constant star-formation history
(SFH). While this is a simplification, it is a common assumption made in high-
redshift galaxy analyses (e.g. Eyles et al., 2007; Stark et al., 2009;
González et al., 2011; Grazian et al., 2015; Kikuchihara et al., 2020) and the
results have been shown to not strongly differ from more flexible analytic SFH
assumptions such as e.g. a delayed exponential SFH (Furtak et al., 2021). In
order to avoid over-fitting due to the relatively low number of photometric
bands available for this work, we limit the fit-parameters to four:

  * •

Photometric redshift $z_{\mathrm{phot}}$ both fixed to the EAZY result first
and then with a uniform prior in the limits of $z_{\mathrm{phot}}\in[0,25]$
for confirming the robustness of the photometric redshift estimate.

  * •

Stellar mass $M_{\star}$ with a log-uniform prior in the limits of
$\log(M_{\star}/\mathrm{M}_{\odot})\in[6,10]$.

  * •

Maximum stellar age $t_{\mathrm{age}}$ with a log-uniform prior in the limits
of $\log(t_{\mathrm{age}}/\mathrm{yr})\in[7,t_{\mathrm{universe}}]$ where
$t_{\mathrm{universe}}$ is the age of the Universe at the redshift of the
galaxy.

  * •

Effective V-band dust attenuation optical depth $\hat{\tau}_{V}$ with a
uniform prior in the limits of $\hat{\tau}_{V}\in[0,0.5]$.

We furthermore fix the stellar metallicity to a constant value of
$Z=0.1\,\mathrm{Z}_{\odot}$ since it has been shown that the broad-band
photometry of high-redshift dropout galaxies is not sensitive to the
metallicity in the very low metallicity regime expected for high-redshift
galaxies (Furtak et al., 2021). Finally, we adopt an SMC dust extinction law
(Pei, 1992) which has been shown to match the observations of high-redshift
galaxies best (Capak et al., 2015; Reddy et al., 2015, 2018a), in particular
in the low-metallicity range that we are probing (Shivaei et al., 2020).

The SED fit is run on all six broad-band filters without taking the
gravitational magnification into account yet, in order to avoid including the
uncertainties of magnification in the SED fit. Instead, we correct the stellar
mass for magnification after the fit following our approach in Furtak et al.
(2021). Note though that the magnifications of our objects are not very high,
$\mu\sim 1-4$.

###  5.2 Results

Table 2: Complete list of our high-redshift candidates identified behind SMACS0723 with their derived parameters. ID | RA | Dec | $z_{\mathrm{phot}}$ | $M_{\mathrm{UV}}$ | $\beta$ | $\log(M_{\star}/\mathrm{M}_{\odot})$ | $\log(t_{\mathrm{age}}/\mathrm{yr})$ | $n$ |  $r_{\rm e}$ [kpc]
---|---|---|---|---|---|---|---|---|---
$z\sim 9-11$ candidates
SMACS_z10a | 7:23:26.252 | -73:26:56.940 | $9.78^{+0.02}_{-0.02}$ | $-18.77\pm 0.20$ |  $-1.72\pm$0.04 | $9.11^{+0.07}_{-0.07}$ | $8.68^{+0.2}_{-0.2}$ | $<1$ | $0.22$
SMACS_z10b | 7:23:22.709 | -73:26:06.183 | $8.88^{+0.02}_{-0.02}$ | $-20.78\pm 0.13$ |  $-1.36\pm$0.19 | $10.20^{+0.03}_{-0.03}$ | $8.74^{+0.2}_{-0.2}$ | $<1$ | $0.65$
SMACS_z10c | 7:23:20.169 | -73:26:04.233 | $9.77^{+0.02}_{-0.02}$ | $-20.19\pm 0.15$ |  $-2.14\pm$0.12 | $9.53^{+0.02}_{-0.02}$ | $8.68^{+0.2}_{-0.2}$ | $<1$ | $0.41$
SMACS_z10d | 7:22:46.696 | -73:28:40.898 | $9.31^{+0.06}_{-0.08}$ | $-19.76\pm 0.18$ |  $-2.22\pm$0.17 | $7.77^{+0.14}_{-0.11}$ | $7.29^{+0.18}_{-0.12}$ | $<1$ | $0.55$
SMACS_z10e | 7:22:45.304 | -73:29:30.557 | $10.89^{+0.16}_{-0.14}$ | $-18.91\pm 0.26$ |  $-2.03\pm$0.19 | $8.51^{+0.22}_{-0.16}$ | $7.49^{+0.26}_{-0.20}$ | $<1$ | $0.33$
SMACS_z11a | 7:22:39.505 | -73:29:40.224 | $11.05^{+0.09}_{-0.08}$ | $-18.55\pm 0.38$ |  $-1.97\pm$0.38 | $8.77^{+0.17}_{-0.24}$ | $8.20^{+0.21}_{-0.32}$ | $<1$ | $0.39$
$z\sim 11-15$ candidates
SMACS_z12a | 7:22:47.380 | -73:30:01.785 | $12.20^{+0.21}_{-0.12}$ | $-19.75\pm 0.23$ |  $-2.69\pm$0.16 | $8.14^{+0.21}_{-0.17}$ | $7.61^{+0.26}_{-0.20}$ | $<1$ | $0.47$
SMACS_z12b | 7:22:52.261 | -73:27:55.497 | $12.26^{+0.17}_{-0.16}$ | $-20.01\pm 0.17$ |  $-2.82\pm$0.21 | $7.91^{+0.26}_{-0.17}$ | $7.56^{+0.30}_{-0.23}$ | $4.0\pm 1.0$ | $1.99$
$z>15$ candidates
SMACS_z16a | 7:23:26.393 | -73:28:04.561 | $15.92^{+0.17}_{-0.15}$ | $-20.59\pm 0.15$ |  $-2.63\pm$0.13 | $8.79^{+0.32}_{-0.33}$ | $7.65^{+0.36}_{-0.39}$ | $1.1\pm 0.2$ | $0.40$
SMACS_z16b | 7:22:39.439 | -73:30:08.185 | $15.32^{+0.16}_{-0.13}$ | $-20.96\pm 0.14$ |  $-2.40\pm$0.34 | $8.80^{+0.44}_{-0.25}$ | $7.46^{+0.52}_{-0.35}$ | $2.8\pm 0.6$ | $1.17$
Figure 6: Cutout images of SMACS_z16a in the JWST bands centered on the
source. The left-most column shows the LW-stack detection image (cf. section
2) and the source partition map. The following columns show the seven NIRCam
and NIRISS science images in the top rows, best-fit model image in the middle
rows and the residual in the bottom rows. The science and model images are
scaled using a linear-stretch normalization, with minimum value equal to zero
and maximum value 10 times the $2\sigma$-clipped standard deviation. The
residuals are scaled linearly between $\pm 7$ times the $2\sigma$-clipped
standard deviation.

Figure 7: Best-fit solution for the SED and photometric redshift of
SMACS_z16a. Upper row: Best-fit SED using the BEAGLE code. Left: Best-fit SED
(black solid curve) with the observed photometric data (blue points) and
expected model photometric points (black points) and associated uncertainties
(pink areas). Right: Triangle plot of the posterior probability distribution
of the four fitted galaxy parameters: redshift, stellar mass, stellar age and
attenuation. Bottom row: Best-fit SED using the EAZY code. Left panel: The
best-fit SED over-plotted over the observed flux densities (in dark squares).
Model flux densities are shown in blue circles. The Lyman-break of the SED of
this galaxy is estimated at $z=15.88$ and the redshift probability
distribution function is shown in the right panel. Both codes agree on a high-
redshift solution with a relatively narrow posterior distribution and which
does not show a secondary peak at lower redshift. Figure 8: Same as Fig. 6,
this time for the $z\sim 12$ candidate galaxy SMACS_z12a.

Figure 9: Same as Fig. 7, here showing the $z\sim 12$ galaxy candidate
SMACS_z12a. Figure 10: Example of a galaxy candidate, SMACS_z10c, at $z\simeq
9.76\pm 0.01$ showing an excess in the F444W-band. The best-fit SED (same
legend as Fig. 7) shows the presence of a Balmer break around 4.5 µm.

The final sample of high-redshift candidates consists of a total of 6 galaxies
in the redshift range $9<z<11$, two galaxies at $11<z<15$, and two galaxies at
$z>15$. All sources are new identifications, since Salmon et al. (2020) only
reported galaxies below $z\sim 8$ in this field. The complete results for our
full sample are reported in Tab. 2.

Our highest redshift galaxy, SMACS_z16a, is consistently identified by the
dropout selection and the two SED fitting procedures, as shown in Fig. 7. The
photometric redshifts derived by both codes (cf. section 5) are in excellent
agreement: $z=15.92^{+0.17}_{-0.15}$ and $z=15.89\pm 0.43$ for BEAGLE and EAZY
respectively. This places the galaxy just $\sim 250$ Myr after the Big Bang.
The galaxy exhibits an extremely blue UV continuum slope of $\beta=-2.63\pm
0.13$ which is consistent with a young and low-attenuation galaxy (cf. Fig.
7). We note in general that no template in EAZY is able to match the very blue
UV continuum, whereas BEAGLE allows for more flexibility in this range of
galaxy templates. To obtain a more accurate estimate of the physical
properties of our galaxies, we ran a second round of SED fitting with BEAGLE
using a prior in redshift from the PDF of EAZY (cf. section 5.1). We derive a
stellar mass of $\log(M_{\star}/\mathrm{M}_{\odot})=8.79^{+0.32}_{-0.33}$ and
a stellar age of $\log(t_{\mathrm{age}}/\mathrm{yr})=7.65^{+0.36}_{-0.39}$ for
SMACS_z16a. This galaxy is also very compact, as can be seen in Fig. 6, with
an estimated size of $\sim 0.4$ kpc based on the model fitting of SE++ (cf.
section 3.1) and after correcting for the lensing distortion.

In the next lower redshift range, one of the two most robust candidates,
SMACS_z12a, is presented in Fig. 8. Again, the source is well detected in all
four filters and shows a clear continuum break. As before, there is a good
agreement between the two photometric redshifts $z=12.20^{+0.21}_{-0.12}$
(BEAGLE) and $z=12.11\pm 0.19$ (EAZY). The PDF from EAZY is broader than the
BEAGLE one, because of the challenge in matching the blue colors
($\beta\sim-2.60\pm 0.32$) of this object (Fig. 9). This source is also among
the highest-redshift candidates identified photometrically in the literature
(e.g. Harikane et al., 2022; Naidu et al., 2022; Castellano et al., 2022),
around 360 Myr after the Big Bang. The derived parameters are very similar to
SMACS_z16a, with a stellar mass of
$\log(M_{\star}/\mathrm{M}_{\odot})=8.14^{+0.21}_{-0.17}$, and an age of
$\log(t_{\mathrm{age}}/\mathrm{yr})=7.61^{+0.21}_{-0.20}$.

In the lowest redshift range, $9<z<11$, some sources around $z\sim 10$ show
indication of a Balmer break as a clear excess in the F444W-band. An example
is shown in Fig. 10, where flux excess is observed relative to a flat
continuum in the bluer bands. This could indicate that an evolved stellar
population was already in place at those early epochs, which is confirmed by
the best-fit SED with a stellar age around $\sim 400$ Myr. An example for our
$9<z<11$ objects is shown in appendix B. Of course, a stronger contribution
from rest-frame optical emission lines is also possible, although this
alternative does not yield the best-fit solution.

##  6 Discussion

Figure 11: Absolute UV magnitude $M_{\mathrm{UV}}$ as a function of redshift
for the high-redshift candidates in this work, along with a compilation of
known galaxies at $z>8$. The candidates in this work are shown as orange
filled hexagons, while empty hexagon marks the less secure candidate. The
empty purple pentagons show the recently identified candidates from JWST in
the GLASS parallel field by Naidu et al. (2022). The empty blue circles show
photometric-redshift galaxies, while the solid circles represent
spectroscopically confirmed galaxies compiled from Bouwens et al. (2022).

While observing through a lensing cluster offers invaluable flux
amplification, it also results in a smaller survey area in the source plane.
The field-of-view in the NIRCam image plane is about $9.7$ arcmin2, which
translates into $7.23$ arcmin2 in the source plane after de-lensing. Although
the area reduction is still reasonable (because one of the NIRCam modules is
less affected by magnification), wide-area surveys are more prone to uncover
rare and bright high-redshift galaxies (e.g. Roberts-Borsani et al., 2016;
Oesch et al., 2016; Bowler et al., 2020; Harikane et al., 2022), while
lensing-assisted surveys are sensitive to the fainter population. While some
of our sources are relatively bright, they all have luminosities around, or
below, the characteristic luminosity $M^{\star}$ at $z\sim 9-10$ (Oesch et
al., 2018; Bouwens et al., 2021). The highest redshift candidate in our
sample, SMACS_z16a, has an intrinsic UV magnitude of
$M_{\mathrm{UV}}=-20.59\pm 0.15$, which is fainter than previous very high
redshift candidates Gn_z11 and HD1, and roughly similar to the recently
discovered candidate GLz13 (Naidu et al., 2022). Moreover, the yield of the
galaxy hunt presented in this paper is surprisingly high, with the discovery
of 2 candidates at $z\sim 16$, 2 candidates at $z\sim 12$, and 6 candidates at
$z=10-11$. In Fig. 11, we put these discoveries in the context of known high-
redshift galaxies and compare our candidates to the published samples of
galaxies at $z>8$. Not only do these results demonstrate the exceptional
capabilities of JWST, but they also suggest a higher number density than
predicted by an extrapolation of the $z\sim 9-10$ luminosity function (e.g.
Bouwens et al., 2021) and theoretical predictions (e.g., Ocvirk et al., 2020;
Dayal et al., 2022). Some early theoretical interpretations suggest that these
observations could be compatible with a total absence of dust in these
galaxies, as their very blue UV slopes might also indicate (Ziparo et al.,
2022). The evolution of the UV luminosity function observed through lensing
clusters will be the scope of a future paper.

The inferred age of SMACS_z16a indicates a recent formation history, which
points to a rapid stellar mass assembly in this galaxy. In the context of
early galaxy formation, one would expect the first galaxies to more resemble
SMACS_z16a in terms of physical parameters, than the bright and uncommon
galaxies discussed earlier. As we are closing in on the formation epoch of the
first galaxies, it will become increasingly challenging to find such evolved
and relatively massive objects. Perhaps upcoming and future deep surveys
through gravitational lenses will offer the best route to unveil these
smallest structures. This first JWST cycle already includes several ERS,
General Observer (GO), and Guaranteed Time (GTO) programs that target lensing
clusters using different observing modes.

Another interesting finding is the small sizes of these galaxies. For
instance, SMACS_z16a has an effective radius of $r_{\rm e}=0.4$ kpc, which is
consistent with recent size measurement of high-redshift galaxies (Kawamata et
al., 2018; Adams et al., 2022). More importantly, the Sérsic index, derived
from the SE++ modeling and source-plane reconstruction, hints to a disk-like
morphology. We measured $n=1.1\pm 0.2$ for SMACS_z16a, which corresponds to a
disk-like profile, although there is a real possibility that we are detecting
only the brightest star-forming region in the galaxy. Such a morphology was
also observed for GLz13 (Naidu et al., 2022), as well as in several other
high-redshift candidates reported in Adams et al. (2022). The existence of
well settled disks at early epochs would provide strong constraints on
theoretical models of galaxy formation.

While the nature of these ERO observations implies potential sources of
uncertainties, we performed photometric quality verifications, and also
combined multiple selection methods to build our sample of high-redshift
sources. The dropout selection and the photometric redshifts from two
independent codes significantly increase the robustness of the candidates. In
the near future, deeper imaging in the short wavelength range will bring
additional confirmation to these candidates, in addition to follow-up
spectroscopy, e.g with NIRSpec, which will provide definitive redshifts.

##  7 Summary

In this paper, we presented our search for very high-redshift galaxies in the
first JWST observations of the lensing cluster SMACS 0723. We report the
discovery of the highest redshift candidate known to date at $z\sim 16$. We
also discover two candidates at $z\sim 12$, and 11 candidates at $z\sim
10-11$.

We have used 6 NIRCam broadband filters, and one NIRISS filter, between $\sim
0.8$ µm and 5 µm to identify Lyman break galaxies in the redshift ranges
$9<z<11$, $11<z<15$, and $16<z<21$. We combined the dropout selection with
SED-fitting to obtain accurate photometric redshifts using two different
codes, BEAGLE and EAZY. We also simultaneously derive the physical properties
of these galaxies. Our final sample is based on all candidates that satisfy
the three selection methods and whose photometric redshifts from the two SED
fits are in agreement. To determine the intrinsic luminosity of galaxies and
account for potential multiple images, we use the one of the latest lensing
models of SMACS 0723 produced from the same JWST data set in a companion paper
(Mahler et al., 2022). The measured amplification factors range from 1 to
$\sim 4$, which gives us a lensing boost up to 1.5 magnitudes in the best
case.

With an intrinsic UV luminosity of $M_{\mathrm{UV}}=-20.59\pm 0.15$, we find
that SMACS_z16a is fainter than typical $z>10$ galaxies in the literature and
close to the characteristic luminosity $M^{\star}$ at $z\sim 10$. The best-fit
SED of SMACS_z16a is compatible with a young age of
$\log(t_{\mathrm{age}}/\mathrm{yr})=7.65^{+0.36}_{-0.39}$ and a stellar mass
of $\log(M_{\star}/\mathrm{M}_{\odot})=8.79^{+0.32}_{-0.33}$. The bulk of its
stellar mass build-up took place in the last 20 Myr. This galaxy, like some
other candidates, shows a very blue UV continuum-slope of $\beta=-2.63\pm
0.13$ consistent with young ages. In general, the galaxy candidates, from
$\sim 10$ to $z\sim 16$, have ages in the range
$\log(t_{\mathrm{age}}/\mathrm{yr})\sim=7.3-8.7$, and stellar masses of
$\log(M_{\star}/\mathrm{M}_{\odot})\sim 7.8-10.2$. Our source extraction and
photometry procedure with SE++ also provide morphological measurements, since
it relies on model-fitting. The measured sizes were reconstructed back to the
lens plane using lenstool to determine the intrinsic effective radius and the
Sersić index of the sources. Our measurements indicate very compact sources,
with $r_{\rm e}$ below 1 kpc for most of the candidates. Some of these sources
also have Sérsic indices compatible with disk-like a morphology. If true, it
is surprising, for a galaxy like SMACS_z16a, to evolve so quickly in the span
of less than $\sim 250$ Myr.

The discovery of one, possibly two, galaxy candidates only $\sim 250$ Myr
after the Big Bang, together with eight new candidates at $9<z<13$, announces
great discoveries awaiting in the deep extra-galactic fields that will be
observed by JWST. The number density, the physical properties and the
morphologies of these galaxies will be extremely valuable in constraining
hydrodynamical simulations of galaxy formation in a cosmological context.

## Acknowledgements

We would like to thank Gabriel Brammer for the grizli reduction notebooks. HA
acknowledges support from CNES. LF and AZ acknowledge support by Grant No.
2020750 from the United States-Israel Binational Science Foundation (BSF) and
Grant No. 2109066 from the United States National Science Foundation (NSF). AZ
acknowledges support by the Ministry of Science & Technology, Israel.

This work is based on observations obtained with the NASA/ESA/CSA JWST and the
NASA/ESA Hubble Space Telescope (HST), retrieved from the Mikulski Archive for
Space Telescopes (MAST) at the Space Telescope Science Institute (STScI).
STScI is operated by the Association of Universities for Research in
Astronomy, Inc. under NASA contract NAS 5-26555. This work has made use of the
CANDIDE Cluster at the Institut d’Astrophysique de Paris (IAP), made possible
by grants from the PNCG and the region of Île de France through the program
DIM-ACAV+.

The JWST Early Release Observations (ERO) and associated materials were
developed, executed, and compiled by the ERO production team: Hannah Braun,
Claire Blome, Matthew Brown, Margaret Carruthers, Dan Coe, Joseph De-
Pasquale, Nestor Espinoza, Macarena Garcia Marin, Karl Gordon, Alaina Henry,
Leah Hustak, Andi James, Ann Jenkins, Anton Koekemoer, Stephanie LaMassa,
David Law, Alexandra Lockwood, Amaya Moro-Martin, Susan Mullally, Alyssa
Pagan, Dani Player, Klaus Pontoppidan, Charles Prof- fitt, Christine Pulliam,
Leah Ramsay, Swara Ravindranath, Neill Reid, Massimo Robberto, Elena Sabbi,
Leonardo Ubeda. The EROs were also made possible by the foundational efforts
and support from the JWST instruments, STScI planning and scheduling, and Data
Management teams.

## Data Availability

The data underlying this article are publicly available on the Mikulski
Archive for Space Telescopes222https://archive.stsci.edu/ (MAST), under
program ID 2736.

## References

  * Adams et al. (2022) Adams N. J., et al., 2022, arXiv e-prints, p. arXiv:2207.11217
  * Allard et al. (2001) Allard F., Hauschildt P. H., Alexander D. R., Tamanai A., Schweitzer A., 2001, ApJ, 556, 357
  * Atek et al. (2011) Atek H., et al., 2011, ApJ, 743, 121
  * Atek et al. (2014a) Atek H., et al., 2014a, ApJ, 786, 60
  * Atek et al. (2014b) Atek H., et al., 2014b, ApJ, 789, 96
  * Atek et al. (2015) Atek H., et al., 2015, ApJ, 800, 18
  * Atek et al. (2018) Atek H., Richard J., Kneib J.-P., Schaerer D., 2018, MNRAS, 479, 5184
  * Bacon et al. (2010) Bacon R., et al., 2010, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 773508, doi:10.1117/12.856027
  * Bertin et al. (2020) Bertin E., Schefer M., Apostolakos N., Álvarez-Ayllón A., Dubath P., Kümmel M., 2020, in Pizzo R., Deul E. R., Mol J. D., de Plaa J., Verkouter H., eds, Astronomical Society of the Pacific Conference Series Vol. 527, Astronomical Data Analysis Software and Systems XXIX. p. 461
  * Bhatawdekar et al. (2019) Bhatawdekar R., Conselice C. J., Margalef-Bentabol B., Duncan K., 2019, MNRAS, 486, 3805
  * Bouwens et al. (2015) Bouwens R. J., et al., 2015, ApJ, 803, 34
  * Bouwens et al. (2017) Bouwens R. J., Oesch P. A., Illingworth G. D., Ellis R. S., Stefanon M., 2017, ApJ, 843, 129
  * Bouwens et al. (2021) Bouwens R. J., et al., 2021, AJ, 162, 47
  * Bouwens et al. (2022) Bouwens R. J., et al., 2022, ApJ, 931, 160
  * Bowler et al. (2020) Bowler R. A. A., Jarvis M. J., Dunlop J. S., McLure R. J., McLeod D. J., Adams N. J., Milvang-Jensen B., McCracken H. J., 2020, MNRAS, 493, 2059
  * Boyer et al. (2022) Boyer M. L., et al., 2022, arXiv e-prints, p. arXiv:2209.03348
  * Brammer et al. (2008) Brammer G. B., van Dokkum P. G., Coppi P., 2008, ApJ, 686, 1503
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, Monthly Notices of the Royal Astronomical Society, 344, 1000
  * Bruzual A. & Charlot (1993) Bruzual A. G., Charlot S., 1993, ApJ, 405, 538
  * Capak et al. (2015) Capak P. L., et al., 2015, Nature, 522, 455
  * Castellano et al. (2022) Castellano M., et al., 2022, arXiv e-prints, p. arXiv:2207.09436
  * Chabrier et al. (2000) Chabrier G., Baraffe I., Allard F., Hauschildt P., 2000, ApJ, 542, 464
  * Chevallard & Charlot (2016) Chevallard J., Charlot S., 2016, MNRAS, 462, 1415
  * Coe et al. (2019) Coe D., et al., 2019, ApJ, 884, 85
  * Conroy & Gunn (2010) Conroy C., Gunn J. E., 2010, ApJ, 712, 833
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * Dayal et al. (2022) Dayal P., et al., 2022, MNRAS, 512, 989
  * De Barros et al. (2019) De Barros S., Oesch P. A., Labbé I., Stefanon M., González V., Smit R., Bouwens R. J., Illingworth G. D., 2019, MNRAS, 489, 2355
  * Doyon et al. (2012) Doyon R., et al., 2012, in Clampin M. C., Fazio G. G., MacEwen H. A., Oschmann Jacobus M. J., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave. p. 84422R, doi:10.1117/12.926578
  * Duncan et al. (2014) Duncan K., et al., 2014, MNRAS, 444, 2960
  * Ebeling et al. (2001) Ebeling H., Edge A. C., Henry J. P., 2001, ApJ, 553, 668
  * Elíasdóttir et al. (2007) Elíasdóttir Á., et al., 2007, preprint,  (arXiv:0710.5636)
  * Euclid Collaboration et al. (2022a) Euclid Collaboration et al., 2022a, arXiv e-prints, p. arXiv:2209.12906
  * Euclid Collaboration et al. (2022b) Euclid Collaboration et al., 2022b, arXiv e-prints, p. arXiv:2209.12907
  * Eyles et al. (2007) Eyles L. P., Bunker A. J., Ellis R. S., Lacy M., Stanway E. R., Stark D. P., Chiu K., 2007, MNRAS, 374, 910
  * Ferland et al. (2013) Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
  * Finkelstein et al. (2015) Finkelstein S. L., et al., 2015, ApJ, 810, 71
  * Fox et al. (2022) Fox C., Mahler G., Sharon K., Remolina González J. D., 2022, ApJ, 928, 87
  * Furtak et al. (2021) Furtak L. J., Atek H., Lehnert M. D., Chevallard J., Charlot S., 2021, MNRAS, 501, 1568
  * Golubchik et al. (2022) Golubchik M., Furtak L. J., Meena A. K., Zitrin A., 2022, arXiv e-prints, p. arXiv:2207.05007
  * González et al. (2011) González V., Labbé I., Bouwens R. J., Illingworth G., Franx M., Kriek M., 2011, ApJ, 735, L34
  * Grazian et al. (2015) Grazian A., et al., 2015, A&A, 575, A96
  * Gutkin et al. (2016) Gutkin J., Charlot S., Bruzual G., 2016, MNRAS, 462, 1757
  * Harikane et al. (2022) Harikane Y., et al., 2022, ApJ, 929, 1
  * Inoue et al. (2014) Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805
  * Ishigaki et al. (2018) Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, ApJ, 854, 73
  * Jakobsen et al. (2022) Jakobsen P., et al., 2022, A&A, 661, A80
  * Jiang et al. (2021) Jiang L., et al., 2021, Nature Astronomy, 5, 256
  * Jullo & Kneib (2009) Jullo E., Kneib J.-P., 2009, MNRAS, 395, 1319
  * Jullo et al. (2007) Jullo E., Kneib J.-P., Limousin M., Elíasdóttir Á., Marshall P. J., Verdugo T., 2007, New Journal of Physics, 9, 447
  * Kauffmann et al. (2022) Kauffmann O. B., et al., 2022, arXiv e-prints, p. arXiv:2207.11740
  * Kawamata et al. (2015) Kawamata R., Ishigaki M., Shimasaku K., Oguri M., Ouchi M., 2015, ApJ, 804, 103
  * Kawamata et al. (2018) Kawamata R., Ishigaki M., Shimasaku K., Oguri M., Ouchi M., Tanigawa S., 2018, ApJ, 855, 4
  * Kikuchihara et al. (2020) Kikuchihara S., et al., 2020, ApJ, 893, 60
  * Kinney et al. (1996) Kinney A. L., Calzetti D., Bohlin R. C., McQuade K., Storchi-Bergmann T., Schmitt H. R., 1996, ApJ, 467, 38
  * Kneib et al. (1996) Kneib J. P., Ellis R. S., Smail I., Couch W. J., Sharples R. M., 1996, ApJ, 471, 643
  * Kümmel et al. (2020) Kümmel M., Bertin E., Schefer M., Apostolakos N., Álvarez-Ayllón A., Dubath P., 2020, in Pizzo R., Deul E. R., Mol J. D., de Plaa J., Verkouter H., eds, Astronomical Society of the Pacific Conference Series Vol. 527, Astronomical Data Analysis Software and Systems XXIX. p. 29
  * Livermore et al. (2017) Livermore R. C., Finkelstein S. L., Lotz J. M., 2017, ApJ, 835, 113
  * Lotz et al. (2017) Lotz J. M., et al., 2017, ApJ, 837, 97
  * Mahler et al. (2022) Mahler G., et al., 2022, arXiv e-prints, p. arXiv:2207.07101
  * McLure et al. (2011) McLure R. J., et al., 2011, MNRAS, 418, 2074
  * Montes (2022) Montes M., 2022, Nature Astronomy, 6, 308
  * Naidu et al. (2022) Naidu R. P., et al., 2022, arXiv e-prints, p. arXiv:2207.09434
  * Ocvirk et al. (2020) Ocvirk P., et al., 2020, MNRAS, 496, 4087
  * Oesch et al. (2016) Oesch P. A., et al., 2016, ApJ, 819, 129
  * Oesch et al. (2018) Oesch P. A., Bouwens R. J., Illingworth G. D., Labbé I., Stefanon M., 2018, ApJ, 855, 105
  * Oke & Gunn (1983) Oke J. B., Gunn J. E., 1983, ApJ, 266, 713
  * Pei (1992) Pei Y. C., 1992, ApJ, 395, 130
  * Polletta et al. (2007) Polletta M., et al., 2007, ApJ, 663, 81
  * Reddy et al. (2015) Reddy N. A., et al., 2015, ApJ, 806, 259
  * Reddy et al. (2018a) Reddy N. A., et al., 2018a, ApJ, 853, 56
  * Reddy et al. (2018b) Reddy N. A., et al., 2018b, ApJ, 869, 92
  * Repp & Ebeling (2018) Repp A., Ebeling H., 2018, MNRAS, 479, 844
  * Richard et al. (2014) Richard J., et al., 2014, preprint,  (arXiv:1405.3303)
  * Rigby et al. (2022) Rigby J., et al., 2022, arXiv e-prints, p. arXiv:2207.05632
  * Roberts-Borsani et al. (2016) Roberts-Borsani G. W., et al., 2016, ApJ, 823, 143
  * Salmon et al. (2020) Salmon B., et al., 2020, ApJ, 889, 189
  * Schaerer & de Barros (2009) Schaerer D., de Barros S., 2009, A&A, 502, 423
  * Schaerer & de Barros (2010) Schaerer D., de Barros S., 2010, A&A, 515, A73
  * Schaerer et al. (2022) Schaerer D., Marques-Chaves R., Oesch P., Naidu R., Barrufet L., Izotov Y. I., Guseva N. G., Brammer G., 2022, arXiv e-prints, p. arXiv:2207.10034
  * Sérsic (1963) Sérsic J. L., 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Shivaei et al. (2020) Shivaei I., et al., 2020, ApJ, 899, 117
  * Silva et al. (1998) Silva L., Granato G. L., Bressan A., Danese L., 1998, ApJ, 509, 103
  * Smit et al. (2014) Smit R., et al., 2014, ApJ, 784, 58
  * Song et al. (2016) Song M., et al., 2016, ApJ, 825, 5
  * Stark et al. (2009) Stark D. P., Ellis R. S., Bunker A., Bundy K., Targett T., Benson A., Lacy M., 2009, ApJ, 697, 1493
  * Stefanon et al. (2021) Stefanon M., Bouwens R. J., Labbé I., Illingworth G. D., Gonzalez V., Oesch P. A., 2021, ApJ, 922, 29
  * Steidel et al. (1996) Steidel C. C., Giavalisco M., Pettini M., Dickinson M., Adelberger K. L., 1996, The Astrophysical Journal, 462, L17
  * Ziparo et al. (2022) Ziparo F., Ferrara A., Sommovigo L., Kohandel M., 2022, arXiv e-prints, p. arXiv:2209.06840
  * de Barros et al. (2014) de Barros S., Schaerer D., Stark D. P., 2014, A&A, 563, A81

##  Appendix A Photometry validation

Figure 12: Magnitude number counts in the HST WFC3/IR F140W and JWST NIRCam
F150W bands. The magnitudes are obtained for objects detected in the shallower
F140W-band and measured simultaneously in F140W and F150W. The number counts
show a good agreement and no photometric offset is present. Figure 13:
Magnitude comparison between HST-WFC3/IR F140W and JWST-NIRCam F150W bands.
The blue line marks the running median offset and the blue envelope the
$\sigma$-clipped standard deviation. The magnitudes of the two bands show a
relatively tight dispersion and an offset of about 0.24 mag, that can be
explained by the difference between the filter curves and pivot wavelengths.

We perform photometry validation by comparing our photometry with the well
known and calibrated HST WFC3/IR images from the RELICS survey. We run SE++ on
a stacked detection image and do photometric measurements simultaneously on
the WFC3/IR F140W and NIRCam F150W bands.

Fig. 12 shows the AUTO magnitude number counts in NIRCam F150W and WFC3/IR
F140W as measured with SE++. The number counts are in excellent agreement,
which is not surprising given the same detection image and with minimal or no
offsets. This confirms the accurate zero-point calibration of our JWST data
reduction, explained in section 2. The comparison of magnitudes is shown in
Fig 2, where on the $y-$axis we show the magnitude difference between the
F140W and F150W bands. There is a relatively tight scatter and a median offset
of about $0.36$ towards fainter F140W magnitudes. Most of this offset can be
attributed to the difference between the F140W and F150W filter curves and
pivot wavelengths. Using Pysynphot, we performed synthetic photometry in the
two filters, assuming a variety of galaxy templates from BC95 (Bruzual A. &
Charlot, 1993) and KC96 (Kinney et al., 1996) libraries, and a redshift range
of $z=0-5$. We found a median offset of $\sim 0.3$ magnitudes.

##  Appendix B Example of a $z\sim 10$ candidate

Figure 14: Same as Fig. 6, this time for the $z\sim 12$ candidate galaxy
SMACS_z11a.

Figure 15: Same as Fig. 7, here showing the galaxy candidate SMACS_z11a.

In addition to the $z\sim 16$ and $z\sim 12$ examples shown in the main text
(cf. section 5.2), we display here an example of a $z\sim 10$ galaxy, the
candidate SMACS_z11a, in Figs 14 and 15.

##  Appendix C Full candidate list

Figure 16: Photometric data and best-fit SED results for candidate SMACS_z10a.

Figure 17: Photometric data and best-fit SED results for candidate SMACS_z10b.

Figure 18: Photometric data and best-fit SED results for candidate SMACS_z10c.

Figure 19: Photometric data and best-fit SED results for candidate SMACS_z10d.

Figure 20: Photometric data and best-fit SED results for candidate SMACS_z10e.

Figure 21: Photometric data and best-fit SED results for candidate SMACS_z12b.

Figure 22: Photometric data and best-fit SED results for candidate SMACS_z16b.