# Can self-training identify suspicious ugly duckling lesions?

Mohammadreza Mohseni Corresponding Author School of Computing Science, Simon
Fraser University  MetaOptima Technology Inc  Jordan Yap  MetaOptima
Technology Inc  William Yolland  MetaOptima Technology Inc  Arash Koochek
Banner Health  M Stella Atkins  School of Computing Science, Simon Fraser
University  Department of Skin Science and Dermatology, University of British
Columbia

###### Abstract

One commonly used clinical approach towards detecting melanomas recognises the
existence of Ugly Duckling nevi, or skin lesions which look different from the
other lesions on the same patient. An automatic method of detecting and
analysing these lesions would help to standardize studies, compared with
manual screening methods. However, it is difficult to obtain expertly-labelled
images for ugly duckling lesions. We therefore propose to use self-supervised
machine learning to automatically detect outlier lesions. We first
automatically detect and extract all the lesions from a wide-field skin image,
and calculate an embedding for each detected lesion in a patient image, based
on automatically identified features. These embeddings are then used to
calculate the L2 distances as a way to measure dissimilarity. Using this deep
learning method, Ugly Ducklings are identified as outliers which should
deserve more attention from the examining physician. We evaluate through
comparison with dermatologists, and achieve a sensitivity rate of 72.1% and
diagnostic accuracy of 94.2% on the held-out test set.

##  1 Introduction

In order to diagnose skin cancers such as malignant melanoma, careful visual
inspection of all the patient’s skin lesions, particularly melanocytic nevi,
must be undertaken. To this end, dermatologists have introduced several
approaches which can be helpful in detecting melanomas, including the well-
known ABCD dermatology criteria (Asymmetry, Border irregularity, Color
irregularity, Diameter$>$6mm)[26, 1] and 7-point checklist[4, 5]. Another
useful clinical indicator is the existence of Ugly Duckling (UD) lesions,
which look different from the other lesions on the same patient[12]. The
existence of a UD lesion is strongly correlated with the existence of melanoma
[11, 33], although there is considerable inter-observer variability in
visually selecting UD lesions from whole body imaging [16, 32]. After
identification of UDs, clinically concerning lesions can then be inspected
with a close-up examination, for example under dermatoscopic inspection.
However, converting the clinical view of a small 1.5mm diameter nevus into a
suitable digital representation for imaging and automatically identifying
suspicious lesions is challenging.

Automatic UD detection can be treated as a form of outlier detection problem
whose data comes from wide-field (clinical) images of the same patient,
typically acquired during total body photography (TBP) also known as full body
imaging. In TBP, overview images are typically taken with cameras placed 20-50
cm away from the skin surface, showing enough context to identify the body
part, mimicking the visual inspection made by a clinician.

However a limitation of UD analysis from TBP images is lesion size. With
current camera resolutions, very small lesions less than 1.5mm diameter cannot
be readily detected in a TBP image of one body part. And just like a physical
examination, a closer-up view is needed for lesion diagnosis.

An example of a typical TBP image taken with a smartphone camera from about
50cm away is shown in Fig 1. The filesize after this image has been compressed
using JPEG is 187 KB.

Figure 1: Typical image of part of a back, 1640 x 1116 pixels, used with
permission from the SD-260 dataset [39]

Converting the clinical view of a small 1.5mm diameter nevus into a suitable
digital representation for imaging and diagnosing suspicious nevi is
challenging. In our experience, the lesion should occupy an area of at least
20x20 pixels to be reliably detected by algorithms and for automated analysis.
Therefore, a small 1.5mm diameter lesion requires resolution of around
0.075mm/pixel at the imaging surface. A wide-field digital image taken at
about 20cm from the subject with a 10 MP smartphone camera, allows digital
resolution of around 0.075mm/pixel [9]; smaller lesions cannot be reliably
detected or analysed [34]. Artifacts from reflections, lack of polarization,
image compression, motion blurring, lighting, and shadows all contribute to
the challenge of both identifying UDs from a single digital wide-field image,
as well as in detecting changes in lesions found in a series of TBP images
taken over time.

Difficulties in obtaining high enough resolution TBP images from personal
cameras such as smartphones for automated analysis of small lesions, has also
led to difficulties in obtaining expertly-labelled digital TBP images. In our
experience, the expert reader cannot zoom the image enough to make an
analysis, unlike in a clinical environment where the physician can look closer
at suspicious lesions using a magnifying dermoscope. For example, dermoscopic
images are shown in Fig 2.The size of these dermoscopic images are 6000 x 4000
pixels, JPEG Size 1 MB, providing much more detail for machine learning
classification algorithms.

Even very small lesions such as in Fig 2(a), where the nevus is about 3mm
wide, occupy about 30,000 pixels.

(a) Dermoscopic nevus

(b) Dermoscopic melanoma

Figure 2: Two dermoscopic images from the same patient, each of size 6000x4000
pixels (publicly available skin images from ISIC [37] ). Note the mm scale
ruler in the top left corner of each image.

A few commercial systems can provide automated analysis of digital TBP images,
where the images are taken with several fixed cameras for super-high
resolution. One is the Canfield 3D Vectra system, where the lesion images are
automatically stitched together, and then suspicious lesions can be viewed
clinically with a dermoscope [36]. Another system which supports TBP using
mobile cameras is provided by Dermengine [8]. An overview of such systems is
provided by the International Society for Digital Imaging of the Skin (ISDIS)
[ISDIS 2021 [18]].

Automated analysis of TBP "wide-field" digital images taken with smartphones
remains challenging, because of the difficulty of classifying and labelling
such small lesions in the images. We hypothesised that outlier and suspicious
UD lesions could be identified using a self-supervised learning method,
leading to identification of suspicious lesions which can then be chosen for
further detailed analysis (by computer algorithm and/or examination by a
board-certified dermatologist).

In this paper we describe our objective, unbiased, and reproducible method to
automatically identify UDs in TBP images, achieving accuracy comparable to
that of expert dermatologists. By using our approach, expert dermatologists
can be helped to identify which lesions deserve further examination in
patients with numerous atypical lesions. These tasks are repetitive and time-
consuming, with poor intra-observer and inter-observer reliability and
consistency. Furthermore, as Schlessinger et al.suggested [31], our method can
also be used in objective tasks needed for various outcome measures in
clinical trials involving TBP images.

Figure 3: UD Detection Pipeline

##  2 Related Work

With the recent advances in deep learning and its promising results, deep
learning found its way to various areas of science and dermatology was no
exception. Esteva et al.showed how promising deep learning can be in skin
cancer classification [10]. Liu et al.extended this success to classification
of a large number of general dermatology and skin cancer classes [24]. The
contribution of deep learning has not been limited to classification. It has
also contributed to segmentation and detection tasks.

One still unsolved challenge in deep learning is outlier detection. There have
been many solutions proposed to overcome this challenge but no final solution
is known yet. In the literature, this problem is known by different names such
as outlier detection, anomaly detection, or one-class classification. A large
body of research has demonstrated that finding a good lower dimension
representation can lead to promising solutions to outlier detection problems
[30, 7, 2]. Perera et al.showed that representations learned for one-class
classification, can be transferred from an unrelated task[28]. Autoencoders
are a family of neural networks which learn low-dimensional representation
from the input domain. Variational autoencoders (VAE) are a subtype of
autoencoders which control the latent space by enforcing a prior [20]. Higgins
et al.showed that VAEs are able to learn disentangled features[14] and
introduced parameter $\beta$ and $\beta$-VAE as a means for disentangled
representation learning[15].

Despite the recent advances in machine learning and deep learning, UD lesion
detection still remains a challenge in computer aided dermatology. Soenksen et
al.recently published a supervised machine learning method to identify UD
lesions using the features of a pretrained network [34, 6]. Three board-
certified dermatologists labelled UDs by rank order in 135 TBP patient images.
Their results show between 83% and 88% agreement between the top 3 UDs ranked
by the dermatologists and the top 3 outlier lesions ranked by the algorithm.

In this paper we show that it is possible to skip the long process of
collecting labelled data and training a classifier. We propose self-training
as an approach which enables researchers to reproduce our results easily; and
comparably does as good job as using the labelled data.

##  3 Method

In order to detect and classify each lesion as an inlier or UD, we developed a
multi-stage pipeline which consists of three main modules (shown in Fig 3).
First, we pass the input TBP image to the detection module. This module
detects all the lesions in the TBP image, crops a fixed sized window centered
on the lesion and passes them to the segmentation module. The segmentation
module segments the lesions from the skin and prepares the segmented images
for the outlier detection module which is the last step.

###  3.1 Detection

The detection module used was based on the Single Shot MultiBox Detector (SSD)
[22]. We used Resnet-18 [13] as our backbone along with a Feature Pyramid
Network (FPN) [21] for improved multi-scale detection performance. The
detection network was trained using 512x512 images as input to the network.

To detect lesions across a large TBP image, first the full-sized image is
split up into 512x512 tiles with a 50% overlap. Tiles are then sent through
the detection module to obtain the predicted locations of lesions. Detections
are then aggregated across all tiles and non-max suppression is performed to
obtain the final set of lesion detections for the entire TBP image.

###  3.2 Segmentation

The segmentation module is a smaller variant of the U-Net architecture[29],
with only 6.2% of the trainable parameters as the original model. We use no
batch normalization[17] in the model, similar to the original paper, and
unlike some of the recent implementations in deep learning libraries[27]. The
model runs on 64x64 RGB images, cropped around detected lesions, and produces
a 64x64 probability map. An optimal threshold is determined on a validation
set for binary lesion/background segmentation.

The model is trained on thousands of manually segmented lesions along with a
collection of skin-only patches all of which are sourced from a set of
internal full-body images spanning various skin-types, and body locations. We
found that without the true negative patches, the model produced false
positive segmentation masks when presented with detection errors (skin-only
patches). We train with random rotations, flips, crops, and color-jitter, and
optimize for a pixel-wise cross-entropy loss using the Adam[19] optimizer with
a fixed step-down schedule and an initial learning rate of $1e^{-3}$.

###  3.3 Outlier Detection

Variational autoencoders tend to have higher reconstruction loss on anomalous
samples[3, 25, 23]. This fact suggests that by training a VAE on the extracted
lesions from one patient, UDs end up with high reconstruction loss value. The
intuition behind this statement is that VAEs try to reconstruct common looking
lesions (the majority) as perfectly as possible. UD lesions (which are the
minority) are less prioritized, thus leading to higher reconstruction losses.

Figure 4: An overview of a Variational Autoencoder (VAE) architecture (image
taken from [38])

In our first experiments we self-trained a VAE on extracted lesion images of
one patient, and by defining a threshold on the reconstructed losses, we
classified each lesion as either common looking or as a UD. In order to train
the VAE, we enforce it to reconstruct the input image in the output as shown
in Fig 4.

This approach’s limitation is that for each TBP image query, the VAE needs to
be trained for almost 130 epochs which takes $\approx 2$ minutes. We tried to
overcome this issue by doing pre-training to speed things up.

###  3.4 Fast Outlier Detection

We observed when training the VAE from scratch, in the initial epochs VAE is
learning the basics of describing the lesions as shown in Fig 5(a). With the
increase in the number of epochs, our model learns more about the input domain
and reconstructs better quality images in the output (as shown in Fig 5(b)).
In order to decrease the time spent for the model to learn the basic features,
we self-trained a base VAE model on 300 TBP images. In the query time when a
TBP image is given, we fine-tune the base model on the lesion images from the
TBP image for a few epochs and then calculate reconstruction loss and features
from lesion images. Using this approach, we were able to identify UDs in just
a few seconds and get results as accurate as before.

(a) Reconstructed lesion images at Epoch 13

(b) Reconstructed lesion images at Epoch 121

Figure 5: Two unordered sets of reconstructed images at different epochs. The
model has more information regarding input domain in epoch 121 compared to
epoch 13.

Additionally, we observed that by incorporating a $\beta$-VAE’s disentangled
latent space and calculating the distance between each lesion’s extracted
features and mean of features belonging to the lesions in the same TBP image,
we still can identify UDs and optionally can skip the finetuning phase.

##  4 Data

We evaluated our UD detection algorithm on 75 TBP images in total. 32 images
were sourced from the SD-198 [35], and SD-260 datasets [39], and an additional
43 were collected internally from various clinics. UD lesions in these images
were labelled by a board certified dermatologist. TBP images used in this
experiment contained a varying number of lesions ranging from 10 to 182.
Overall, we extracted a total of 4628 lesion images across all 75 TBP images.
53 of the images contained at least one lesion labelled as UD, with an average
of 1.44 UDs on each TBP image. The validity of data was manually checked after
the segmentation module.

Figure 6: Automatic analysis of all the lesions in the image, with an ordered
list of all the distance rankings, highest at top left hand side. The outliers
which had a distance higher than our defined threshold of 38.63, are
highlighted in blue on the image. The mean and standard deviation of the
distances were 21.60 and 17.03 respectively.

##  5 Results

The UD detection problem can be considered either as ranking prediction
problem or binary classification problem. In order to evaluate our algorithm’s
performance, we used ranking evaluation metrics and binary classification
evaluation metrics. In order to find the rankings, first we calculated an
embedding for each lesion and then calculated the L2 distance between each
embedding to the mean of embeddings. By sorting these distances, we generated
the rankings. In the binary classification setting, we defined a threshold on
each lesion in order to determine whether a lesion should be called UD or not.
Using the following formula, we obtained the threshold for each TBP image:

$threshold=mean(distances)+\min\begin{cases}mean(distances)\\\
std(distances)\end{cases}$ (1)

Fig 6 shows the results of automated analysis of all the lesions on the image
in Fig 1. The 7 lesions outlined in blue rectangles are the outliers according
to our threshold.

###  5.1 Ranking Evaluation Metrics

We used average precision (AP) and reciprocal rank (RR) as our primary
metrics. Additionally we used top-3 and top-7 agreement (Agr.) from [34]. We
calculated each metric for all TBP images and report the average over all TBP
images. In order to calculate top-3 agreement, for each lesion labelled as UD,
we found its rank in our results. If the ranking was less than or equal to 3,
we called success (TBP image counted as 1), otherwise we called it failure
(TBP image counted as 0). We calculate and report the average of success over
TBP images. Top-7 agreement is also calculated in a similar manner. Since not
every TBP image contains a UD lesion, we only calculate top-3 and top-7
agreement for the TBP images containing at least one lesion labelled as UD.
Results are shown in Table 1.

MAP | MRR | Top-3 Agr. | Top-7 Agr.
---|---|---|---
0.659 | 0.721 | 86.79% | 94.34%
Table 1: UD Ranking Evaluation

###  5.2 Binary Classification Metrics

Binary classification predictions can be obtained by applying a threshold on
the distances calculated for each lesion.Results are shown in Table 2. On
average our model predicts $4.25$ UDs per TBP image. In order to evaluate our
method, we measured accuracy, sensitivity, and specificity over all extracted
lesions and also averaged per TBP image.

Metric Evaluation Type | Over All Lesions (Micro Avg.) | Averaged Over TBPs (Macro Avg.)
---|---|---
Accuracy | 94.16% | 94.23%
Sensitivity | 72.07% | 71.91%
Specificity | 94.70% | 94.95%
Table 2: UD Binary Classification Evaluation. It should be noted that when
calculating average sensitivity over TBP images, only images with at least one
UD lesion were considered.

(a) Before segmentation

(b) After segmentation

Figure 7: Ranking predicted before and after incorporating segmentation
module. Following the same rules as 6

##  6 Discussion

It is sometimes difficult to distinguish pigmented nevi from other skin
conditions, such as lentigos and ephelides (freckles) which often occur on
sun-exposed areas of the upper back, face, back of hands, and forearms. Also,
differentiating small nevi (those less than 2 mm in diameter) from ephelides
is not easy even clinically. Depending on the image resolution, our algorithm
automatically detects all the skin lesions larger than around 1.5mm diameter,
and ranks them in order of similarity to each other, providing an objective
measure of the outlier nature of each lesion. Our algorithm provides an
opportunity to perform quick and unbiased screening for ugly duckling lesions
in TBP images, especially useful for patients with many lesions. Additionally,
our method can work on all types of skin lesions, which means it is not
limited to pigmented lesions.

Our algorithm also benefits much from incorporation of segmentation module.
Before incorporating segmentation module our results were heavily influenced
by noise. Some removed noises were shadow, skin color around a lesion, and
proximity to the body edge. With the help of the segmentation module, our
model made much better predictions as it became more robust to noise and
avoided many false positives. Fig 7 shows how our ranking was improved with
the help of segmentation module. It is important to mention that segmentation
module also resolves the need for additional normalization methods to fix the
lighting. Although normalization methods are mostly helpful, but sometimes
they can drastically decrease the quality of images in some edge-cases.

Limitations of our method include difficulties in detecting UD lesions when
the skin around the lesion has useful information for making the decision
(e.g. scars after excision), as we are not able to incorporate additional
domain knowledge in our decision making. We also remove hairs during detection
and segmentation, which means a loss of a clue such as a hair growing out of a
lesion, which indicates it might have been a congenital nevus. Since our model
does not have domain knowledge, it can not leverage these facts to make better
decisions.

##  7 Conclusion and Future Work

We have developed an algorithm to quickly identify outlier lesions in total
body photography images, at even the relatively low resolution obtained by
imaging backs with mobile smartphone cameras. Higher resolution images of
lesions on arms, etc. can be analysed with even more accuracy. This opens up
the opportunity to use this as a screening aid for patients to submit their
photographs using teledermatology.

Future work will involve further evaluation of the method against a larger
image set of expertly labelled images. We may also use different distance
measures other than the L2 to determine outliers. Also we hope the result of
this UD outlier detection project can improve classification of skin lesions
extracted from TBP images. We hope introduction of our pipeline opens the path
for future research in the domain of TBP images.

## References

  * [1] Naheed R Abbasi, Helen M Shaw, Darrell S Rigel, Robert J Friedman, William H McCarthy, Iman Osman, Alfred W Kopf, and David Polsky.  Early diagnosis of cutaneous melanoma: revisiting the abcd criteria.  JAMA, 292(22):2771–2776, 2004.
  * [2] Charu C. Aggarwal.  Outlier Analysis.  Springer International Publishing, 2 edition, 2017.
  * [3] Jinwon An and Sungzoon Cho.  Variational autoencoder based anomaly detection using reconstruction probability.  Special Lecture on IE, 2(1):1--18, 2015.
  * [4] Giuseppe Argenziano, Gabriella Fabbrocini, Paolo Carli, Vincenzo De Giorgi, Elena Sammarco, and Mario Delfino.  Epiluminescence microscopy for the diagnosis of doubtful melanocytic skin lesions: comparison of the abcd rule of dermatoscopy and a new 7-point checklist based on pattern analysis.  Archives of Dermatology, 134(12):1563–1570, 1998.
  * [5] G. Betta, Giuseppe DI Leo, Gabriella Fabbrocini, Alfredo Paolillo, and M. Scalvenzi.  Automated application of the 7-point checklist diagnosis method for skin lesions: Estimation of chromatic and shape parameters.  Proceedings of the IEEE Instrumentation and Measurement Technology Conference, page 1818–1822, 2005.
  * [6] Judith S. Birkenfeld, Jason M. Tucker-Schwartz, Luis R. Soenksen, José A. Avilés-Izquierdo, and Berta Marti-Fuster.  Computer-aided classification of suspicious pigmented lesions using wide-field images.  Computer Methods and Programs in Biomedicine, 195:105631, Oct. 2020\.
  * [7] Raghavendra Chalapathy, Aditya Krishna Menon, and Sanjay Chawla.  Robust, Deep and Inductive Anomaly Detection.  In Michelangelo Ceci, Jaakko Hollmén, Ljupčo Todorovski, Celine Vens, and Sašo Džeroski, editors, Machine Learning and Knowledge Discovery in Databases, Lecture Notes in Computer Science, pages 36--51, Cham, 2017. Springer International Publishing.
  * [8] Dermengine.  Dermengine total body photography.  = https://www.dermengine.com/en-ca/total-body-photography/.  Accessed: 2021-03-11.
  * [9] Bogdan Dugonik, Aleksandra Dugonik, Maruška Marovt, and Marjan Golob.  Image Quality Assessment of Digital Image Capturing Devices for Melanoma Detection.  Applied Sciences, 10(8):2876, Apr. 2020.
  * [10] Andre Esteva, Brett Kuprel, Roberto A. Novoa, Justin Ko, Susan M. Swetter, Helen M. Blau, and Sebastian Thrun.  Dermatologist-level classification of skin cancer with deep neural networks.  Nature, 542(7639):115--118, Feb. 2017.
  * [11] Caroline Gaudy-Marqueste, Yanal Wazaefi, Yvane Bruneu, Raoul Triller, Luc Thomas, Giovanni Pellacani, Josep Malvehy, Marie-Françoise Avril, Sandrine Monestier, Marie-Aleth Richard, Bernard Fertil, and Jean-Jacques Grob.  Ugly Duckling Sign as a Major Factor of Efficiency in Melanoma Detection.  JAMA Dermatology, 153(4):279--284, 04 2017.
  * [12] J J Grob and J J Bonerandi.  The ’ugly duckling’ sign: identification of the common characteristics of nevi in an individual as a basis for melanoma screening.  Archives of dermatology, 134(1):103--4, Jan 1998.
  * [13] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  CoRR, abs/1512.03385, 2015.
  * [14] Irina Higgins, Loic Matthey, Xavier Glorot, Arka Pal, Benigno Uria, Charles Blundell, Shakir Mohamed, and Alexander Lerchner.  Early visual concept learning with unsupervised deep learning.  arXiv preprint arXiv:1606.05579, 2016.
  * [15] Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner.  beta-vae: Learning basic visual concepts with a constrained variational framework.  International Journal on Learning Representations, 2016.
  * [16] Muneeb Ilyas, Collin M. Costello, Nan Zhang, and Amit Sharma.  The role of the ugly duckling sign in patient education.  Journal of the American Academy of Dermatology, 77(6):1088--1095, 2017.
  * [17] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In International conference on machine learning, pages 448--456. PMLR, 2015.
  * [18] International Society for Digital Imaging of the Skin ISDIS.  Isdis.  https://isdis.org/total-body-photography/.  Accessed: 2021-03-11.
  * [19] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [20] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [21] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection, 2017.
  * [22] Wei Liu, Dragomir Anguelov, Dumitru Erhan, Christian Szegedy, Scott Reed, Cheng-Yang Fu, and Alexander C. Berg.  Ssd: Single shot multibox detector.  Lecture Notes in Computer Science, page 21–37, 2016.
  * [23] W. Liu, R. Li, M. Zheng, S. Karanam, Z. Wu, B. Bhanu, R. J. Radke, and O. Camps.  Towards Visually Explaining Variational Autoencoders.  In 2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 8639--8648, June 2020.  ISSN: 2575-7075.
  * [24] Yuan Liu, Ayush Jain, Clara Eng, David H. Way, Kang Lee, Peggy Bui, Kimberly Kanada, Guilherme de Oliveira Marinho, Jessica Gallegos, Sara Gabriele, Vishakha Gupta, Nalini Singh, Vivek Natarajan, Rainer Hofmann-Wellenhof, Greg S. Corrado, Lily H. Peng, Dale R. Webster, Dennis Ai, Susan J. Huang, Yun Liu, R. Carter Dunn, and David Coz.  A deep learning system for differential diagnosis of skin diseases.  Nature Medicine, 26(6):900--908, June 2020.
  * [25] Yuchen Lu and Peng Xu.  Anomaly Detection for Skin Disease Images Using Variational Autoencoder.  arXiv:1807.01349 [cs, stat], July 2018.  arXiv: 1807.01349.
  * [26] Franz Nachbar, Wilhelm Stolz, Tanja Merkle, Armand B. Cognetta, Thomas Vogt, Michael Landthaler, Peter Bilek, Otto Braun-Falco, and Gerd Plewig.  The abcd rule of dermatoscopy: high prospective value in the diagnosis of doubtful melanocytic skin lesions.  J of American Academy of Derm., 30(4):551--559, 1994.
  * [27] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024--8035. Curran Associates, Inc., 2019.
  * [28] P. Perera and V. M. Patel.  Learning Deep Features for One-Class Classification.  IEEE Transactions on Image Processing, 28(11):5450--5463, Nov. 2019\.
  * [29] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Nassir Navab, Joachim Hornegger, William M. Wells, and Alejandro F. Frangi, editors, Medical Image Computing and Computer-Assisted Intervention -- MICCAI 2015, pages 234--241, Cham, 2015. Springer International Publishing.
  * [30] Thomas Schlegl, Philipp Seeböck, Sebastian M. Waldstein, Ursula Schmidt-Erfurth, and Georg Langs.  Unsupervised Anomaly Detection with Generative Adversarial Networks to Guide Marker Discovery.  In Marc Niethammer, Martin Styner, Stephen Aylward, Hongtu Zhu, Ipek Oguz, Pew-Thian Yap, and Dinggang Shen, editors, Information Processing in Medical Imaging, Lecture Notes in Computer Science, pages 146--157, Cham, 2017. Springer International Publishing.
  * [31] Daniel I. Schlessinger, Guillaume Chhor, Olivier Gevaert, Susan M. Swetter, Justin Ko, and Roberto A. Novoa.  Artificial intelligence and dermatology: opportunities, challenges, and future directions.  Seminars in Cutaneous Medicine and Surgery, 38(1):E31--37, Mar. 2019\.
  * [32] Alon Scope, Stephen W. Dusza, Allan C. Halpern, Harold Rabinovitz, Ralph P. Braun, Iris Zalaudek, Giuseppe Argenziano, and Ashfaq A. Marghoob.  The “Ugly Duckling” Sign: Agreement Between Observers.  Archives of Dermatology, 144(1):58--64, 01 2008.
  * [33] Alon Scope and Ashfaq A. Marghoob.  The “Ugly Duckling” Sign: The “ugly duckling” sign: An early melanoma recognition tool for clinicians and the public.  Melanoma Letters, 25:1--3, 2007.
  * [34] Luis R. Soenksen, Timothy Kassis, Susan T. Conover, Berta Marti-Fuster, Judith S. Birkenfeld, Jason Tucker-Schwartz, Asif Naseem, Robert R. Stavert, Caroline C. Kim, Maryanne M. Senna, José Avilés-Izquierdo, James J. Collins, Regina Barzilay, and Martha L. Gray.  Using deep learning for dermatologist-level detection of suspicious pigmented skin lesions from wide-field images.  Science Translational Medicine, 13(581):eabb3652, Feb. 2021.
  * [35] Xiaoxiao Sun, Jufeng Yang, Ming Sun, and Kai Wang.  A Benchmark for Automatic Visual Classification of Clinical Skin Disease Images.  In Bastian Leibe, Jiri Matas, Nicu Sebe, and Max Welling, editors, Computer Vision – ECCV 2016, Lecture Notes in Computer Science, pages 206--222, Cham, 2016. Springer International Publishing.
  * [36] Canfield Scientific Imaging Systems.  Canfield.  https://www.canfieldsci.com/imaging-systems/vectra-wb360-imaging-system/.  Accessed: 2021-03-11.
  * [37] Philipp Tschandl, Cliff Rosendahl, and Harald Kittler.  The HAM10000 dataset, a large collection of multi-source dermatoscopic images of common pigmented skin lesions.  Scientific Data, 5(1):180161, Aug. 2018.
  * [38] Lilian Weng.  From autoencoder to beta-vae.  lilianweng.github.io/lil-log, 2018.
  * [39] J. Yang, X. Wu, J. Liang, X. Sun, M.-M. Cheng, P. L. Rosin, and L. Wang.  Self-Paced Balance Learning for Clinical Skin Disease Recognition.  IEEE Transactions on Neural Networks and Learning Systems, 31(8):2832--2846, Aug. 2020.


# Two Remarkably Luminous Galaxy Candidates at $z\approx 10-12$ Revealed by
JWST

Rohan P. Naidu Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA  MIT Kavli Institute for Astrophysics and
Space Research, 77 Massachusetts Ave., Cambridge, MA 02139, USA  Pascal A.
Oesch Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290
Versoix, Switzerland  Cosmic Dawn Center (DAWN), Niels Bohr Institute,
University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark  Pieter
van Dokkum Astronomy Department, Yale University, 52 Hillhouse Ave, New Haven,
CT 06511, USA  Erica J. Nelson Department for Astrophysical and Planetary
Science, University of Colorado, Boulder, CO 80309, USA  Katherine A. Suess
Department of Astronomy and Astrophysics, University of California, Santa
Cruz, 1156 High Street, Santa Cruz, CA 95064 USA  Kavli Institute for Particle
Astrophysics and Cosmology and Department of Physics, Stanford University,
Stanford, CA 94305, USA  Gabriel Brammer Cosmic Dawn Center (DAWN), Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, København N, DK-2200,
Denmark  Katherine E. Whitaker Department of Astronomy, University of
Massachusetts, Amherst, MA 01003, USA  Cosmic Dawn Center (DAWN), Denmark
Garth Illingworth Department of Astronomy and Astrophysics, University of
California, Santa Cruz, CA 95064, USA  Rychard Bouwens Leiden Observatory,
Leiden University, NL-2300 RA Leiden, Netherlands  Sandro Tacchella Kavli
Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK  Cavendish Laboratory, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge, CB3 0HE, UK  Jorryt Matthee Department of Physics, ETH
Zürich, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland  Natalie Allen
Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen,
Jagtvej 128, København N, DK-2200, Denmark  Rachel Bezanson Department of
Physics and Astronomy and PITT PACC, University of Pittsburgh, Pittsburgh, PA
15260, USA  Charlie Conroy Center for Astrophysics $|$ Harvard & Smithsonian,
60 Garden Street, Cambridge, MA 02138, USA  Ivo Labbe Centre for Astrophysics
and Supercomputing, Swinburne University of Technology, Melbourne, VIC 3122,
Australia  Joel Leja Department of Astronomy & Astrophysics, The Pennsylvania
State University, University Park, PA 16802, USA  Institute for Computational
& Data Sciences, The Pennsylvania State University, University Park, PA, USA
Institute for Gravitation and the Cosmos, The Pennsylvania State University,
University Park, PA 16802, USA  Ecaterina Leonova GRAPPA, Anton Pannekoek
Institute for Astronomy and Institute of High-Energy Physics, University of
Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands  Dan Magee
UCO/Lick Observatory, University of California, Santa Cruz, CA, 95064  Sedona
H. Price Max-Planck-Institut für extraterrestrische Physik (MPE),
Giessenbachstr. 1, D-85748 Garching, Germany  David J. Setton Department of
Physics and Astronomy and PITT PACC, University of Pittsburgh, Pittsburgh, PA
15260, USA  Victoria Strait Cosmic Dawn Center (DAWN), Niels Bohr Institute,
University of Copenhagen, Jagtvej 128, København N, DK-2200, Denmark  Mauro
Stefanon Departament d’Astronomia i Astrofìsica, Universitat de València, C.
Dr. Moliner 50, E-46100 Burjassot, València, Spain  Unidad Asociada CSIC
”Grupo de Astrofísica Extragaláctica y Cosmología” (Instituto de Física de
Cantabria - Universitat de València)  Sune Toft Cosmic Dawn Center (DAWN),
Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N,
DK-2200, Denmark  John R. Weaver Department of Astronomy, University of
Massachusetts, Amherst, MA 01003, USA  Andrea Weibel Department of Astronomy,
University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland

###### Abstract

The first few hundred Myrs at $z>10$ mark the last major uncharted epoch in
the history of the Universe, where only a single galaxy (GNz11 at $z\approx
11$) is currently spectroscopically confirmed. Here we present a search for
luminous $z>10$ galaxies with JWST/NIRCam photometry spanning $\approx
1-5\mu$m and covering 49 arcmin2 from the public JWST Early Release Science
programs (CEERS and GLASS). Our most secure candidates are two
$M_{\rm{UV}}\approx-21$ systems: GLASS-z12 and GLASS-z10. These galaxies
display abrupt $\gtrsim 1.8$ mag breaks in their spectral energy
distributions, consistent with complete absorption of flux bluewards of
Lyman-$\alpha$ that is redshifted to $z=12.4^{+0.1}_{-0.3}$ and
$z=10.4^{+0.4}_{-0.5}$. Lower redshift interlopers such as quiescent galaxies
with strong Balmer breaks would be comfortably detected at $>5\sigma$ in
multiple bands where instead we find no flux. From SED modeling we infer that
these galaxies have already built up $\sim 10^{9}$ solar masses in stars over
the $\lesssim 300-400$ Myrs after the Big Bang. The brightness of these
sources enable morphological constraints. Tantalizingly, GLASS-z10 shows a
clearly extended exponential light profile, potentially consistent with a disk
galaxy of $r_{\rm{50}}\approx 0.7$ kpc. These sources, if confirmed, join
GNz11 in defying number density forecasts for luminous galaxies based on
Schechter UV luminosity functions, which require a survey area $>10\times$
larger than we have studied here to find such luminous sources at such high
redshifts. They extend evidence from lower redshifts for little or no
evolution in the bright end of the UV luminosity function into the cosmic dawn
epoch, with implications for just how early these galaxies began forming.
This, in turn, suggests that future deep JWST observations may identify
relatively bright galaxies to much earlier epochs than might have been
anticipated.

High-redshift galaxies (734), Galaxy formation (595), Galaxy evolution (594),
Early universe (435)

††facilities: JWST, HST††software:  IPython (Pérez & Granger, 2007),
matplotlib (Hunter, 2007), numpy (Oliphant, 2015), scipy (Virtanen et al.,
2020), jupyter (Kluyver et al., 2016), Astropy (Astropy Collaboration et al.,
2013, 2018), grizli (v1.5.0; Brammer 2018; Brammer et al. 2022), MIST (Choi et
al., 2017), Prospector (Leja et al., 2017, 2019; Johnson et al., 2021), FSPS
(Conroy et al., 2009, 2010; Conroy & Gunn, 2010a, b; Foreman-Mackey et al.,
2014), EAZY (Brammer et al., 2008), SExtractor (Bertin & Arnouts, 1996),
GALFIT (Peng et al., 2002, 2010) ††thanks: NASA Hubble Fellow

##  1 Introduction

When and how the first galaxies formed remains one of the most intriguing
questions of extragalactic astronomy and observational cosmology (see Dayal &
Ferrara, 2018; Robertson, 2021, for recent reviews). Although deep
observations with the Hubble Space Telescope (HST) have pushed our cosmic
horizon to within the first 400 Myr of the Big Bang, galaxies at $z\gtrsim 12$
cannot be observed with HST due to the limit of its wavelength coverage at 1.6
µm.

With the advent of JWST, we now have an unprecedented view of the Universe at
$\sim 2-5$ µm thanks to the extremely sensitive NIRCam instrument (see, e.g.,
Rieke et al., 2005). The extended wavelength coverage enables the study of
rest-frame optical wavelengths up to $z\sim 10$ and allows for rest-frame UV
selections of galaxies out to much higher redshifts.

Here we present first results from a search for particularly luminous $z>10$
sources across the two JWST Early Release Science deep fields. The most
luminous galaxies are of particular importance. They may trace overdensities
and thus pinpoint where galaxy formation first started in the early Universe
(e.g., Leonova et al., 2021; Endsley et al., 2021; Larson et al., 2022).
Furthermore, they provide the most stringent constraints on early galaxy
build-up and promise rich scientific returns.

One particular example of this is provided by GN-z11 (Oesch et al., 2016) that
was detected with HST. Its discovery in the two CANDELS/GOODS fields that
cover a search volume of only $\sim 10^{6}$ Mpc3 was initially quite
surprising. Theoretical and empirical models of early galaxy formation
predicted that a 10-100$\times$ larger survey would have been required to find
one such bright galaxy at $z=11$ (e.g. Waters et al., 2016; Mutch et al.,
2016). This highlights the potential of the brightest galaxies at the cosmic
frontier to set unique constraints on the physics of galaxy formation (see
also Behroozi & Silk, 2018). In particular, the number density of such bright
sources, i.e. the bright end cutoff of the UV luminosity function, is a very
powerful tool to test the efficiency of star-formation and potential feedback
mechanisms in the very early Universe (Bowler et al., 2014; Tacchella et al.,
2018; Bowler et al., 2020).

These results have been extended over the last few years, and evidence is
emerging for a differential evolution of the galaxy population during the
reionization epoch at $z>6$. While the number densities of fainter galaxies
continue to decline with redshift, the most UV-luminous sources seem to be in
place rather early (e.g., Stefanon et al., 2019; Bowler et al., 2020;
Morishita et al., 2020; Harikane et al., 2022; Bagley et al., 2022).
Furthermore, several authors found evidence for pronounced Balmer breaks in
bright $z\sim 8-10$ galaxies, which would indicate a very early formation
epoch with intense star-formation (e.g., Hashimoto et al., 2018; Roberts-
Borsani et al., 2020; Laporte et al., 2021). However, others find very young
ages for the average population (Stefanon et al., 2022a, b). These inferences,
at the moment, are also highly sensitive to the prior adopted on the star-
formation history (e.g., Tacchella et al., 2022; Whitler et al., 2022). Timing
the onset of first star-formation in bright galaxies out to $z\sim 10$ is thus
still highly uncertain, and probing the number density of the most luminous
sources at even higher redshifts is the most direct way of addressing this.

The brightest galaxies are also the most amenable to follow-up studies through
spectroscopy. For instance, GN-z11 had its redshift confirmed both through
grism spectroscopy with HST and with emission lines from ground-based Keck
observations (see Oesch et al., 2016; Jiang et al., 2021). With the
combination of NIRCam and NIRSpec, we enter a new era – the rest-frame optical
features of galaxies in the Epoch of Reionization ($z\approx 6-9$, e.g., Mason
et al. 2019) will come fully into view. However, the strongest emission lines
of galaxies at $z>10$ will still remain out of reach. Spectroscopically
confirming such sources will require measuring order of magnitude fainter
lines or strong continuum breaks. Identifying luminous systems at $z>10$ to
facilitate these measurements is therefore a crucial step in fulfilling JWST’s
mission of charting cosmic dawn.

This paper is organized as follows – §2 describes the datasets analyzed in
this work, §3 presents our sample selection, and in section §4 we show our
results. A discussion on the implications follows in §5, before we conclude
with a summary and an outlook in §6.

Magnitudes are in the AB system (e.g., Oke & Gunn, 1983). For summary
statistics we report medians along with 16th and 84th percentiles. We adopt a
Planck Collaboration et al. (2015) cosmology.

Figure 1: Summary of photometry and redshift solution for GL-z12. Top:
4.5”$\times$4.5” images spanning $\approx 0.9-4.5\mu$m centered on the
$z\approx 12$ candidate highlighted with white crosshairs. The source is well-
detected ($>20\sigma$) in F200W and all redder bands, and abruptly drops out
in the bluer filters. Bottom left: Photometry for the source is shown in
purple, with upper limits for non-detections plotted at the 1-$\sigma$ level.
The best-fit spectral energy distribution (SED) template from EAZY is shown in
dark orange – a Lyman-break galaxy (LBG) at $z=12.2$. The best-fit SED from
EAZY constrained to lie at $z<6$ is plotted in silver, which corresponds to a
quiescent galaxy at $z\approx 3.5$ whose Balmer break produces a drop-off
across F200W and F150W. However, such a quiescent galaxy is predicted to be
detected ($>5\sigma$) in bluer bands, and is at odds with the dramatic $>1.8$
mag break observed. Bottom right: Probability distributions for the source
redshift derived using EAZY (solid orange) and Prospector (dashed orange). We
adopt a flat prior across the redshift range depicted ($z=0-20$). The derived
distributions are in excellent agreement and suggest a redshift of $z\approx
12$, with negligible (EAZY) or no (Prospector) support for solutions at $z<10$

.

Figure 2: Summary of photometry and redshift solution for GL-z10, similar to
Figure 1. Top: GL-z10 is well-detected in all but the two bluest bands. Bottom
left: The best-fit low-$z$ solution (quiescent galaxy at $z\approx 2.5$) is
disfavored by the F115W image, where a $>5\sigma$ detection is expected. In
addition to the JWST data (dark purple), we measure HST photometry (light
purple) for this source from data acquired by the BUFFALO program (Steinhardt
et al., 2020). The HST data are fully consistent with the JWST data as well as
the best-fit SED. Bottom right: The EAZY and Prospector posteriors agree on a
$z\approx 10$ galaxy. Figure 3: Absolute UV magnitude vs. Redshift for a
representative sample of known galaxies in the first billion years of the
Universe. Galaxies with photometric redshifts, sourced from Bouwens et al.
(2022), are shown as points, and those with spectroscopic redshifts compiled
from the literature as squares. The candidates presented in this work are
depicted as purple stars, and populate a hitherto unoccupied region of
parameter space. The brightness of these sources present a unique opportunity
to efficiently extend the spectroscopic frontier to the first few hundred Myrs
after the Big Bang.

##  2 Data

###  2.1 Early Extragalactic JWST Observations

Our analysis is based on some of the first JWST/NIRCam datasets that have been
observed and released over extragalactic fields. In particular, we analyze the
two Early Release Science programs GLASS and CEERS.

The first is the NIRCam parallel field of the ERS program GLASS (PID: 1324,
Treu et al. 2022). This program obtained a single NIRCam pointing in seven
wide filters F090W, F115W, F150W, F200W, F277W, F356W, and F444W, observed for
3.3, 3.3, 1.7, 1.5, 1.5, 1.7, and 6.6 hrs, respectively.

We also analyzed the first four NIRCam pointings from the ERS program CEERS
(PID: 1345, Finkelstein et al., in prep.) that have been observed to date. The
CEERS images include F115W+F277W, F115W+F356W, F150W+F410M, and F200W+F444W
short- (SW) and long-wavelength (LW) exposure sets, for a typical integration
time of 0.8 hr per filter, except for F115W that obtained double this
integration time.

The combined area of these five NIRCam fields used in our analysis amounts to
49 arcmin2 ($\approx 10$ arcmin2 in GLASS, and $\approx 40$ arcmin2 in CEERS)
reaching an unprecedented 5$\sigma$ depth that ranges between 28.6 and 29.6 AB
mag at 4 $\micron$ as measured in 0$\farcs$32 diameter apertures (see Table
1).

###  2.2 Data Reduction

The NIRCam images were reduced with the standard JWST pipeline up to stage 2
using the reference files from jwst_0942.pmap, as well as our own sky flats.
Additional, chip-dependent zeropoint offsets were applied based on our own
reductions of the flux standard star J1743045 as well as observations of the
Large Magellanic Cloud. For some filters, the individual corrections in
different modules of the NIRCam detector can reach up to 30% 111for more
information see: https://github.com/gbrammer/grizli/pull/107. The resulting
images were then processed with the
grizli222https://github.com/gbrammer/grizli/ pipeline for proper alignment
onto a common WCS that was matched to the Gaia DR3 catalogs. To do this,
grizli re-computes the traditional SIP distortion headers that were common for
HST data. This allows us to use drizzlepac as with HST images to combine the
individual frames and produce distortion-corrected mosaics. Additionally,
grizli mitigates the 1/f noise and masks ‘snowballs’ that are most prominent
in short-wavelength filters.

The pipeline was used to drizzle images at 20mas / 40mas pixels for the SW and
LW data, respectively. In the following, our analysis is based on pixel-
matched 40mas images, however. For more information on the grizli processing
see Brammer et al., (in prep).

###  2.3 Multi-Wavelength Catalogs

After processing the new NIRCam images, we produced photometric catalogs in
all fields using the SExtractor software (Bertin & Arnouts, 1996). Sources
were detected in dual mode with two different detection images: F200W or a
weighted combination of all the LW filters. Fluxes were measured in small
circular apertures of 0$\farcs$32 diameter and were corrected to total using
the AUTO flux measurement from the detection image. An additional correction
of typically a few percent only was applied for remaining flux outside of this
aperture based on the predicted encircled energy for the JWST point-spread
functions. Flux uncertainties were estimated based on sigma-clipped histograms
of circular apertures placed throughout the images in random sky positions.
These were then used to rescale the drizzled rms maps. Thus, our uncertainties
are as close to the data as possible. The 5$\sigma$ depths derived in this way
are listed in Table 1.

###  2.4 Quality Control

Given that JWST is a completely new facility for which calibration is still
ongoing, it is important to test the resulting images for any issues. Indeed,
the NIRCam data revealed several features that are not accounted for in the
standard pipeline. In particular, the SW data suffer from significant
scattered light, if there are bright stars in the vicinity. This is
particularly pronounced in the GLASS parallel field where a 10th mag star just
outside of the field seems to cause artificial images across the field. This
is particularly pronounced in the F090W filter in the B4 detector. However,
also other detectors and filters seem to be affected, albeit to a lower
extent. While this issue could be overcome with improved processing, or with
additional observations at a different roll angle, this does not significantly
affect the current analysis. Since we are searching for very high redshift
galaxies that disappear at shorter wavelengths, these data issues only
introduce some level of incompleteness. Most importantly, the areas of the two
candidate sources presented later in this paper are not affected.

Where available, we compared our JWST photometry in the shorter wavelength
filters with existing HST data to check both for issues with residual
distortion or with magnitude zeropoint offsets. In particular, for the GLASS
parallel field, we used HST images made available by Kokorev et al. (2022)
from the BUFFALO survey (Steinhardt et al., 2020). However, only a small
portion of the field is covered by both HST and JWST. For the CEERS data, we
make use of re-reductions from the imaging taken by the CANDELS survey (e.g.,
Koekemoer et al., 2012). No significant offsets or issues were detected.
Zeropoint corrections remained small ($<15\%$; see next section).

Table 1: 5$\sigma$ Depth of JWST Data in this Analysis

Band | GLASS-ERS | CEERS-ERS
---|---|---
$F$090$W$ | 29.0 |  $\cdots$
$F$115$W$ | 29.0 | 28.8
$F$150$W$ | 28.9 | 28.6
$F$200$W$ | 29.1 | 28.8
$F$277$W$ | 29.4 | 28.9
$F$356$W$ | 29.4 | 29.1
$F$410$M$ |  $\cdots$ | 28.4
$F$444$W$ | 29.6 | 28.7

Note. — Measured in 0$\farcs$32 diameter circular apertures.

Table 2: Photometry in units of nJy

Band | GL-z10 | GL-z12
---|---|---
$F$090$W$ | 4$\pm$4 | 5$\pm$4
$F$115$W$ | 1$\pm$3 | 8$\pm$5
$F$150$W$ | 52$\pm$3 | 13$\pm$3
$F$200$W$ | 67$\pm$2 | 71$\pm$2
$F$277$W$ | 68$\pm$2 | 59$\pm$2
$F$356$W$ | 71$\pm$2 | 51$\pm$2
$F$444$W$ | 98$\pm$2 | 57$\pm$2

Note. — We set an error floor of $10\%$ on our measured fluxes for EAZY and
Prospector fits to account for systematic uncertainty not reflected in the
errors stated above.

Table 3: Summary of properties.

| GL-z10 | GL-z12
---|---|---
R.A. | +0:14:02.86 | +0:13:59.76
Dec. |  $-$30:22:18.7 |  $-$30:19:29.1
Redshift $z_{\rm{Prospector}}$ | $10.4^{+0.4}_{-0.5}$ | $12.4^{+0.1}_{-0.3}$
Redshift $z_{\rm{EAZY}}$ | $10.4^{+0.2}_{-0.7}$ | $12.2^{+0.1}_{-0.2}$
Stellar Mass $\log$($M_{\rm{\star}}/M_{\rm{\odot}}$) | $9.6^{+0.2}_{-0.4}$ | $9.1^{+0.3}_{-0.4}$
UV Luminosity ($M_{\rm{UV}}$) | $-20.7^{+0.2}_{-0.2}$ | $-21.0^{+0.1}_{-0.1}$
UV Slope ($\beta$; $f_{\rm{\lambda}}\propto\lambda^{\beta}$) | $-1.9^{+0.2}_{-0.1}$ | $-2.3^{+0.2}_{-0.2}$
Dust Attenuation ($A_{\rm{5500\AA}}$) | $0.3^{+0.4}_{-0.2}$ | $0.1^{+0.2}_{-0.1}$
Dust Attenuation ($A_{\rm{1500\AA}}$) | $0.8^{+0.7}_{-0.5}$ | $0.3^{+0.4}_{-0.2}$
Age ($t_{\rm{50}}$/Myr) | $163^{+20}_{-133}$ | $111^{+26}_{-83}$
SFR${}_{\rm{50\ Myr}}$ ($M_{\rm{\odot}}$/yr) | $10^{+17}_{-5}$ | $6^{+5}_{-2}$
$r_{\mathrm{eff}}$ [kpc] | 0.7 | 0.5
Sersic Index $n$ | 0.8 | 1.0

Note. — SED fitting assumes a continuity prior on the star-formation history
and a Chabrier (2003) IMF.

##  3 Sample Selection & Methods

Photometric redshifts form the basis of our search for bright $z>10$ galaxies.
We fit redshifts using EAZY (Brammer et al., 2008) adopting the
“tweak_fsps_QSF_12_v3” template set derived from FSPS (Conroy et al., 2009,
2010; Conroy & Gunn, 2010a). The allowed range is $0.1<z<20$ adopting a flat
luminosity prior, after applying modest ($<15\%$) zero-point corrections that
are derived iteratively. Candidates of interest are selected to have best-fit
redshifts $z>10$ along with $>84\%$ of their derived redshift probability
distribution function, $p(z)$, lying at $z>10$.

We further require $>10\sigma$ detections in both F356W and F444W. These bands
sample the rest-frame UV at $z>10$ and are critical in establishing the flux
levels with respect to which we seek strong Lyman breaks.

We inspect images of every candidate source for data quality issues (e.g.,
contamination from neighbors, diffraction spikes, location on the edge of the
detector). In tandem, we examine plausible low-$z$ solutions for the
candidates by running EAZY constrained to $z<6$ (primarily dusty, quiescent
galaxies with Balmer breaks masquerading as Lyman breaks). We inspect low-$z$
solutions with the understanding that the errors on fluxes in the dropout band
may be underestimated in some cases (in e.g., image areas with residual
striping), which can have an important impact on the $p(z)$. For instance, we
find one bright source in CEERS with a confident $z_{\rm{EAZY}}\approx 17$ and
$z_{\rm{Prospector}}\approx 17$ whose $z>10$ solution assumes a secure non-
detection in F150W – unfortunately, the source falls in a low-SNR region of
the F150W image and it is difficult to judge the reality of the non-detection.

We find 5 plausible $z>10$ candidates that survive all our conservative checks
– 3 in CEERS, and 2 in GLASS. Of these, the two identified in GLASS – GL-z10
($z_{\rm{EAZY}}=10.4^{+0.2}_{-0.7}$) and GL-z12
($z_{\rm{EAZY}}=12.4^{+0.2}_{-0.2}$) – stand out as being particularly
luminous and secure. No other objects when fit with the Prospector SED-fitting
code (see following Section) have a $p(z)$ with all modes contained at $z>10$.
Further, the GLASS candidates are among the most luminous found – GL-z10, in
particular, is by far the brightest of all sources (by $>2\times$ in F444W).
For the rest of this work, we focus on these two particularly luminous
candidates and we defer the rest of the sources to future papers that present
an analysis of the full $z>10$ galaxy population in these fields.

##  4 Results

###  4.1 Two Luminous $z>10$ Galaxy Candidates

We confirm the photometric redshifts for the two GLASS candidates and derive
stellar population properties using the Prospector SED fitting code (Leja et
al., 2017, 2019; Johnson et al., 2021). The SED parameter space explored by
Prospector is more expansive than EAZY’s linear template combinations, and
therefore it acts as an important check on our derived redshifts. We use FSPS
(Conroy et al., 2009, 2010; Conroy & Gunn, 2010a) with the MIST stellar models
(Choi et al., 2017). We adopt the 19-parameter physical model and parameter
choices described in Tacchella et al. (2022) that fits for the redshift,
stellar and gas-phase metallicities, stellar mass, star-formation history,
dust properties, AGN emission, and scaling of the IGM attenuation curve. We
make slight modifications to their setup – in particular we explore a broader
redshift range of $z=0.1-20$ and keep two bins fixed at lookback times of
$0-5$ Myr and $5-10$ Myrs in the star-formation history following Whitler et
al. (2022) to capture recent bursts that may be powering extreme nebular
emission expected to occur generically at the redshifts of interest (e.g.,
Labbé et al., 2013; De Barros et al., 2019; Endsley et al., 2019). We adopt a
“continuity” prior on the star-formation history, which limits the amount of
variance across consecutive time-bins resulting in smooth histories (Leja et
al., 2019; Tacchella et al., 2022). For further details we direct readers to
Table 1 and §3.4 of Tacchella et al. (2022).

The redshift fits from Prospector are in excellent agreement with EAZY – we
find $z_{\rm{Prospector}}=10.4^{+0.4}_{-0.5}$ for GL-z10 and
$z_{\rm{Prospector}}=12.4^{+0.1}_{-0.3}$ for GL-z12. The photometry and
redshift inference for these sources are summarized in Figures 1 and 2, with
fluxes listed in Table 2. We confirm that no significant data quality issues
affects the $z>10$ candidacy of the sources in the imaging. We derive
$p(z>10)\approx 100\%$ for GL-z12 and $p(z>9.4)\approx 100\%$ for GL-z10, with
their dramatic $\gtrsim 1.8$ mag breaks explained by total absorption of
photons bluewards of Lyman-$\alpha$ by neutral Hydrogen in the intergalactic
medium.

Both galaxies are detected at very high significance in all filters longward
of their break, by virtue of our selection. While they appear very luminous in
the JWST data, these sources have UV absolute magnitudes
($M_{\mathrm{UV}}\approx-21$) that correspond to $L_{\mathrm{UV}}^{*}$ at
$z\sim 8-10$ (see, e.g., Bouwens et al., 2021). This also makes them 1 mag
fainter than GN-z11 and even 2.5 mag fainter than the possible $z\sim 13$
galaxy candidate HD1 (Harikane et al., 2022). Hence, these sources are not
really extreme outliers (see also Fig. 3). Nevertheless, it is interesting
that the first few images with JWST already reveal two such bright sources. We
will discuss their implications on the UV LF in a later section.

###  4.2 Possible Lower Redshift Contamination

The non-detections of both sources in deep, shorter wavelength images
essentially rules out a lower redshift solution. Nevertheless, it is
interesting to explore the nature of possible contaminants. We thus rerun our
photometric redshift codes and force them to find lower redshift fits. The
best $z<6$ solutions in our low-$z$ EAZY runs for these sources are $\approx
10^{8-9}M_{\rm{\odot}}$ quiescent galaxies at $z\approx 3.4$ ($z\approx 2.5$)
with Balmer breaks straddling the dropout filter (silver SEDs in Figs. 1, 2).
However note that Balmer breaks, even in the most pathological cases (e.g.,
4000Å falls just redward of the dropout filter in a super-solar metallicity
galaxy as old as the age of the Universe at $z\approx 2-3.5$), can only
produce drops of $\lesssim 1.5$ mag (assuming no attenuation). The best-fit
low-$z$ solutions have $A_{\rm{V}}\approx 0.1$ – stronger attenuation that
deepens the break is disfavored by the blue continuum slope at wavelengths
longer than the break. In other words, the best-fit low-$z$ solutions predict
$>5\sigma$ detections in bands where we find no flux, and continuum slopes
redder than we observe.

In order to allow for possible systematic effects in the new JWST data, we
perform further testing. We re-fit redshifts to multiple versions of
photometry for these sources – e.g., by adding PSF corrections using WebbPSF,
by increasing the error floor on the photometry, by extracting photometry
using different apertures and detection bands. The only test that produces
viable low-$z$ solutions is when we set a 10 nJy error-floor on all photometry
– this is roughly the level in the SW filters at which the strongest Balmer
breaks at $z\approx 2-3$ can no longer be ruled out (see open silver squares
in bottom-left panels of Figures 1 and 2). This test is a vivid demonstration
of why the sensitivity of JWST/NIRCam is required to identify objects like
GL-z10 and GL-z12 with confidence.

###  4.3 Physical Properties – Billion $M_{\rm{\odot}}$ Galaxies within
$\approx 400$ Myrs of the Big Bang

While the discovery of GN-z11 has already demonstrated that the formation of
billion solar mass galaxies was well underway at $\sim$400 Myr after the Big
Bang, the discovery of these two new sources allows us to derive further
constraints on the physical properties of galaxies at this very early epoch of
the Universe. The Prospector results are summarized in Table 3. In order to
efficiently sample the redshift range of interest, we assume a tighter
redshift prior (a Gaussian centred on the EAZY $p(z)$ with width set to the
84th \- 16th percentile) than in our previous runs when fitting for the
redshift.

The stellar mass for both objects is constrained to be $\approx
10^{9}M_{\rm{\odot}}$, comparable to GNz11 (Oesch et al., 2016; Johnson et
al., 2021; Tacchella et al., 2022). We have verified the stellar mass is
stable to changes in the star-formation history prior by also testing the
“bursty” prior from Tacchella et al. (2022) which allows more rapid
fluctuations in the SFH from time-bin to time-bin than the fiducial model. The
star-formation rates averaged over the last 50 Myrs (SFR50) are typical for
galaxies of comparable mass at $z\approx 7-10$ (e.g., Stefanon et al., 2022a).
The SEDs are consistent with negligible dust attenuation and have blue UV
slopes, $\beta\lesssim-2$. We note that all these derived properties from the
SED fits are collectively consistent with a $z>10$ interpretation for these
galaxies.

Figure 4: Results of the GALFIT morphology analysis for our two sources
(GL-z10 top and GL-z12 bottom). The different columns from left to right
correspond to the original data (in the F444W filter), the model, and the
residual. The sizes and Sersic profiles of both sources are well constrained.
GL-z10 shows some clear extension, consistent with a disk galaxy of 0.7 kpc at
$z\sim 11$. GL-z12 appears quite compact with an estimated size of 0.5 kpc.

###  4.4 Galaxy-galaxy Lensing

Galaxy-galaxy lensing may be particularly important at the redshift frontier
where flux-limited surveys may be preferentially sampling magnified sources
(e.g., Wyithe et al., 2011). Here we make a simple estimate of how lensed our
sources are by assuming their neighbors are singular isothermal spheres (e.g.,
Fort & Mellier, 1994; Schneider et al., 2006; Treu, 2010) following e.g.,
McLure et al. (2006); Oesch et al. (2014); Matthee et al. (2017). For this
estimate, galaxy redshifts and stellar masses are based on our EAZY fits.
Velocity dispersions that trace the underlying dark matter halos are inferred
from the stellar mass by extrapolating the empirical scaling relation in Zahid
et al. (2016) that is fit to $z<0.7$ quiescent galaxies that span
$M_{\rm{\star}}\approx 10^{9}-10^{12}M_{\rm{\odot}}$. We choose a local
relation to cover the low masses relevant to the most massive neighbors at
$<10^{\prime\prime}$ that are likely to produce significant magnification,
while noting that the redshift evolution of such relations at least for
$M_{\rm{\star}}\gtrsim 10^{11}M_{\rm{\odot}}$ is expected to be gradual –
e.g., $\approx 20\%$ higher dispersion at fixed stellar mass at $z\approx 2$,
(e.g., Mason et al., 2015b).

For both GL-z10 and GL-z12 we find negligible lensing ($\mu<1.1$) from all
foreground sources within 10$\farcs$. GL-z12 has two relatively massive
$M_{\rm{\star}}\approx 10^{9}M_{\rm{\odot}}$ neighbors apparent in the bottom-
left quadrant of the stamps in Fig. 1. Even for these two nearby neighbors
($z\approx 2$ at a separation of 0$\farcs$8, and $z\approx 3$ at 2$\farcs$0)
the lensing is modest. We further note that GL-z13 has a compact morphology
(§4.5) that does not show elongation along any particular direction that would
hint at strong magnification. Based on these considerations we conclude that
the observed luminosities of GL-z10 and GL-z12 are likely to be their
intrinsic luminosities.

###  4.5 The Sizes of Luminous $z\approx 10-12$ Galaxies

We fit the sizes of both candidates in the F444W imaging
($\lambda_{\mathrm{rest}}\sim 3500\mathrm{\AA}$) using GALFIT (Peng et al.,
2010). We create 100-pixel cutouts around each galaxy, then use photutils and
astropy to create a segmentation map to identify nearby galaxies. We
simultaneously fit any sources with magnitudes (estimated from the
segmentation map) up to 2.5 magnitudes fainter than the target galaxy that
have centers within 3″of the galaxy; we mask fainter or more distant galaxies.
In our fits, we constrain the center of the target galaxy to be within 10
pixels (0$\farcs$4) of the input value, the Sersic index $n$ to be between
0.01 and 8, the magnitude to be between 0 and 45, and the half-light radius
$r_{e}$ to be between 0.3 and 200 pixels (0$\farcs$012 - 8$\farcs$0). We
calculate and subtract off a scalar sky background correction from each
cutout, estimated from the masked, sigma-clipped cutout, then fix the sky
background component in GALFIT to zero. We use a theoretical PSF model
generated from WebbPSF at our 0$\farcs$04 pixel scale; we oversample the PSF
by a factor of 9 in order to minimize artifacts as we rotate the PSF to the
GLASS observation angle calculated from the APT file, then convolve with a 9x9
pixel square kernel and downsample to the mosaic resolution.

We find reliable Sersic fits for both galaxies, with half-light radii of 0.5
and 0.7 kpc, respectively, and disk-like profiles ($n=1$ and $n=0.8$,
respectively). The models are shown in Fig. 4, and the size and Sersic profile
estimates are listed in Table 3.

The resulting sizes of 0.5 and 0.7 kpc are typical for luminous $L^{*}$
galaxies at $z\sim 6-9$, where measurements have been possible to date (e.g.,
Holwerda et al., 2015; Shibuya et al., 2015; Bowler et al., 2017; Kawamata et
al., 2018; Yang et al., 2022). They are also consistent with expectations from
simulations for $z>9$ galaxies (e.g., Roper et al., 2022; Marshall et al.,
2022). However, at $z\sim 7$, the most luminous sources often break up in
multiple clumps (Bowler et al., 2017). This is not the case for these two
sources, at least down to the resolution limit of order 500 pc for the F444W
bandpass. Interestingly, GL-z10 even shows tantalizing evidence for being an
ordered disk galaxy at $z\sim 10$, based on the exponential light profile and
elliptical morphology. If we interpret GL-z10’s projected axis ratio of 0.65
using a sample of randomly oriented axisymmetric oblate rotators (following
e.g., Holden et al. 2012, Chang et al. 2013, van der Wel et al. 2014) and
adopt $c/a\leq 0.4$ as a threshold for disks, we find that the observed axis
ratio implies P(disk) $\sim 0.5$. Our analysis shows the unparalleled power of
JWST to provide accurate profile measurements of early Universe galaxies.

##  5 Discussion

###  5.1 Caveats

The key caveat, as well as the key animating spirit of this work, is that
these data are among the first deep extragalactic fields collected by a new
Great Observatory. Systematic uncertainties (e.g., zero-point corrections,
treatment of artefacts) can still be significant. We have tested for zero-
point offsets by comparing HST and JWST photometry for brighter sources where
possible and have not found any major issues at the $\lesssim 10\%$ level.
Nevertheless, we have attempted to account for remaining uncertainties with
conservative choices – e.g., a $10\%$ error floor on all fluxes and focusing
on bright galaxies whose $>2$ mag breaks are robust to even major
uncertainties.

A next caveat applies to the SED models underpinning the stellar population
parameters and photo-$z$ fitting. Important details about the nebular emission
and nature of massive stars at low metallicities, which dominate the light in
these few hundred Myr old star-forming systems, remain unconstrained (e.g.,
Stanway et al., 2020). These uncertainties directly translate to the
parameters we recover from SED fitting and the sources for which we are able
to fit high quality redshifts. Fortunately, for our redshift range of
interest, from the perspective of continuum fitting for the stellar population
analyses, extreme nebular emission from strong rest-frame optical lines is
shifted out of all NIRCam bands and the most important feature for the
photo-$z$s in the bright galaxies we study is the Lyman break. Further, there
exists a range of plausible, but hitherto unconstrained physical ingredients
that are unaccounted for in our models (e.g., primordial AGN, top-heavy IMFs,
super-luminous Pop III stars; Windhorst et al. 2018; Pacucci et al. 2022;
Steinhardt et al. 2022). Some of these ingredients may potentially produce
large UV luminosities in the absence of substantial stellar mass.
Spectroscopic follow-up is therefore essential to confirm the nature of these
sources.

Figure 5: Constraints on the bright end of the UV LF at $z\sim 10-13$. The
current JWST data allow us to derive a first estimate of the number density of
galaxies with $M_{\mathrm{UV}}\sim-21$ at these redshifts (purple star). While
this estimate lies a factor $\sim 10\times$ above the extrapolation of
Schechter function constraints to $z=11.5$ from Bouwens et al. (2021, dashed
black line), they are in very good agreement with extrapolated double-power
law LFs from Bowler et al. (2020, black solid line). In fact, GN-z11 (orange
star) is also consistent with the double-power law LF. Other LF estimates and
upper limits at $z\sim 10$ are shown as open symbols (see legend for
references). Simulated predictions are shown from the UniverseMachine models
at $z\sim 11$ and $z\sim 12$ (dotted lines) and from the Delphi model at
$z\sim 11$ (darker gray shaded region).

###  5.2 Implication: The Number Density of Luminous Galaxies in the Early
Universe

These two galaxies enable a first estimate of the number densities of
relatively luminous sources at $z=10-13$. Given the large uncertainties in
redshift, we conservatively estimate a selection volume across this full
redshift range, i.e., $\Delta z=3$. Given that the depth of all the fields we
studied is far deeper than the bright magnitude of these two sources, these
galaxies would have been identified over the full area, without foreground
galaxies. This amounts to a volume of $2.2\times 10^{5}$ Mpc3. The detection
of two galaxies with $M_{\mathrm{UV}}=-21$ thus results in an estimated UV LF
point of $\log\phi[\mathrm{Mpc}^{-3}\mathrm{mag}^{-1}]=-5.05^{+0.37}_{-0.45}$.
In §4.4 we conclude that this UV LF point likely reflects the intrinsic
luminosity of the sources, and is unaffected by lensing.

In Figure 5, we compare this estimate with previous UV LF determinations and
with extrapolations from lower redshifts. In particular, we show that an
extrapolation of the Schechter function trends estimated at $z=3-10$ results
in an LF at $M_{\mathrm{UV}}=-21$ that is a factor $>$10$\times$ lower than
our estimate. Interestingly, however, when extrapolating the trends in the
double-power law LFs from Bowler et al. (2020) to $z\sim 11.5$ (the mean
redshift of our sources), we find relatively good agreement. In fact, GN-z11
also lies on this extrapolated LF. However, this would indicate very little
evolution in the bright galaxy population at $z>8$. Indeed, our estimate is in
good agreement with previous $z\sim 10$ UV LF determinations and constraints
at the bright end from Oesch et al. (2018); Bouwens et al. (2019); Stefanon et
al. (2019); Morishita et al. (2020); Finkelstein et al. (2022); Bagley et al.
(2022); Leethochawalit et al. (2022).

Finally, we also briefly compare our estimates with simulated LFs from the
UniverseMachine model (Behroozi et al., 2019) and from the Delphi model (Dayal
et al., 2014, 2022). While our estimate is in rough agreement with these
prediction at $z\sim 11$, the model LFs evolve very rapidly at these early
times, such that the $z\sim 12$ LF is already $>30\times$ below our estimate.
This is a general trend of model predictions: a relatively rapid evolution of
the LF at $z>10$, driven by the underlying evolution of the dark matter halo
mass function (see also Oesch et al., 2018; Tacchella et al., 2018; Bouwens et
al., 2021). However, the handful of bright galaxies that have been found at
$z\sim 10-13$ to date appear to oppose this trend.

It is still unclear what the physical reason for this might be. Combined with
the discovery of GN-z11 (Oesch et al., 2016), and taking our SED fits at face
value (see §5.1 for caveats), evidence is mounting that the star-formation
efficiency in the early Universe may be much higher than expected (e.g.,
Tacchella et al. 2013; Mason et al. 2015a; Tacchella et al. 2018) in at least
a few sources, thus resulting in the early appearance of UV-luminous galaxies
with stellar masses as high as $\approx 10^{9}M_{\rm{\odot}}$ already a few
hundred Myrs after the Big Bang. The existence of these massive galaxies at
such early times raises interesting questions about just how early such
galaxies began forming, potentially earlier than current expectations. Wider
area datasets will be required to increase the search volume, for more
reliable constraints on the number densities of luminous sources.

##  6 Summary & Outlook

This paper presented a search for luminous $z>10$ galaxies across the two JWST
Early Release Science programs in extragalactic fields. We find the following
–

  * •

We identify two particularly luminous sources in the GLASS ERS program. These
sources, GL-z10 and GL-z12, have continuum magnitudes of $\sim 27$ at 2
$\micron$ and display dramatic $>1.8$ mag breaks in their SEDs that are best
fit as Lyman breaks occurring at redshifts of $z\approx 11$ and $z\approx 13$
respectively. [Fig. 1, Fig. 2, §4.1]

  * •

SED modeling of these sources shows they have properties (e.g., $\beta$
slopes, specific star-formation rates) expected of $z>10$ galaxies. These
systems are a billion solar mass galaxies, having built up their mass only
$<300-400$ Myrs after the Big Bang. [Table 3, §4.3]

  * •

The brightness of these objects present a unique opportunity for detailed
spectroscopic and morphological follow-up at $z>10$. As a demonstration, we
model the morphology of both galaxies finding that they are well-described by
disk-like profiles with small sizes (half-light radii $\sim 0.6$ kpc). GL-z10,
in particular, shows an extended exponential light profile, that may be
tantalizing evidence for a disk already in place at $z\approx 10$. [Fig. 4,
§4.5]

  * •

These two objects already place novel constraints on galaxy evolution in the
cosmic dawn epoch. They indicate that the discovery of GN-z11 was not simply a
matter of good fortune, but that there is likely a population of UV luminous
sources with very high star-formation efficiencies capable of compiling
$>10^{9}M_{\rm{\odot}}$ at $z>10$. [Fig. 3, §5.2]

  * •

The inferred number-density of $M_{\rm{UV}}\approx-21$ sources from our search
($\log\phi[\mathrm{Mpc}^{-3}\mathrm{mag}^{-1}]=-5.05^{+0.37}_{-0.45}$)
strongly supports a significant deviation from the Schechter UV luminosity
function at the bright end, and is consistent with the double-power law
evolution reported at lower redshifts. The physical mechanisms driving this
departure are yet to be definitively established. These luminous sources
highly conducive to NIRSpec spectroscopy may hold the key. [Fig. 5, §5.2]

If these candidates are confirmed spectroscopically, and indeed two $z\approx
10-12$ candidates lie awaiting discovery in every $\sim$50 arcmin2
extragalactic field, it is clear that JWST will prove highly successful in
pushing the cosmic frontier all the way to the brink of the Big Bang.

We thank the referee for insightful comments that strengthened this
manuscript. We are grateful to the CEERS and GLASS teams for planning these
early release observations. RPN acknowledges funding from JWST programs
GO-1933 and GO-2279. Support for this work was provided by NASA through the
NASA Hubble Fellowship grant HST-HF2-51515.001-A awarded by the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Incorporated, under NASA contract
NAS5-26555. We acknowledge support from: the Swiss National Science Foundation
through project grant 200020_207349 (PAO, AW). The Cosmic Dawn Center (DAWN)
is funded by the Danish National Research Foundation under grant No. 140. RJB
and MS acknowledge support from NWO grant TOP1.16.057. Cloud-based data
processing and file storage for this work is provided by the AWS Cloud Credits
for Research program. This work is based on observations made with the
NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the
Mikulski Archive for Space Telescopes at the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These
observations are associated with programs # 1324 and # 1345.

## References

  * 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
  * Bagley et al. (2022) Bagley, M. B., Finkelstein, S. L., Rojas-Ruiz, S., et al. 2022, arXiv e-prints, arXiv:2205.12980.  https://arxiv.org/abs/2205.12980
  * Behroozi & Silk (2018) Behroozi, P., & Silk, J. 2018, MNRAS, 477, 5382, doi: 10.1093/mnras/sty945
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 1134, doi: 10.1093/mnras/stz1182
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
  * Bouwens et al. (2019) Bouwens, R. J., Stefanon, M., Oesch, P. A., et al. 2019, ApJ, 880, 25, doi: 10.3847/1538-4357/ab24c5
  * Bouwens et al. (2021) Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, arXiv e-prints, arXiv:2102.07775.  https://arxiv.org/abs/2102.07775
  * Bouwens et al. (2022) Bouwens, R. J., Smit, R., Schouws, S., et al. 2022, ApJ, 931, 160, doi: 10.3847/1538-4357/ac5a4a
  * Bowler et al. (2017) Bowler, R. A. A., Dunlop, J. S., McLure, R. J., & McLeod, D. J. 2017, MNRAS, 466, 3612, doi: 10.1093/mnras/stw3296
  * Bowler et al. (2020) Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059, doi: 10.1093/mnras/staa313
  * Bowler et al. (2014) Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2014, MNRAS, 440, 2810, doi: 10.1093/mnras/stu449
  * Brammer (2018) Brammer, G. 2018, Gbrammer/Grizli: Preliminary Release, 0.4.0, Zenodo, Zenodo, doi: 10.5281/zenodo.1146905
  * Brammer et al. (2022) Brammer, G., Strait, V., Matharu, J., & Momcheva, I. 2022, grizli, 1.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.6672538
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Chabrier (2003) Chabrier, G. 2003, Publications of the Astronomical Society of the Pacific, 115, 763, doi: 10.1086/376392
  * Chang et al. (2013) Chang, Y.-Y., van der Wel, A., Rix, H.-W., et al. 2013, ApJ, 773, 149, doi: 10.1088/0004-637X/773/2/149
  * Choi et al. (2017) Choi, J., Conroy, C., & Byler, N. 2017, ApJ, 838, 159, doi: 10.3847/1538-4357/aa679f
  * Conroy & Gunn (2010a) Conroy, C., & Gunn, J. E. 2010a, ApJ, 712, 833, doi: 10.1088/0004-637X/712/2/833
  * Conroy & Gunn (2010b) —. 2010b, FSPS: Flexible Stellar Population Synthesis, Astrophysics Source Code Library, record ascl:1010.043.  http://ascl.net/1010.043
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Conroy et al. (2010) Conroy, C., White, M., & Gunn, J. E. 2010, ApJ, 708, 58, doi: 10.1088/0004-637X/708/1/58
  * Dayal & Ferrara (2018) Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1, doi: 10.1016/j.physrep.2018.10.002
  * Dayal et al. (2014) Dayal, P., Ferrara, A., Dunlop, J. S., & Pacucci, F. 2014, MNRAS, 445, 2545, doi: 10.1093/mnras/stu1848
  * Dayal et al. (2022) Dayal, P., Ferrara, A., Sommovigo, L., et al. 2022, MNRAS, 512, 989, doi: 10.1093/mnras/stac537
  * De Barros et al. (2019) De Barros, S., Oesch, P. A., Labbé, I., et al. 2019, MNRAS, 489, 2355, doi: 10.1093/mnras/stz940
  * Endsley et al. (2019) Endsley, R., Behroozi, P., Stark, D. P., et al. 2019, arXiv e-prints, arXiv:1907.02546.  https://arxiv.org/abs/1907.02546
  * Endsley et al. (2021) Endsley, R., Stark, D. P., Chevallard, J., & Charlot, S. 2021, MNRAS, 500, 5229, doi: 10.1093/mnras/staa3370
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M., Song, M., et al. 2022, ApJ, 928, 52, doi: 10.3847/1538-4357/ac3aed
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Sick, J., & Johnson, B. 2014, python-fsps: Python bindings to FSPS (v0.1.1), doi: 10.5281/zenodo.12157
  * Fort & Mellier (1994) Fort, B., & Mellier, Y. 1994, A&A Rev., 5, 239, doi: 10.1007/BF00877691
  * Harikane et al. (2022) Harikane, Y., Inoue, A. K., Mawatari, K., et al. 2022, ApJ, 929, 1, doi: 10.3847/1538-4357/ac53a9
  * Hashimoto et al. (2018) Hashimoto, T., Laporte, N., Mawatari, K., et al. 2018, Nature, 557, 392, doi: 10.1038/s41586-018-0117-z
  * Holden et al. (2012) Holden, B. P., van der Wel, A., Rix, H.-W., & Franx, M. 2012, ApJ, 749, 96, doi: 10.1088/0004-637X/749/2/96
  * Holwerda et al. (2015) Holwerda, B. W., Bouwens, R., Oesch, P., et al. 2015, ApJ, 808, 6, doi: 10.1088/0004-637X/808/1/6
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jiang et al. (2021) Jiang, L., Kashikawa, N., Wang, S., et al. 2021, Nature Astronomy, 5, 256, doi: 10.1038/s41550-020-01275-y
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22, doi: 10.3847/1538-4365/abef67
  * Kawamata et al. (2018) Kawamata, R., Ishigaki, M., Shimasaku, K., et al. 2018, ApJ, 855, 4, doi: 10.3847/1538-4357/aaa6cf
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Koekemoer et al. (2012) Koekemoer, A. M., Ellis, R. S., McLure, R. J., et al. 2012, ArXiv e-prints, 1212.1448
  * Kokorev et al. (2022) Kokorev, V., Brammer, G., Fujimoto, S., et al. 2022, arXiv e-prints, arXiv:2207.07125.  https://arxiv.org/abs/2207.07125
  * Labbé et al. (2013) Labbé, I., Oesch, P. A., Bouwens, R. J., et al. 2013, ApJ, 777, L19, doi: 10.1088/2041-8205/777/2/L19
  * Laporte et al. (2021) Laporte, N., Meyer, R. A., Ellis, R. S., et al. 2021, MNRAS, 505, 3336, doi: 10.1093/mnras/stab1239
  * Larson et al. (2022) Larson, R. L., Finkelstein, S. L., Hutchison, T. A., et al. 2022, ApJ, 930, 104, doi: 10.3847/1538-4357/ac5dbd
  * Leethochawalit et al. (2022) Leethochawalit, N., Roberts-Borsani, G., Morishita, T., Trenti, M., & Treu, T. 2022, arXiv e-prints, arXiv:2205.15388.  https://arxiv.org/abs/2205.15388
  * Leja et al. (2019) Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, ApJ, 876, 3, doi: 10.3847/1538-4357/ab133c
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Leonova et al. (2021) Leonova, E., Oesch, P. A., Qin, Y., et al. 2021, arXiv e-prints, arXiv:2112.07675.  https://arxiv.org/abs/2112.07675
  * Marshall et al. (2022) Marshall, M. A., Wilkins, S., Di Matteo, T., et al. 2022, MNRAS, 511, 5475, doi: 10.1093/mnras/stac380
  * Mason et al. (2019) Mason, C. A., Naidu, R. P., Tacchella, S., & Leja, J. 2019, MNRAS, 489, 2669, doi: 10.1093/mnras/stz2291
  * Mason et al. (2015a) Mason, C. A., Trenti, M., & Treu, T. 2015a, ApJ, 813, 21, doi: 10.1088/0004-637X/813/1/21
  * Mason et al. (2015b) Mason, C. A., Treu, T., Schmidt, K. B., et al. 2015b, ApJ, 805, 79, doi: 10.1088/0004-637X/805/1/79
  * Matthee et al. (2017) Matthee, J., Sobral, D., Darvish, B., et al. 2017, MNRAS, 472, 772, doi: 10.1093/mnras/stx2061
  * McLure et al. (2006) McLure, R. J., Cirasuolo, M., Dunlop, J. S., et al. 2006, MNRAS, 372, 357, doi: 10.1111/j.1365-2966.2006.10864.x
  * Morishita et al. (2020) Morishita, T., Stiavelli, M., Trenti, M., et al. 2020, ApJ, 904, 50, doi: 10.3847/1538-4357/abba83
  * Mutch et al. (2016) Mutch, S. J., Liu, C., Poole, G. B., et al. 2016, MNRAS, 463, 3556, doi: 10.1093/mnras/stw2187
  * Oesch et al. (2018) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105, doi: 10.3847/1538-4357/aab03f
  * Oesch et al. (2014) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2014, ApJ, 786, 108, doi: 10.1088/0004-637X/786/2/108
  * Oesch et al. (2016) Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129, doi: 10.3847/0004-637X/819/2/129
  * Oke & Gunn (1983) Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, doi: 10.1086/160817
  * Oliphant (2015) Oliphant, T. E. 2015, Guide to NumPy (Continuum Press)
  * Pacucci et al. (2022) Pacucci, F., Dayal, P., Harikane, Y., Inoue, A. K., & Loeb, A. 2022, MNRAS, 514, L6, doi: 10.1093/mnrasl/slac035
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H. 2002, AJ, 124, 266, doi: 10.1086/340952
  * Peng et al. (2010) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097, doi: 10.1088/0004-6256/139/6/2097
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Planck Collaboration et al. (2015) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015, ArXiv e-prints.  https://arxiv.org/abs/1502.01589
  * Rieke et al. (2005) Rieke, M. J., Kelly, D., & Horner, S. 2005, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5904, Cryogenic Optical Systems and Instruments XI, ed. J. B. Heaney & L. G. Burriesci, 1–8, doi: 10.1117/12.615554
  * Roberts-Borsani et al. (2020) Roberts-Borsani, G. W., Ellis, R. S., & Laporte, N. 2020, MNRAS, 497, 3440, doi: 10.1093/mnras/staa2085
  * Robertson (2021) Robertson, B. E. 2021, arXiv e-prints, arXiv:2110.13160.  https://arxiv.org/abs/2110.13160
  * Roper et al. (2022) Roper, W. J., Lovell, C. C., Vijayan, A. P., et al. 2022, MNRAS, 514, 1921, doi: 10.1093/mnras/stac1368
  * Schneider et al. (2006) Schneider, P., Kochanek, C. S., & Wambsganss, J. 2006, in Saas-Fee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro, ed. G. Meylan, P. Jetzer, & P. North
  * Shibuya et al. (2015) Shibuya, T., Ouchi, M., & Harikane, Y. 2015, The Astrophysical Journal Supplement Series, 219, 15, doi: 10.1088/0067-0049/219/2/15
  * Stanway et al. (2020) Stanway, E. R., Chrimes, A. A., Eldridge, J. J., & Stevance, H. F. 2020, MNRAS, 495, 4605, doi: 10.1093/mnras/staa1166
  * Stefanon et al. (2022a) Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2022a, arXiv e-prints, arXiv:2206.13525.  https://arxiv.org/abs/2206.13525
  * Stefanon et al. (2022b) —. 2022b, ApJ, 927, 48, doi: 10.3847/1538-4357/ac3de7
  * Stefanon et al. (2019) Stefanon, M., Labbé, I., Bouwens, R. J., et al. 2019, ApJ, 883, 99, doi: 10.3847/1538-4357/ab3792
  * Steinhardt et al. (2020) Steinhardt, C. L., Jauzac, M., Acebron, A., et al. 2020, ApJS, 247, 64, doi: 10.3847/1538-4365/ab75ed
  * Steinhardt et al. (2022) Steinhardt, C. L., Sneppen, A., Mostafa, B., et al. 2022, ApJ, 931, 58, doi: 10.3847/1538-4357/ac62d6
  * Tacchella et al. (2018) Tacchella, S., Bose, S., Conroy, C., Eisenstein, D. J., & Johnson, B. D. 2018, ApJ, 868, 92, doi: 10.3847/1538-4357/aae8e0
  * Tacchella et al. (2013) Tacchella, S., Trenti, M., & Carollo, C. M. 2013, ApJ, 768, L37, doi: 10.1088/2041-8205/768/2/L37
  * Tacchella et al. (2022) Tacchella, S., Finkelstein, S. L., Bagley, M., et al. 2022, ApJ, 927, 170, doi: 10.3847/1538-4357/ac4cad
  * Treu (2010) Treu, T. 2010, ARA&A, 48, 87, doi: 10.1146/annurev-astro-081309-130924
  * Treu et al. (2022) Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110, doi: 10.3847/1538-4357/ac8158
  * van der Wel et al. (2014) van der Wel, A., Chang, Y.-Y., Bell, E. F., et al. 2014, ApJ, 792, L6, doi: 10.1088/2041-8205/792/1/L6
  * 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
  * Waters et al. (2016) Waters, D., Wilkins, S. M., Di Matteo, T., et al. 2016, MNRAS, 461, L51, doi: 10.1093/mnrasl/slw100
  * Whitler et al. (2022) Whitler, L., Stark, D. P., Endsley, R., et al. 2022, arXiv e-prints, arXiv:2206.05315.  https://arxiv.org/abs/2206.05315
  * Windhorst et al. (2018) Windhorst, R. A., Timmes, F. X., Wyithe, J. S. B., et al. 2018, ApJS, 234, 41, doi: 10.3847/1538-4365/aaa760
  * Wyithe et al. (2011) Wyithe, J. S. B., Yan, H., Windhorst, R. A., & Mao, S. 2011, Nature, 469, 181, doi: 10.1038/nature09619
  * Yang et al. (2022) Yang, L., Leethochawalit, N., Treu, T., et al. 2022, MNRAS, 514, 1148, doi: 10.1093/mnras/stac1236
  * Zahid et al. (2016) Zahid, H. J., Geller, M. J., Fabricant, D. G., & Hwang, H. S. 2016, ApJ, 832, 203, doi: 10.3847/0004-637X/832/2/203

##  Appendix A Comparison with Initial Reduction and Calibration

In Figure 6 we compare the results presented in this paper with fluxes, SEDs,
and $p(z)$ based on the initial NIRCam calibrations (e.g., zero-points, flats,
darks) available with the ERS data release in July 2022. Relative to results
based on these calibrations, the $p(z)$ distribution has shifted slightly
towards lower redshifts for both galaxies. GL-z10 now shows a stronger hint of
rest-optical lines or a Balmer break in its $F$444$W$-$F$356$W$ color, along
with a slightly redder $\beta_{\rm{UV}}$ slope ($-1.9$ vs. $-2.1$) which does
not require as much damping by the IGM to explain the F150W flux. GL-z12 is
now detected at higher significance in $F$150$W$ due to better overall
handling of background features (e.g., wisps) in the SW filters, ruling out
the secondary peak in the Prospector $p(z)$ distribution centered at $z\approx
14$.

Figure 6: Comparison of SEDs (left) and $p(z)$ distributions (right) derived
based on initial NIRCam calibrations available at the time of the ERS data
release (light blue) with results based on updated calibrations used in this
paper (orange). Photometry from the initial calibrations is shown in the left
panels in transparent purple, and from the updated calibrations in solid
purple. The SEDs are broadly similar, with the updated calibrations implying
slightly lower redshifts for both sources. The high redshift nature of these
sources remains secure.