# EPOCHS IV: SED Modelling Assumptions and their impact on the Stellar Mass
Function at $6.5~{}\leq z\leq~{}13.5$ using PEARLS and public JWST
observations

Thomas Harvey Jodrell Bank Centre for Astrophysics, University of Manchester,
Oxford Road, Manchester M13 9PL, UK  Christopher J. Conselice Jodrell Bank
Centre for Astrophysics, University of Manchester, Oxford Road, Manchester M13
9PL, UK  Nathan J. Adams Jodrell Bank Centre for Astrophysics, University of
Manchester, Oxford Road, Manchester M13 9PL, UK  Duncan Austin Jodrell Bank
Centre for Astrophysics, University of Manchester, Oxford Road, Manchester M13
9PL, UK  Ignas Juodžbalis Jodrell Bank Centre for Astrophysics, University of
Manchester, Oxford Road, Manchester M13 9PL, UK  Kavli Institute for
Cosmology, University of Cambridge, Cambridge, Madingley Road, Cambridge, CB3
0HA  James Trussler Jodrell Bank Centre for Astrophysics, University of
Manchester, Oxford Road, Manchester M13 9PL, UK  Qiong Li Jodrell Bank Centre
for Astrophysics, University of Manchester, Oxford Road, Manchester M13 9PL,
UK  Katherine Ormerod Jodrell Bank Centre for Astrophysics, University of
Manchester, Oxford Road, Manchester M13 9PL, UK  Astrophysics Research
Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool, L3
5RF  Leonardo Ferreira Department of Physics & Astronomy, University of
Victoria, Finnerty Road, Victoria, British Columbia, V8P 1A1, Canada  Qiao
Duan  Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford
Road, Manchester M13 9PL, UK  Lewi Westcott  Jodrell Bank Centre for
Astrophysics, University of Manchester, Oxford Road, Manchester M13 9PL, UK
Honor Harris  Jodrell Bank Centre for Astrophysics, University of Manchester,
Oxford Road, Manchester M13 9PL, UK  Rachana Bhatawdekar European Space Agency
(ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n,
28692 Villanueva de la Cañada, Madrid, Spain  Dan Coe AURA for the European
Space Agency (ESA), Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Seth H. Cohen School of Earth and Space Exploration,
Arizona State University, Tempe, AZ 85287-1404  Joseph Caruana Department of
Physics, University of Malta, Msida MSD 2080, Malta; & Institute of Space
Sciences & Astronomy, University of Malta, Msida MSD 2080, Malta  Cheng Cheng
Chinese Academy of Sciences South America Center for Astronomy, National
Astronomical Observatories, CAS, Beijing 100101, People’s Republic of China  2
CAS Key Laboratory of Optical Astronomy, National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100101, People’s Republic of China  Simon
P. Driver International Centre for Radio Astronomy Research (ICRAR) and the
International Space Centre (ISC), The University of Western Australia, M468,
35 Stirling Highway, Crawley, WA 6009, Australia  Brenda Frye University of
Arizona, Department of Astronomy/Steward Observatory, 933 N Cherry Ave,
Tucson, AZ85721  Lukas J. Furtak Physics Department, Ben-Gurion University of
the Negev, P. O. Box 653, Be’er-Sheva, 8410501, Israel  Norman A. Grogin Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Nimish P. Hathi Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Benne W. Holwerda Department of Physics and
Astronomy, University of Louisville, 102 Natural Sciences Building,
Louisville, KY 40292, USA  Rolf A. Jansen School of Earth and Space
Exploration, Arizona State University, Tempe, AZ 85287-1404  Anton M.
Koekemoer Space Telescope Science Institute, 3700 San Martin Drive, Baltimore,
MD 21218, USA  Christopher C. Lovell Institute ofCosmology and Gravitation,
University ofPortsmouth, Burnaby Road, Portsmouth, PO1 3FX, UK  Madeline A.
Marshall National Research Council of Canada, Herzberg Astronomy &
Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7,
Canada; & ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions
(ASTRO 3D), Australia  Mario Nonino INAF-Osservatorio Astronomico di Trieste,
Via Bazzoni 2, I-34124 Trieste, Italy  Ian Smail Centre for Extragalactic
Astronomy, Department of Physics, Durham University, South Road, Durham DH1
3LE, UK  Aswin P. Vijayan Cosmic Dawn Center (DAWN)  Astronomy Centre,
Department of Physics and Astronomy, University of Sussex, Brighton, BN1 9QH,
UK  DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs.
Lyngby, Denmark  Stephen M. Wilkins Astronomy Centre, Department of Physics
and Astronomy, University of Sussex, Brighton, BN1 9QH, UK  Rogier Windhorst
School of Earth and Space Exploration, Arizona State University, Tempe, AZ
85287-1404  Christopher N. A. Willmer Steward Observatory, University of
Arizona, 933 N Cherry Ave, Tucson, AZ, 85721-0009  Haojing Yan Department of
Physics and Astronomy, University of Missouri, Columbia, MO 65211  Adi Zitrin
Department of Physics, Ben-Gurion University of the Negev, P.O. Box 653,
Be’er-Sheva 84105, Israel

###### Abstract

We utilize deep JWST NIRCam observations for the first direct constraints on
the Galaxy Stellar Mass Function (GSMF) at $z>10$. Our EPOCHS v1 sample
includes 1120 galaxy candidates at $6.5<z<13.5$ taken from a consistent
reduction and analysis of publicly available deep JWST NIRCam data covering
the PEARLS, CEERS, GLASS, JADES GOOD-S, NGDEEP, and SMACS0723 surveys,
totalling 187 arcmin2. We investigate the impact of SED fitting methods,
assumed star formation histories (SFH), dust laws, and priors on galaxy masses
and the resultant GSMF. Whilst our fiducial GSMF agrees with the literature at
$z<13.5$, we find that the assumed SFH model has a large impact on the GSMF
and stellar mass density (SMD), finding a 0.75 dex increase in the SMD at
$z=10.5$ between a flexible non-parametric and standard parametric SFH.
Overall, we find a flatter SMD evolution at $z\geq 9$ than some studies
predict, suggesting a rapid buildup of stellar mass in the early Universe. We
find no incompatibility between our results and those of standard cosmological
models, as suggested previously, although the most massive galaxies may
require a high star formation efficiency. We find that the ‘Little Red Dot’
galaxies dominate the $z=7$ GSMF at high-masses, necessitating a better
understanding of the relative contributions of AGN and stellar emission. We
show that assuming a theoretically motivated top-heavy IMF reduces stellar
mass by 0.5 dex without affecting fit quality, but our results remain
consistent with existing cosmological models with a standard IMF.

galaxies: photometry, high-redshift, evolution, statistics

††facilities: JWST (STScI), HST (STScI)

tablewidth xtablewidth

##  1 Introduction

The James Webb Space Telescope (JWST) has pushed backed the redshift frontier
when searching for the earliest galaxies. The highly sensitive Near Infrared
Camera (NIRCam) on JWST has led to an influx of high redshift galaxy
candidates through photometry, in surveys such as CEERS, GLASS, PEARLS, NGDEEP
& JADES (Adams et al., 2023; Castellano et al., 2022; Windhorst et al., 2023;
Austin et al., 2023; Hainline et al., 2023; Bagley et al., 2023a). The
wavelength coverage from 0.6 to 5$\mu$m, enables identification of Lyman-break
galaxies (LBGs) at redshifts $z\geq 6.5$ by their photometry. In the first 18
months of science operations, tens of candidates above $z\geq 10$ have been
identified (Adams et al., 2023; Donnan et al., 2022; Atek et al., 2023; Naidu
et al., 2022, 2022; Harikane et al., 2023; Castellano et al., 2022;
Finkelstein et al., 2023a; Austin et al., 2023; Hainline et al., 2023; Furtak
et al., 2023; Labbé et al., 2023; Pérez-González et al., 2023; Finkelstein et
al., 2023b; Leung et al., 2023; McLeod et al., 2023; Willott et al., 2023a),
including a spectroscopically confirmed galaxy at $z=13.27$ (Curtis-Lake et
al., 2023) and candidates at $z\geq 16$ (e.g. Atek et al., 2023; Yan et al.,
2023; Furtak et al., 2023; Hainline et al., 2023).

An immediate result was a potential overabundance of high-redshift galaxies
compared to theoretical predictions and extrapolations of Hubble Space
Telescope (HST)/Spitzer results (Naidu et al., 2022; Haslbauer et al., 2022;
Mauerhofer & Dayal, 2023; Mason et al., 2023). The mere existence of some of
the galaxies at the inferred redshifts and stellar masses, ($\geq 10^{10.5}$
M⊙ at $z\geq 7.5$) seem to be in tension with standard $\Lambda$CDM cosmology
given the small areas currently probed with JWST (Labbé et al., 2023; Lovell
et al., 2023; Boylan-Kolchin, 2023). However, spectroscopic observations of
some of the highest-mass candidates in Labbé et al. (2023) have reduced their
redshifts and stellar masses, or hinted at the presence of an Active Galactic
Nuclei (AGN), and hence reduced possible $\Lambda$CDM tension (Kocevski et
al., 2023; Fujimoto et al., 2022).

Initial overlap of galaxy candidates between independent studies was poor, but
has since improved due to agreement on photometric calibration and reduction
techniques (Adams et al., 2023; Rieke et al., 2022). Spectroscopic
confirmation exists only for a fraction of potential candidates, but most
spectroscopic programs have had a high success rate, along with a few notable
low-$z$ interlopers (Robertson et al., 2023b; Curtis-Lake et al., 2023;
Arrabal Haro et al., 2023; Tang et al., 2023; Laseter et al., 2023; Fujimoto
et al., 2023; Bunker et al., 2023a; Wang et al., 2023a).

The combination of NIRCam’s high sensitivity and infrared (IR) wavelength
coverage allows characterisation of the rest-frame optical emission of
$0.5\leq z\leq 10$ galaxies, which is crucial for accurate stellar mass
estimates. A more complete census of the high-redshift Universe is also
possible, as galaxies without a strong Lyman-break (so called ‘HST-dark’
galaxies) were often missed in the Ultraviolet (UV) selected samples of pre-
JWST studies (Pérez-González et al., 2023). Observations with Spitzer IRAC
were available only for the brightest sources due to low sensitivity and
angular resolution. Intrinsically UV-faint galaxies are often found to be
dusty or evolved systems, and accurate characterisation of this population is
essential when measuring the total buildup of stellar mass in the Universe,
which is typically done by measuring the Galaxy Stellar Mass Function (GSMF).

Stellar masses are typically estimated from fitting spectral energy
distributions to broadband photometry, and inferring a star formation history
along with other physical parameters such as metallicity, dust obscuration and
ionization state of the gas (e.g., Brinchmann & Ellis, 2000; Bundy et al.,
2006; Mortlock et al., 2011; Duncan et al., 2014; Mortlock et al., 2015;
Grazian et al., 2015; Song et al., 2016; Bhatawdekar et al., 2019). This
approach has generally been found to be reliable in the local Universe. At
high-redshift there are a number of complicating factors that must be
considered, which have different effects but overall act to increase the
uncertainty in stellar masses beyond the statistical uncertainty from the
fitting (Lower et al., 2020; Wang et al., 2023b). There is a growing consensus
that the star formation histories of many high-$z$ galaxies are stochastic and
characterised by repeated cycles of a short burst of rapid star formation
followed by a temporary period of quiescence (Faucher-Giguère, 2018; Looser et
al., 2023; Asada et al., 2023; Dressler et al., 2023; Dome et al., 2024).
Young, bright stars created in the most recent bursts of star formation can
dominate the spectral energy distributions of galaxies and obscure older
populations, known as ’outshining’ and lead to the stellar mass being
underestimated by up to 1 dex (Giménez-Arteaga et al., 2023; Endsley et al.,
2021; Pérez-González et al., 2023; Papovich et al., 2022; Jain et al., 2024;
Giménez-Arteaga et al., 2024). For extremely stochastic SFHs, information loss
of the first periods of star formation may occur, leading to large
uncertainties in stellar mass (Narayanan et al., 2023; Shen et al., 2023; van
Mierlo et al., 2023). Flexible ‘non-parametric’ SFHs, such as those presented
in Leja et al. (2019); Tacchella et al. (2022); Robertson et al. (2023b), are
able to reproduce these bursty SFHs more accurately at high-$z$, typically
finding systematically larger stellar masses, (Giménez-Arteaga et al., 2023;
Jain et al., 2024), and may produce more reliable stellar mass estimates when
compared to traditional parametric SFHs, which typically vary smoothly
(Carnall et al., 2019, e.g. delayed exponential).

Another assumption is that of a possible universal initial mass function
(IMF), which predicts the number of stars as a function of stellar mass
(Salpeter, 1955). The presence of low-mass stars, which dominate the stellar
mass, is inferred entirely from the shape of the assumed IMF in most galaxies.
A universal IMF (e.g. Salpeter, 1955; Chabrier et al., 2000; Kroupa, 2001),
which is constant across time and space, has long been assumed in the majority
of galaxy studies. Models of the the physics of high-$z$ star formation
suggest that it may have deviated from the universal IMF above redshift $z=8$
(Hopkins et al., 2005; Jermyn et al., 2018; Steinhardt et al., 2021, 2023).
Low metallicity, high star formation rates, an increasing CMB temperature and
high cosmic ray density could all contribute to heating of star formation
regions and lead to a top-heavy IMF at high redshift (Gunawardhana et al.,
2011; Clauwens et al., 2016; Papadopoulos et al., 2011). Observations of local
elliptical galaxies have also found them to be inconsistent with the universal
IMF, and instead find evidence for an IMF gradient, with evidence for a
different IMF between younger and older stellar populations (e.g. Weidner et
al., 2013, & references therein). Studies such as Steinhardt et al. (2023);
Sneppen et al. (2022); Woodrum et al. (2023) have introduced temperature-
dependent IMF models for use in SED-fitting, which can decrease stellar mass
estimates by up to 1 dex at constant redshift. Constraining the IMF is
extremely difficult, but recent studies are beginning to find possible
evidence for a top-heavy IMF in the early Universe (e.g. Katz et al., 2022;
Cameron et al., 2023; Mowla et al., 2024). Ultimate conclusions on this are
however far from certain.

In the EPOCHS paper series we have presented an independent and consistent
reduction of deep JWST observations from available GTO, GO and ERS data
including the CEERS, GLASS, SMACS-0723, JADES and PEARLS fields. We found 1165
robust galaxy candidates above redshift $z=6.5$, with a total area of 187
arcmin2 (Conselice et al., in prep; Adams et al., 2023).

In this paper we present a detailed examination of the inferred physical
parameters of our high-redshift sample, with a particular focus on the stellar
mass of our galaxy candidates. We derive a galaxy stellar mass function at
$6.5\leq z\leq$13.5, and estimate the stellar mass density in order to trace
the buildup of stellar mass in the early Universe. Importantly, we explore the
impact of some of the key assumptions used in deriving stellar masses at high
redshifts, such as star formation histories and the IMF.

In § 2 we present a brief overview of the data products used in this work and
detail our data reduction procedure. § 3 describes our process for catalogue
creation and robust sample selection of high-$z$ galaxy candidates. We detail
our SED-fitting procedure and the impact of different assumptions on the
inferred properties of high-$z$ galaxies in § 4. In § 5 we use our stellar
mass estimates to build on the UV luminosity function presented in Adams et
al. (2023) and construct a galaxy stellar mass function at redshifts from 6.5
$\leq z\leq$ 13.5. We discuss our findings and make comparisons to the
literature in § 6. Finally § 7 summarises the findings of this work and looks
at possibilities for future studies.

We assume a standard $\Lambda\mathrm{CDM}$ cosmology with $H_{0}=70$ km s-1
Mpc-1, $\Omega_{\rm M}=0.3$ and $\Omega_{\Lambda}=0.7$. All magnitudes listed
follow the AB magnitude system (Oke, 1974; Oke & Gunn, 1983). All stellar
masses measured in this work use a Kroupa (2002) unless otherwise indicated.

##  2 Products and Data Reduction

This section briefly details the JWST programs and data products used in this
analysis. For further details on the fields used please see Austin et al.
(2023), Adams et al. (2023) and Conselice et al. (in prep). Table 1 shows the
available unmasked areas, JWST NIRCam filters and 5$\sigma$ depths of each
dataset used.

Table 1: Table showing the unmasked areas and depths of the observations used in this work. Depths are given at 5$\sigma$ in AB magnitudes, measured in 0$\farcs$16 radius apertures. Depths are calculated by placing non-overlapping apertures in empty regions of the image, as measured by the SExtractor segmentation maps and our image masks. The nearest 200 apertures are used to calculate the Normalised Mean Absolute Deviation (NMAD) to derive local depths for each individual source. Where depths are tiered across mosaics (e.g. HST Advanced Camera for Surveys Wide Field Channel (ACS/WFC) observations in the Hubble Ultra Deep Field (HUDF) Parallel 2) we have listed the depths and areas separately. The four spokes of the NEP-TDF and ten CEERS pointings have uniform depths (within 0.1 mags) with the exception of CEERS P9 which we list separately. Areas are given in arcmin2 and measured from the mask to account for the masked areas of the image and unused cluster modules. Fields with a ‘*’ indicate that we have excluded the NIRCam module containing a lensing cluster from our analysis. | Area | HST/ACS_WFC | JWST/NIRCam
---|---|---|---
Field | (arcmin2) | F606W | F814W | F090W | F115W | F150W | F200W | F277W | F335M | F356W | F410M | F444W
NEP | 57.32 | 28.74 | - | 28.50 | 28.50 | 28.50 | 28.65 | 29.15 | - | 29.30 | 28.55 | 28.95
El Gordo* | 3.90 | - | - | 28.23 | 28.25 | 28.18 | 28.43 | 28.96 | - | 29.02 | 28.45 | 28.83
MACS-0416* | 12.3 | - | - | 28.67 | 28.62 | 28.49 | 28.64 | 29.16 | - | 29.33 | 28.74 | 29.07
CLIO* | 4.00 | - | - | 28.12 | - | 28.07 | 28.21 | 28.675 | - | 28.91 | - | 28.71
CEERS | 66.40 | 28.6 | 28.30 | - | 28.70 | 28.60 | 28.89 | 29.20 | - | 29.30 | 28.50 | 28.85
CEERSP9 | 6.08 | 28.31 | 28.32 | - | 29.02 | 28.55 | 28.78 | 29.20 | - | 29.22 | 28.50 | 29.12
SMACS-0723* | 4.31 | - | - | 28.75 | - | 28.81 | 28.95 | 29.45 | - | 29.55 | - | 29.28
GLASS | 9.76 | - | - | 29.14 | 29.11 | 28.86 | 29.03 | 29.55 | - | 29.61 | - | 29.84
NGDEEP HST-S | 1.28 | 29.20 | 28.80 | - | 29.78 | 29.52 | 29.48 | 30.28 | - | 30.22 | - | 30.22
NGDEEP HST-D | 4.03 | 30.30 | 30.95 | - | 29.78 | 29.52 | 29.48 | 30.28 | - | 30.22 | - | 30.22
JADES Deep GS | 22.98 | 29.07 | - | 29.58 | 29.78 | 29.68 | 29.72 | 30.21 | 29.58 | 30.17 | 29.64 | 29.99

###  2.1 PEARLS

We incorporate NIRCam observations from the proprietary GTO survey Prime
Extragalactic Areas for Reionization Science (PEARLS, PI: R. Windhorst &
H.Hammel, PID: 1176 & 2738, Windhorst et al., 2023). We use observations of
three fields targeting gravitationally lensed clusters, and one blank field
consisting of a mosaic of 8 JWST NIRCam pointings. The gravitationally lensed
clusters consist of MACS J0416.1-2403 (hereafter referred to as MACS-0416), El
Gordo (Menanteau et al., 2012, $z\sim 0.87$, ACT-CL J0102-4915 in the Atacama
Cosmology survey, [) and Clio ($z\sim 0.42$, Designation GAMA100050 in the
GAMA Galaxy Group Catalog v6+, Robotham et al., 2011). El Gordo and Clio have
been visited once with JWST/NIRCam, with the cluster centered in one NIRCam
module and the other observing a neighbouring blank field $\sim 3\arcmin$ away
(Griffiths et al., 2018). MACS-0416 has been observed 3 times, resulting in 3
separate parallel observations at different position angles. The PEARLS blank
field is the North Ecliptic Pole Time Domain Field (Jansen & Windhorst, 2018,
(NEP-TDF,). The NEP-TDF is positioned so it can be observed throughout the
year, making it ideal for time-domain science and constructing a large deep
field. Observations of the NEP-TDF consist of four pairs of overlapping NIRCam
pointings (8 pointings total), with each of these four pairs orientated at 90
degree intervals like spokes on a windmill. NIRCam observations of the NEP-
TDF, El Gordo and MACS0416 use the standard 8 photometric bands; F090W, F115W,
F150W, F200W, F277W, F356W, F410M and F444W. Clio uses 6 of the 8 previous
bands, but lacks F115W and F410M. Within the NEP-TDF field we incorporate HST
Advanced Camera for Surveys Wide Field Channel (ACS/WFC) imaging in the F606W
filter, collected as part of the GO-15278 (PI: R. Jansen) and GO-16252/16793
(PIs: R. Jansen & N. Grogin) between October 1 2017 and October 31 2022.
Mosaics of these data, astrometrically aligned to Gaia/DR3 and resampled on
0$\farcs$03 pixels, were made available pre-publication by R. Jansen & R.
O’Brien (private comm. O’Brien et al., 2024). For full details of the PEARLS
program please see Windhorst et al. (2023).

###  2.2 ERS and GO Data

We incorporate Early Release Science (ERS) and public General Observer (GO)
data from SMACS-0723 (PID: 2736, PI: K. Pontoppidan, Pontoppidan et al.,
2022), the Cosmic Evolution Early Release Science Survey (CEERS, PID: 1345,
PI: S. Finkelstein, see also Bagley et al., 2023b), the Grism Lens Amplified
Survey from Space survey (GLASS, PID: 1324, PI: T. Treu, Treu et al., 2022)
and the Next Generation Deep Extragalactic Exploratory Public Survey (NGDEEP,
PID: 2079, PIs: S. Finkelstein, Papovich and Pirzkal, Bagley et al., 2023a).
We incorporate HST ACS/WFC observations of the Extended Groth Strip (EGS Davis
et al., 2007) into our CEERS dataset in the F606W and F814W filters. This was
obtained as part of the Cosmic Assembly Near-infrared Deep Extragalactic
Legacy Survey (CANDELS, Grogin et al., 2011; Koekemoer et al., 2011), with
updated astrometric alignment to Gaia EDR3 (Brown et al., 2021) by the CEERS
team111https://ceers.github.io/hdr1.html and released as Hubble Data Release
1. The addition of these observations compensates for the lack of F090W
observations in the CEERS survey.

We also incorporate NIRCam imaging of the Great Observatories Origins Deep
Survey South (GOODS-South) field collected as part of JWST Advanced Deep
Extragalactic Survey (JADES, PID:1180, PI: D. Eisenstein, Eisenstein et al.,
2023a) and released publicly as JADES DR1 (Rieke et al., 2023). In the JADES
and NGDEEP fields, which lie on the GOODS-South footprint, we add in existing
HST data from F606W and F814W from the most recent mosaic (v2.5) from the
Hubble Legacy Fields team (Illingworth et al., 2016; Whitaker et al., 2019).

###  2.3 JWST NIRCam Data Reduction

We have uniformly reprocessed all lower-level JWST data products following our
modified version of the official JWST pipeline. This is a similar process to
that used in Ferreira et al. (2022), Adams et al. (2023), Austin et al.
(2023), and in particular Adams et al. (2023) but with updates based on new
flat-fielding and techniques for dealing with NIRCam imaging artefacts.

We use version 1.8.2 of the official STScI JWST
Pipeline222https://github.com/spacetelescope/jwst (Bushouse et al., 2022) and
Calibration Reference Data System (CRDS) v1084, which contains the most up-to-
date NIRCam calibrations at the time of writing and includes updated flat-
field templates for the LW detectors, resulting in improved average depths
across a single pointing of $\sim$ 0.2 dex in F444W. Next we subtract
templates of wisps, artefacts present in the F150W and F200W imaging, between
stage 1 and stage 2 of the pipeline. After stage 2 of the pipeline, we apply
the $1/f$ noise correction derived by Chris
Willott333https://github.com/chriswillott/jwst, which removes linear features
caused by read-noise from the images. We do not use the sky subtraction step
included in stage 3 of the pipeline, instead performing background subtraction
on individual NIRCam frames between stage 2 and stage 3 (‘cal.fits’ files),
consisting of an initial uniform background subtraction followed by a
2-dimensional background subtraction using photutils (Bradley et al., 2022).
This allows for quicker assessment of the background subtraction performance
and immediate fine-tuning of configuration parameters. After stage 3 of the
pipeline, we align the final F444W image onto a Gaia-derived World Coordinate
System (WCS) (Gaia Collaboration et al., 2018; Vallenari et al., 2022) using
tweakreg, part of the DrizzlePac python
package444https://github.com/spacetelescope/drizzlepac, and then match all
remaining filters to this derived WCS, ensuring the individual images are
aligned to one another. In some cases (NEP, CEERS), we match to a WCS-frame
derived from other space or ground based imaging with a larger FOV, given the
low number of Gaia stars in some individual NIRCam pointings. We then pixel-
match the images to the F444W image with the use of astropy, reproject
(Astropy Collaboration et al., 2013, 2022a; Hoffmann et al.,
2021).555https://reproject.readthedocs.io/en/stable/ The final resolution of
the drizzled images is 0$\farcs$03/pixel. Comparison of our reduction to the
official PEARLS reduction pipeline (Windhorst et al., 2023) is given in Adams
et al. (2023), finding excellent agreement in both observed fluxes (within
0.03 (0.01) magnitudes in the blue (red) NIRCam photometric filters) and
astrometry (within 2 pixels (0$\farcs$07)).

##  3 Catalogue Creation and Sample Selection

Full details of our catalogue creation and sample selection pipeline, called
GALFIND, is available in Conselice et. al. (in prep). We briefly summarise the
procedure used here.

###  3.1 Catalogue Creation

We use the code SExtractor (Bertin & Arnouts, 1996) for source identification
and photometric measurements. We use an inverse-variance weighted stack of the
NIRCam F277W, F356W and F444W images for source detection in order to reliably
identify faint sources and then carry out forced aperture photometry in all
photometric bands. This photometry is calculated within 0.32 arcsecond
diameter circular apertures, correcting for the aperture size with an aperture
correction derived from simulated WebbPSF point spread functions (PSFs) for
each band used (Perrin et al., 2012, 2014). This diameter was chosen to
enclose the central/brightest $70-80$ per cent of the flux of a point source
without a large amount of contamination from neighbouring sources. This
reduces the reliance on a potentially uncertain PSF model whilst still using
the brightest pixels when calculating fluxes.

As SExtractor requires all images to be on the same pixel grid, for aligned
HST imaging on a different pixel grid we use photutils to perform forced
aperture photometry in the same diameter apertures (Bradley et al., 2022).

We next produce masks for our images by eye to cover diffraction spikes, any
remaining snowballs, the cross pattern between SW detectors, image edges
(including a $\sim 50-100$ pixel border around detector edges) and any large
foreground galaxies. The total amount of unmasked area used in this study is
listed alongside the average 5$\sigma$ depths of each field in Table 1.

Following the generation of source catalogues and segmentation maps for each
image, we calculate local depths for each source in each filter. This accounts
for variation in background and noise across the image. Apertures of
0$\farcs$32 diameter are placed in empty regions of the image, as calculated
from the segmentation map to be $\geq 1\arcsec$ from the nearest source. For
each source the nearest 200 apertures are used to calculate the Normalised
Mean Absolute Deviation (NMAD) of the fluxes measured in the apertures, which
corresponds to the 1$\sigma$ flux uncertainty. We convert this to a 5$\sigma$
depth, displaying the average depth in AB magnitudes for each field in Table
1. Where fields consist of mosaics of multiple pointings we display the
average depth across the entire field, but we note that 9/10 CEERS pointings
and 4/4 NEP-TDF have depths consistent within 0.1 magnitudes. The exception in
CEERS is pointing 9 (P9), which has an additional exposure in F115W and F444W
resulting in increased depths of $\sim 0.2$ mags (Adams et al., 2023). Due to
correlated noise, the flux uncertainties calculated by SExtractor are
underestimated (Schlawin et al., 2020). We replace these uncertainty estimates
with the local-depth derived flux errors (Adams et al., 2023).

###  3.2 Sample Selection

To select a sample of high-redshift galaxies we introduce selection criteria
based primarily on photometric SED-fitting with EAZY-py (Brammer et al.,
2008). We aim to select a robust sample of galaxies above z$\geq 6.5$ , where
the Lyman break is within the NIRCam F090W filter.

We use the default EAZY templates (tweak_fsps_QSF_12_v3), along with Set 1 and
Set 4 of the SED templates generated by Larson et al. (2022). These additional
templates were developed to have bluer rest-frame UV colors than the default
templates, as well as stronger emission lines, both of which have been
observed in high-redshift galaxies (Finkelstein et al., 2022; Nanayakkara et
al., 2022; Cullen et al., 2023; Withers et al., 2023). These templates have
young stellar populations, low metallicities and active star formation. Larson
et al. (2022) has shown that they improve the accuracy of photo-$z$ estimates
at high redshift.

We run EAZY-py initially with a uniform redshift prior of $0\leq z\leq 25$,
but then repeat the fitting with a reduced upper redshift limit of $z\leq 6$.
This allows us to compare the goodness of fit of both a high and low-redshift
solution for all galaxies in our sample. We use a minimum flux uncertainty of
10% to account for uncertainties in flux calibrations and aperture corrections
(Rieke et al., 2022), as well as potential differences between the synthetic
templates and our galaxies.

For reproducibility our selection criteria are designed to be based as much as
possible on specific cuts in computed quantities, rather than individual
inspection of candidates which can introduce hard to measure bias and
incompleteness. To ensure robustness in our sample, our final selection
criteria includes a visual review of the cutouts and SED-fitting solutions for
all sources by authors TH, DA, NA, & QL but we reject less than 5% of our
original sample by eye at this stage, which is much lower than comparable
studies (reaching $\geq$50% in some cases, e.g. Hainline et al., 2023).

Our selection criteria for robust high redshift galaxies are as follows.

  1. 1.

We require that the bandpass of the lowest wavelength photometric band must be
entirely below the Lyman break, given the primary photo-$z$ solution. This
sets a lower limit of $z\approx 6.5$ in most of our fields.

  2. 2.

We require a $\leq 3\sigma$ detection in band(s) blueward of the Lyman break.

  3. 3.

We require a $\geq 5\sigma$ detection in the 2 bands directly redward of the
Lyman break, and $\geq 2\sigma$ detection in all other redward bands,
excluding observations in NIRCam medium-bands (e.g. F335M, F410M). If the
galaxy appears only in the long wavelength NIRCam photometry (i.e. a F200W or
higher dropout), we increase the requirement on the 1st band to 7$\sigma$
detection.

  4. 4.

The integral of the photo-$z$ PDF is required to satisfy $\int^{1.10\times
z_{\textrm{p}hot}}_{0.90\times z_{\textrm{p}hot}}\ P(z)\ dz\ \geq\ 0.6$ to
ensure the majority of the redshift PDF is located within the primary peak,
and that the peak is sufficiently narrow to provide a strong constraint on the
redshift. $z_{\textrm{p}hot}$ refers to the redshift with maximum likelihood
from the EAZY-py redshift posterior.

  5. 5.

We require the best-fitting EAZY-py SED to satisfy $\chi^{2}_{red}<3(6)$ to be
classed as a robust (good) fit.

  6. 6.

We require a difference of $\Delta\chi^{2}\geq$ 4 between the high-$z$ and
low-$z$ EAZY-py runs. This ensures that the high-$z$ solution is much more
statistically probable.

  7. 7.

If the half-light radius (FLUX_RADIUS parameter in SExtractor) is smaller than
the FWHM of the PSF in the F444W band, then we require that
$\Delta\chi^{2}\geq$ 4 between the best-fitting high-$z$ galaxy solution and
the best-fitting brown dwarf template. This requirement is designed to remove
Milky Way brown dwarf contaminants and is discussed further in § 5.2.2.

  8. 8.

We require the 50% encircled flux radius from SExtractor to be $\geq$1.5
pixels in the long-wavelength wideband NIRCam photometry (F277W, F356W,
F444W). This avoids the selection of oversampled hot pixels in the LW
detectors as F200W dropouts.

Given our requirement to observe the Lyman break, the lowest redshift at which
we select ‘robust’ galaxies with NIRCam photometry only is $\sim$6.5, where
the break falls within the NIRCam F090W filter. In the fields where HST
ACS/WFC imaging is available, we can robustly identify the Lyman break at
lower redshifts.

Our selection criteria are similar to other high-$z$ galaxy studies, such as
Hainline et al. (2023); Finkelstein et al. (2023a, b); Naidu et al. (2022);
Castellano et al. (2022), who also fit galaxies using EAZY-py, and select
their samples from the rest-UV SNR of their candidates and the resultant
redshift PDF and $\Delta\chi^{2}$ from the SED fitting. SNR requirements vary
between studies, and are also somewhat dependent on the size of the extraction
aperture used. Our 0$\farcs$16 radius apertures are larger than the 0$\farcs$1
radius apertures used in Finkelstein et al. (2023b), resulting in a number of
their candidates being removed from our sample because they do not meet our
5$\sigma$ threshold, despite agreeing on the photo-$z$ solution. Other
studies, such as as Harikane et al. (2023); Yan et al. (2022), use a
combination of color-color cuts and SED-fitting. A comprehensive comparison of
the EPOCHS v1 galaxy sample to other studies is given in Appendix A of Adams
et al. (2023).

Spectroscopy of high redshift galaxies with NIRSpec (e.g. Wang et al., 2023a;
Fujimoto et al., 2023; Bunker et al., 2023a; Curtis-Lake et al., 2023) has
shown that many high redshift galaxies have strong emission lines, with
H$\beta$+[OIII] reaching observed equivalent widths of up to 2000-3000Å
(Withers et al., 2023). For our photometric observations this can result in an
excess in the band covering these emission lines of up to $\sim$ 1 dex. The
emission line modelling in the SED fitting codes used span only a limited
parameter space in equivalent width and line ratios, so this can result in
high $\chi^{2}$ values even if the model is a good fit to the rest of the
photometry. To avoid removing these galaxies from our sample we introduce a
secondary group, referred to as ‘good’ galaxies, which have $3<\chi_{r}^{2}<6$
but meet the rest of our criteria. This applies to a very small fraction of
our total sample, with only 23 galaxies meeting our other selection criteria
but falling within this $\chi^{2}$ range.

##  4 Galaxy Properties from SED Fitting

The basic properties, redshift distribution, number counts and UVJ colors of
the EPOCHS v1 sample are described in detail in Conselice et al. (in prep).
The UV properties ($M_{\mathrm{UV}}$, $\beta$ slopes) of this sample are
explored in Austin et al. (in prep). The UV luminosity function is presented
in Adams et al. (2023), Here we briefly summarise the basic statistics of the
EPOCHS v1 sample, following the selection criteria described in § 3.2.

Across the fields described in § 2, 1214 galaxies pass our initial selection
criteria with $z\geq 6.5$ over a total unmasked area of 187.27 arcmin2,
comprising one of the largest samples of high-$z$ JWST-selected galaxies. Our
visual inspection removes 49 completely, leaving 1165 galaxies, which we refer
to as the EPOCHS v1 sample. Of these 1165, 1054 are classed as “certain” and
111 “uncertain” by our visual inspection. In this work we choose to include
the visually “uncertain” candidates, in order to avoid potentially
underestimating the stellar mass function. Filtering the sample to the
$6.5<z<13.5$ redshift range used in this work, our fiducial sample consists of
1120 galaxies.

The next sections detail the Bagpipes and Prospector SED fitting performed for
all galaxies in the EPOCHS v1 sample. We perform this SED fitting only for our
EAZY-py-selected sample. The purpose of this SED-fitting is to analyse the
properties of the stellar and nebular components of these high-$z$ galaxies.
We note that whilst the majority of galaxies in the EPOCHS v1 sample do not
appear significantly extended beyond our SExtractor extraction aperture, all
stellar masses quoted in the following results have been corrected by the
ratio between the SExtractor FLUX_AUTO_F444W and the aperture corrected
FLUX_APER_F444W fluxes, if it exceeded unity, to account for any residual flux
outside the aperture. We use the longest wavelength band to correct stellar
masses as this is most representative of the rest-frame optical emission.

###  4.1 Bagpipes

Bagpipes (Bayesian Analysis of Galaxies for Physical Inference and Parameter
EStimation, Carnall et al., 2018, 2019)) is a Python package which uses
Bayesian methods to fit galaxy SEDs to photometry. Our fiducial Bagpipes run
uses the Bruzual & Charlot (2003) 2016 stellar population synthesis models,
and a Kroupa (2001) IMF.

Bagpipes can construct and fit SEDs with a variety of SFH models, dust models
and priors. We perform multiple fits for each galaxy in order to test the
consistencies in derived galaxy parameters. For simplicity our approach is to
define an initial ‘fiducial’ model, and then swap out individual model
components or priors for other choices. Our model parameters, priors and
hyper-parameters are detailed in Table 2. This is similar to the approach of
Furtak et al. (2021), who define a reference model which is then compared to
alternative models. In the rest of this section we further explain our choices
of models and priors.

For the star formation history we test both parametric and non-parametric
models, which have been shown to impact the stellar mass estimates (Leja et
al., 2019; Tacchella et al., 2022). For our fiducial model, we use a
parametric lognormal SFH which allows us to recreate the rising SFHs we expect
in the early Universe. We compare this SFH with another commonly used
parametrization, the delayed exponential. Details of the implementation of
these parametric SFHs are available in Carnall et al. (2019). We also test a
non-parametric “continuity” SFH model similar to the model added to Prospector
in Leja et al. (2019). This SFH model fits the star formation rates in fixed
time bins, with $\Delta\log$SFR between bins linked by a Student’s
t-distribution. We recreate the methodology of Tacchella et al. (2022), by
fitting both a “continuity” model, where the Student’s t-distribution has
hyper-parameters $\sigma$ = 0.3 and $\nu$= 2, which weights against rapid
changes in star formation rate, as well as a “continuity bursty” model, with
$\sigma$ = 1 and $\nu$= 2, which allows more stochastic star formation, with
rapid bursts and quenching more similar to the SFH inferred at high-$z$. As is
done in Tacchella et al. (2022) we fit 6 SFH bins for both models, with the
first bin fixed to a lookback time of 0-10 Myr, the last bin ending at z = 20,
and the other 4 bins are equally log-spaced in lookback time. Leja et al.
(2019) showed that this model is relatively insensitive to the number of bins
used, as long as the number exceeds 4. Like Tacchella et al. (2022), we assume
there is no star formation at $z>20$ (Tacchella et al., 2022; Bowman et al.,
2018; Jaacks et al., 2019). It is important to note that this model allows but
does not require ‘bursty’ SFHs, and smooth or quenched SFHs are also possible
if favoured during the fitting.

We include emission lines and nebular continuum based on CLOUDY v17.03
(Ferland et al., 2017). We regenerate CLOUDY models in order to probe a wider
range of the ionisation parameter U between -3 $\leq\log_{10}U\leq$ -1 using a
CLOUDY configuration file distributed with Bagpipes. For the ionization
parameter we use a uniform prior (in log space).

We use the Calzetti et al. (2000) prescription for dust. Bowler et al. (2023)
finds that UV-selected high-$z$ galaxies in the ALMA REBELS survey follow the
local Calzetti-like IRX-$\beta$ relation, so we do not fit a more complex dust
law in our fiducial model. The allowed stellar metallicity ranges from -2.3
$\leq\log_{10}Z_{\star}/Z_{\odot}\leq$ 0.70 with a logarithmic prior, as these
galaxies are expected to have low metallicity, but theoretically could enrich
their local environments quickly (Langeroodi et al., 2023; Curti et al.,
2023).

In order to constrain the redshift parameter space we fix the redshift prior
to the PDF from our EAZY-py SED-fitting, which we approximate as a Gaussian.
The redshift prior draws are capped at $\pm 3\sigma$. We use the default
sampling method, using the Python package PyMultiNest (Feroz et al., 2009;
Buchner, 2016).

Table 2: Summary of parameters, hyper-parameters and priors for our Bagpipes
SED-fitting. The text color corresponds to a specific Bagpipes run, where
black corresponds to our default ‘fiducial’ setting for each parameter.
Parameters and priors for other iterations can be assumed to be the same as
given for the ‘fiducial’ bagpipes run unless otherwise specified. The top
section of the table lists parameters that are common to all of our Bagpipes
models, whereas the lower section gives the model-specific parameters for each
of our chosen configurations.

Common Parameters |  |
---|---|---
Parameter | Prior/Value (Min, Max) | Description
$z_{\textrm{phot}}$ |  EAZY-py Posterior PDF ($\pm 3\sigma$) | Redshift
SPS Model |  Bruzual & Charlot (2003);  BPASS v2.2.1 | Stellar population synthesis model
IMF |  Kroupa (2001);  default BPASS IMF | Stellar Initial Mass Function
Dust Law parametrization |  Calzetti et al. (2000); Salim et al. (2018) | Dust law
| Charlot & Fall (2000) |
A${}_{\textrm{V}}$ | Log-uniform: ($10^{-3}$, 10); uniform (0, 6) | V-band attenuation (all stars)
SFH | lognormal; “continuity bursty”; delayed-$\tau$ | Star formation history
$\log_{10}(M_{\star}/M_{\odot})$ | uniform: (5, 12) | Surviving stellar mass
Z${}_{\star}/\textrm{Z}_{\odot}$ | log-uniform: (0.005, 5); uniform (0, 3) | stellar metallicity
Z${}_{\textrm{gas}}/\textrm{Z}_{\odot}$ | Fixed to Z⋆ | gas-phase metallicity
$\log_{10}\textrm{U}$ | uniform: (-3, -1) | Ionization Parameter

Model Specific Parameters |  |  |
---|---|---|---
Model | Parameter | Prior/Value (Min, Max) | Description
Fiducial | t${}_{\textrm{max}}$ | uniform: (10 Myr, 15 Gyr) | Age of Universe at peak SFR
| FWHM | uniform: (10 Myr, 15 Gyr) | FWHM of SFH
delayed-$\tau$ SFH | $\tau$ | uniform: (10 Myr, 15 Gyr) | e-folding timescale
| Age | log-uniform: (10 Myr, t${}_{\rm univ}(z_{\rm phot}))$ | Time since SF began
“‘continuity bursty”’ | N${}_{\textrm{bins}}$ | 6 bins (5 fitted parameters) | First bin 0 – 10 Myr, SF begins at z = 20,
non-parametric SFH |  |  | others distributed equally in log10 lookback time
| d${}_{\log_{10}\textrm{SFR}}$ | Student’s-t: $\nu=2,\sigma=1.0$ | Ratio of $\log_{10}$SFR in adjacent bins, coupled by $\sigma$
Charlot & Fall (2000) | n | clipped normal: $\mu=0.7,\sigma=0.3$ | Power-law slope of attenuation curve (A $\propto\lambda^{-n}$)
Dust Law |  | (0.3, 2.5) | For Calzetti et al. (2000) n $\approx 0.7$
| $\eta$ | clipped normal: $\mu=2,\sigma=0.3$ | AV,<10Myr/AV ratio between young and old stars
|  | (1, 3) |
Salim et al. (2018) Dust Law | $\delta$ | clipped normal: $\mu=0,\sigma=0.1$ | Deviation from Calzetti et al. (2000) slope
|  | (-0.3, 0.3) |
| $\beta$ | uniform (0, 5) | Strength of 2175Å bump
BPASS SPS Model |  |  | No additional components
Uniform A${}_{\textrm{V}}$ Prior |  |  | No additional components
Uniform Z⋆ Prior |  |  | No additional components

###  4.2 Prospector

In addition to our Bagpipes SED-fitting we also fit our sample using the
Prospector package (Johnson et al., 2021) in order to compare the results of
these two commonly used SED fitting tools. Prospector allows greater
flexibility and control of model parameters than Bagpipes at the expense of
computational time. Prospector uses Bayesian inference to determine galaxy
stellar population properties and star formation histories. Prospector is
built on Flexible Stellar Population Synthesis (FSPS, Conroy & Gunn, 2010)),
using python-fsps (Johnson et al., 2023) and the Modules for Experiments in
Stellar Astrophysics Isochrones and Stellar Tracks (MIST) stellar isochrones
(Choi et al., 2016; Dotter, 2016). Prospector allows for very flexible star
formation histories, including the use of non-parametric SFHs, which have
become increasingly popular in the JWST era. Non-parametric SFHs are modelled
by splitting star-formation history into a number of temporal bins, with the
SFH in each bin linked by a statistical distribution. We generally follow the
prescription of Tacchella et al. (2022) for our Prospector model. We test both
a traditional parametric star formation history as well as the “continuity
bursty” non-parametric SFH used in Bagpipes. We allow the stellar mass to vary
between 6 $\leq\log_{10}M/M_{\odot}\leq$ 12 with a uniform logarithmic prior.
For our parametric star formation history we use the ‘delayed-$\tau$’ model,
where $\textrm{SFR}(t)\propto(t-t_{s})e^{-(t-t_{s})\tau}$.

We note that by default Prospector provides the total stellar mass formed,
rather than the surviving stellar mass. Post-fitting we recalculate the return
fraction within Prospector for the full posterior in order to derive a
surviving stellar mass distribution for each galaxy. We allow the V-band
optical attenuation due to dust to vary between 0 and 6 magnitudes, with a
uniform prior, assuming a Calzetti et al. (2000) law also used in our
‘fiducial’ Bagpipes run.

We model IGM attenuation following Madau (1995). Following Tacchella et al.
(2022) and due to possible line-of sight variations in the optical IGM
attenuation, we allow the IGM opacity to vary with a clipped Gaussian prior
distribution centered on 1, clipped at 0 and 2, and with a dispersion of 0.3.
The stellar metallicity is allowed to vary between -4
$\leq\log_{10}Z/\mathrm{Z}_{\odot}\leq$ 0.16 with a uniform prior. We do not
link the gas-phase metallicity, which is also free to vary with the same range
and prior. Gas-phase metallicity is expected to differ from stellar
metallicity at some stages of galaxy evolution (e.g. following significant gas
accretion into a galaxy, lowering gas-phase metallicity) and decoupling them
permits more flexibility in the stellar and nebular emission line modelling.

We use a Kroupa (2001) IMF for consistency with our Bagpipes SED fitting, but
as Prospector allows tabulated IMFs we also test the impact of a top-heavy IMF
on the derived stellar masses. Top-heavy IMFs are predicted at high-$z$, and
there is preliminary evidence for a top-heavy IMF in extreme nebular dominated
galaxies at high-$z$ (Cameron et al., 2023). The results for different choices
of IMF are detailed in § 4.5.

We run Prospector using nested sampling with dynesty (Speagle, 2020), using
the default sampling settings with the exception of switching from uniform
sampling to the more robust random walk, to efficiently probe the multi-
dimensional parameter space. Since Prospector is computationally expensive to
run, we fit only a subset of our galaxy sample, prioritising those galaxies
which are high-redshift or massive. Specifically we fit all galaxies with
$z_{\textrm{phot}}$ (from EAZY-py) above $z\geq 8.5$, or with a fiducial
Bagpipes stellar mass of $\log_{10}(M_{\star}/M_{\odot})\geq 9.0$.

###  4.3 Stellar Mass Comparisons

Figure 1: Top) Comparison of the observed stellar mass offset between our
fiducial Bagpipes model and alternative models for the entire EPOCHS v1
sample, as a function of fiducial stellar mass. The alternative models
considered vary priors on the metallicity and dust attenuation, as well as the
assumed dust law. Marker shapes show the redshift bin for each galaxy, based
on the fiducial redshift, and each Bagpipes model considered is shown in a
different color. The colored lines show the Locally Weighted Scatterplot
Smoothing fit (LOWESS, Cleveland, 1979). The right-hand plots show the stellar
mass offset as a function of $\Delta\chi^{2}$, between the two fits. Bottom)
Same as the upper plot, but the alternative models vary the assumed star
formation history or stellar population synthesis model instead.

Figure 1 shows a comparison of the derived galaxy properties of our fiducial
Bagpipes run for all sources in our sample to our other Bagpipes results. The
details of the models compared are shown in Table 2. As discussed in § 4.1, we
compare models with different priors, SPS models and SFH parametrizations. The
top plot shows the systematic mass offset found when varying the dust law and
priors, while the second row shows the same offset for parametric and non-
parametric SFHs, as well as a BPASS SPS model (Stanway & Eldridge, 2018). For
each Bagpipes model we show the Locally Weighted Scatterplot Smoothing fit
(LOWESS, Cleveland, 1979). The LOWESS estimator is a non-parametric fitter for
noisy data, and shows the overall trend between the mass discrepancy and
fiducial stellar mass.

We allow the redshift to vary within the EAZY-py posterior PDF for each
Bagpipes fit, resulting in only small changes in redshift between individual
Bagpipes results. 99.6% of photo-$z$ estimates fall within 15% of our fiducial
Bagpipes redshift.

In the following sections we compare the consistency of the derived stellar
masses between our fiducial and alternative Bagpipes models, as shown in
Figure 1. We also include some limited discussion of other inferred galaxy
properties, such as star formation rate, or dust content where they vary
significantly from our fiducial model.

In Figure 3 and Figure 4 we show examples of the photometry, best-fitting
Bagpipes & Prospector SEDs, and posterior redshift and stellar mass estimates
for a selection of high-mass and or discrepant galaxies in our sample.

####  4.3.1 Impact of Priors

The first comparison we make is a substitution of the default logarithmic
prior on the V-band dust extinction, A${}_{\textrm{V}}$, to a uniform prior. A
uniform prior favors higher dust extinctions, and the largest difference to
the fiducial model is seen in galaxies with A${}_{\textrm{V, fid}}<0.2$ mag,
where the dust content is poorly constrained. We see good agreement in
photo-$z$, along with a few extreme outliers, which is unsurprising given the
informative redshift prior from EAZY-py we use. The stellar mass offset is
shown in light blue in the top plot of Figure 1 and shows good agreement
within the posterior uncertainties. Galaxies with higher levels of dust are
found to be significantly younger and more star-forming with a uniform dust
prior, suggesting the inferred star formation histories are dependent on the
dust prior. Comparison of the best-fitting $\chi^{2}$ values shows that both
models are equally well-fitted to the photometry.

We also test the impact of our metallicity prior, which is logarithmic in our
fiducial Bagpipes model. This favours low metallicity, which we expect in the
early Universe. Here we exchange this prior for a uniform distribution which
favours higher metallicity, shown in orange in the top plot of Figure 1. We
see little overall impact on the stellar mass from the metallicity prior, with
individual galaxies scattering up to 0.5 dex and the majority consistent with
results from our fiducial Bagpipes run.

Overall the impact of the dust and metallicity priors alone appears to have
only a small systematic effect on the derived galaxy masses. However, in a
small number of individual cases masses can scatter by $\sim$0.5 dex with
little difference in the goodness of fit.

####  4.3.2 The Assumed Dust Law

Our fiducial Bagpipes model assumes a simple one-component Calzetti et al.
(2000) dust law. The slope of the dust law is known to vary in some galaxies,
and numerous alternative models (e.g. Charlot & Fall, 2000; Salim et al.,
2018) including additional parameters have been suggested. Charlot & Fall
(2000) fit a two-component dust law, with different amounts of extinction for
young ($\leq 10$ Myr) and old stellar populations, to account for dust in
stellar birth clouds. The resultant deviations in stellar mass are shown in
magenta in the top plot of Figure 1, with a significant deviation at the
highest stellar masses, which can be seen in the LOWESS trend. For galaxies
with a fiducial stellar mass of $\geq 10^{10}$ M⊙, the majority are found to
have $\geq$1 dex larger stellar masses with the Charlot & Fall (2000) dust
law. However for these galaxies the goodness of fit is considerably poorer,
with $\Delta\chi^{2}\geq 10$ with the Charlot & Fall (2000) model.

Despite the same prior, the dust attenuation (A${}_{\textrm{V}}$) posteriors
are quite different between the two models, due to the degeneracy between
A${}_{\textrm{V}}$ and slope. Galaxies with moderate dust extinction in the
fiducial results ($\leq 1$ mag) are typically found to have very little dust
extinction on the old stellar population, but the most dust-obscured galaxies
(A${}_{\textrm{V}}>2$) have significantly greater extinction (1.5$\times$ more
on average). The majority of galaxies favor a steeper slope ($n$) for the
attenuation power-law than given by Calzetti et al. (2000), in this mode with
a distribution centered on n$\approx 1.3$, which is 2$\sigma$ from the prior
value of 0.70.

Salim et al. (2018) allow a deviation in slope compared to the Calzetti et al.
(2000) and an additional UV bump at $2175~{}\mathrm{\AA}$. The UV bump in the
Salim et al. (2018) model is driven primarily by Polycyclic Aromatic
Hydrocarbon (PAH) emission, which is expected when emission is strongly
reprocessed. We replace the Calzetti et al. (2000) dust law with the Salim et
al. (2018) model, and fit for these additional components, as shown in Table
2. We find that the inferred stellar mass can increase by $>$1 dex in some
cases, particularly at the highest fiducial stellar mass, but the average
offset (see the LOWESS fit in Figure 1 is smaller than for the Charlot & Fall
(2000) case. In a number of cases the goodness of fit for the Salim et al.
(2018) model is considerably better than our fiducial model, with
$\delta\chi^{2}\leq-10$ in a number of galaxies.

####  4.3.3 Comparison to ‘delayed’ SFH

To test the consistency of galaxy stellar mass estimates with different
parametric star formation histories, we replace our ‘lognormal’ SFH with
another commonly used SFH; a delayed-$\tau$ SFH. We find systematically
slightly lower stellar masses above 108 M⊙, but reasonable agreement at the
lowest stellar masses. This is likely due to the log-uniform prior on age used
in the delayed-$\tau$ model, which results in a younger stellar population.
Comparison of the goodness of fit via the $\chi^{2}$ parameter suggests the
models typically have slightly poorer fits than the fiducial model, but in
some cases the stellar masses are reduced by $\sim$ 0.5 dex with very little
impact on the goodness of fit.

####  4.3.4 Comparison to “continuity bursty” SFH

Here we replace our fiducial Bagpipes ‘lognormal’ SFH with a “continuity
bursty” SFH model described in detail in § 4.1. We reproduce the result of
Tacchella et al. (2022) that this SFH increases the galaxy stellar mass,
finding an average increase of 0.2 dex. The stellar mass discrepancy between
the two models is shown in purple in the lower figure of 1. In individual
cases the increase in stellar mass can reach $\approx$1 dex, with only small
changes in photo-$z$. The largest offsets in stellar mass are typically seen
for galaxies with a fiducial stellar mass of $10^{8}$ M⊙, with the highest
fiducial stellar mass galaxies ($\geq 10^{10}$ M⊙) seeing considerably smaller
increases.

We see the largest discrepancies in stellar mass between the two models when
the $\chi^{2}$ significantly favours the fiducial model, suggesting the higher
stellar mass estimate of the “continuity bursty” model may not be accurate in
these cases. However there are a small number of galaxies in which the offset
exceeds 1 dex with minimal $\chi^{2}$ difference, and even a few galaxies
where the offset exceeds 0.5 dex and the “continuity bursty” model is
significantly preferred. An example of the discrepancy between the stellar
mass PDFs of our fiducial and the “continuity bursty” model can be seen for
the individual galaxy SED of JADES-Deep-GS:9075 shown in Figure 4, where there
is a 1.3 dex difference between the stellar mass estimates, with only a
difference of $\Delta\chi^{2}$ between the two solutions.

Star-formation rates (SFRs) are typically higher, however the model does not
reproduce the highest SFR estimates of the fiducial model. Individual
inspection of these highly-star forming galaxies show that the “continuity
bursty” model struggles to reproduce SFRs high enough to match the measured
H$\beta$+[OIII] equivalent widths inferred from the photometry. Despite the
greater flexibility (and number of fitted parameters), the “continuity bursty”
model has a higher $\chi^{2}$ than our fiducial model for the majority of
galaxies. This may suggest that the fitting procedure is struggling to
accurately constrain the SFH in the non-parametric case.

####  4.3.5 Comparison to BPASS SPS Models

By default Bagpipes uses the 2016 version of Bruzual & Charlot (2003)’s (BC03)
stellar population synthesis (SPS) models. However Bagpipes can also employ
Binary Population and Spectral Synthesis (BPASS, Stanway & Eldridge, 2018)
models. Specifically we use models generated with v2.2.1 of BPASS, with the
default IMF (slope of 1.35, 300 M⋆). The IMF parametrization of this SPS model
differs slightly from the Kroupa (2001) IMF used in the BC03 models, causing
an intrinsic offset in stellar mass. However as shown in orange in the lower
plot of Figure 1, the comparison of the mass estimates is more complex. At the
lowest and highest stellar masses ($\log_{10}M_{\star}/M_{\odot}\leq 7.5$ and
$\log_{10}M_{\star}/M_{\odot}\geq 10$), we see systematically lower stellar
masses than the fiducial Bagpipes model, which is what we would expect based
the IMF difference alone. However for the majority of the sample, which falls
between these two mass extremes, we see a significantly larger offset in
stellar masses, which increases with redshift, as shown by the LOWESS trend.

The SFRs and redshifts are broadly correlated, with some outliers. However the
ages show a large scatter, and higher mass-weighted ages (MWA) are preferred
in the majority of cases when using BPASS, which may be more reasonable that
the young ages ($\langle$MWA$\rangle\sim 20$Myr) inferred with our fiducial
Bagpipes results. In particular, there are a significant subset of galaxies
found to have a $\geq$2 dex shift in mass-weighted age.

The best-fit $\chi^{2}$ shows some scatter, where fits with low $\chi^{2}$ in
the fiducial Bagpipes results typically having a similar or worse fit, but
some galaxies with higher $\chi^{2}$ having significantly improved. For these
galaxies with improved $\chi^{2}$ ($\Delta\chi^{2}>5$), the main difference is
that the best-fitting SED reproduces the observed rest-UV fluxes more closely.
Interestingly, this subset of galaxies with significantly improved $\chi^{2}$
also typically have considerably higher MWA using BPASS than BC03, with their
recreated SFHs suggesting a constant SFH, rather than the recent burst
preferred when fitting with the BC03 SPS models.

####  4.3.6 Other Comparisons

Numerous other Bagpipes models were tested, and we have presented in detail
the results of a subset of them above. Here we summarise the effects of a few
other variations which we do not include in Figure 1.

1\. Fixed redshift to EAZY-py max $P(z)$: Little overall effect on stellar
masses or star formation rates, with some individual scatter. This is the
expected behaviour as the majority of redshifts are consistent within the
EAZY-py posterior uncertainty, given that we use this as a prior in our
fiducial model.

2\. “continuity” non-parametric SFH model’: This variation tests the other
non-parametric SFH model introduced in Leja et al. (2019) and Tacchella et al.
(2022), which more tightly constrains the SFR in neighbouring time bins to
force more smoothly varying star formation histories than the “continuity
bursty” model. We find overall similar behaviour to the “continuity bursty”
model, with systematically higher masses and higher mass-weighted age. The
models typically have worse $\chi^{2}$ statistics than the “continuity bursty”
model fits potentially indicating that more stochastic star formation
histories are preferred for the majority of galaxies, as suggested in the
literature (Faucher-Giguère, 2018; Looser et al., 2023; Asada et al., 2023;
Endsley et al., 2023a).

###  4.4 Comparison to Prospector

In this section we compare our fiducial Bagpipes model to our results from
Prospector. Figure 2 shows a comparison of the stellar masses derived from
Bagpipes and Prospector. We compare our ‘fiducial’ Bagpipes model, as
described in § 4.1, to our parametric SFH Prospector model. These models are
generally similar with a few key differences, but both provide a baseline for
comparison to our other models. We also compare our non-parametric SFH models
to each other, which both employ the same “continuity bursty” SFH, with the
same time bins and priors. We use the same Kroupa (2001) IMF for both Bagpipes
and Prospector, so we do not expect any different in stellar mass estimates
due to the IMF parametrization.

Figure 2: (top) Comparison of redshift estimates between Bagpipes and
Prospector, both with the same Gaussian photo-$z$ prior from EAZY-py.
Photo-$z$ estimates are systematically larger in Prospector when comparing the
parametric SFH models. Median photo-$z$ offset for each model is shown on the
plot. (bottom) Comparison of derived stellar mass estimates between Prospector
and Bagpipes for both parametric and non-parametric SFH models. Median stellar
mass offset is shown on the plot, as is the average uncertainty for both
stellar mass estimates.

For the comparison of the parametric SFH models (delayed-$\tau$ for
Prospector, lognormal for Bagpipes), we see systematically larger photo-$z$
and stellar mass estimates with Prospector. Whilst in the Prospector model we
allow the IGM attenuation to vary, which could impact photo-$z$ estimates, we
do not see the same offset with the non-parametric SFH model, where the IGM
attenuation is also allowed to vary. Both SED-fitting tools are given the same
redshift prior from EAZY-py. We can see for both example galaxy SEDs in Figure
4 that Prospector instead prefers higher-$z$ solutions than Bagpipes because
it is inferring the presence of Lyman-$\alpha$ emission. Whilst Lyman-$\alpha$
emitters have been found at $z\geq 9$ (e.g. Bunker et al., 2023b), we do not
expect to observe Lyman-$\alpha$ from the majority of galaxies at these
redshifts due to the attenuation from neutral hydrogen during the Epoch of
Reionization, and a systematic photo-$z$ offset from Prospector is likely for
galaxies without Lyman-$\alpha$ emission. However given that we see a
photo$-z$ offset primarily for the ’parametric SFH’ only, Lyman-$\alpha$
emission is unlikely to be the only cause of the offset as the Prospector SEDs
in Figure 4 show Lyman-$\alpha$ emission for both SFH models. This offset will
generally cause a slight increase in stellar mass estimates, as a more distant
galaxy must be intrinsically brighter. We see a median increase of 0.43 dex in
stellar mass.

For the “continuity bursty” SFH models we see better agreement in redshift,
with a median offset of only $\delta z=$ 0.1. Stellar mass estimates are also
more consistent on average, although we see large scatter, with individual
mass differences reaching $\approx$1.5 dex.

We typically see comparable $\chi^{2}$ for both the parametric and non-
parametric Prospector SFH models, in contrast to the result with Bagpipes. It
is possible that the nested sampling with dynesty in Prospector provides a
more robust constraint on the binned SFH than the nested sampling in Bagpipes,
and may warrant further investigation.

Figure 3: Bagpipes (labelled BP) and Prospector SEDs (labelled PR) for two of
the galaxies in the EPOCHS v1 sample. The top row shows the galaxy labelled
7463_CEERSP1 with one of the largest fiducial stellar masses in the sample,
which is a ’Little Red Dot’ also identified by Labbé et al. (2023); Kokorev et
al. (2024). The bottom row shows 11955_NEP-4, one of the more massive galaxies
found at $z\geq 10$, although the range of stellar mass estimates shows how
difficult it is to constrain stellar mass at these redshifts. We label the
different SED fits by how they differ to our fiducial Bagpipes SED-fitting
prescription, with a lognormal SFH, Calzetti et al. (2000) dust model,
logarithmic priors on age, dust and metallicity, Bruzual & Charlot (2003) SPS
model and a Kroupa (2001) IMF. For each alternative model we also give the
redshift, stellar mass and $\chi^{2}$ for the fit. On the right of each SED
plot we show the posterior redshift and stellar mass PDF distributions and the
best-fit SFH history is shown as an inset axes on the upper left.  Figure 4:
The same as Figure 3 except these galaxies are shown to demonstrate the
stellar mass discrepancies observed between different Bagpipes and Prospector
SED fits, depending on the chosen star formation history, SED fitting tool,
dust law, and IMF. JADES-Deep-GS:9075 was observed with NIRSpec as part of
Bunker et al. (2023a) (NIRSpec ID 00002773), finding $z_{\rm spec}=12.63$,
meaning that Bagpipes and Prospector slightly underestimate the redshift in
most cases.

###  4.5 Top-heavy IMF

As discussed in § 1, one possibility to explain the high inferred stellar
masses of high-$z$ galaxies is a ’top heavy’ IMF, which results in the
production of more high mass stars compared to a local IMF, and lowers the
inferred mass-to-light ratios in high-$z$ galaxies. We use Prospector to
investigate the impact of varying the initial mass function on the inferred
stellar mass of a subset of the EPOCHS v1 sample. We implement the modified
Kroupa (2001) IMF suggested by Steinhardt et al. (2023), which assumes a gas
temperature evolution T${}_{gas}\propto(1+z)^{2}$. We produce two modified
IMFs, one with T${}_{gas}=45$K, which we use for $8\leq z\leq 12$ , and one
with T${}_{gas}=60$K, which we use at $z\geq 12$. A standard Kroupa (2001) IMF
would have T${}_{gas}=20$K in this parametrization, and given the broken-
power-law shape this modification results in an increasingly top-heavy IMF
with increasing gas temperature. Whilst it is theoretically possible to
produce a unique IMF for each galaxy by assuming a $z-T_{gas}$ relationship,
we avoid doing this for simplicity. We otherwise leave unchanged the
Prospector configuration in order to directly compare the impact of the IMF.

Figure 5: (top) Effect of modified top-heavy IMF on galaxy masses using a
parametric delayed-$\tau$ SFH, for both the HOT 45K IMF (at $8<z<12$) and HOT
60K IMF (at $z>12$) from Steinhardt et al. (2023). Overlaid on the plot are
the median stellar mass and $\chi^{2}$ offsets, showing that these IMF models
systematically reduce the stellar mass with no impact on the goodness of fit.
(bottom) Effect of modified IMF on galaxy masses using a non-parametric
“continuity bursty” SFH for both HOT models.

In Figure 5 we show the results of modifying the IMF on the stellar mass
estimates. We compare the masses derived with a standard Kroupa (2001) IMF to
the two modified IMF models, finding significant decreases in stellar mass
with very little change in the quality of the fits. Noting that due to the
computational intensity of Prospector fitting we only fit a subset of our full
sample, and calculate the median decrease in mass for both SFH models and HOT
IMFs. For the $z>12$ sample where we use the HOT 60K IMF we include our full
parent sample, however for the $z<12$ model we fit only galaxies with either
$z\geq 8.5$ or fiducial Bagpipes stellar mass
$\log_{10}M_{\star,\mathrm{fid}}>9.5$. In terms of numbers, there are 221
galaxies in the HOT 45K IMF group, and 21 galaxies in the HOT IMF 60K group.
Given that this IMF model is predicted to be applicable in the region of 8
$\leq z\leq$12 (Steinhardt et al., 2023), and the lack of mass-dependence in
our results, we do not expect the excluded lower mass galaxies at $8\leq z\leq
8.5$ to significantly impact our findings.

For the non-parametric “delayed” SFH model, we find a median decrease in
stellar mass of 0.33 dex for the $z\leq 12$ (HOT 45K IMF) galaxy sample. For
the $z\geq 12$ sample, we find a median decrease of 0.46 dex. For the
“continuity bursty” SFH model, we find a median stellar mass decrease of 0.26
dex for the $z\leq$12 sample (HOT 45K IMF) and a decrease of 0.36 dex for the
$z\geq$12 sample (HOT 60K IMF). As expected, we see a larger decrease in
stellar mass for the HOT 60K model, which is more ‘top-heavy’ than the HOT 45K
model.

For the galaxies with the largest fiducial Bagpipes stellar masses, which have
the most tension with $\Lambda$CDM, we show the possible decrease in stellar
mass using a top-heavy HOT IMF in Figure 12.

##  5 Galaxy Stellar Mass Functions

The Bagpipes and Prospector fitting we perform in § 4.3 explores the
consistency of stellar mass estimates on an individual galaxy basis. In order
to explore the overall effect of a particular choice of SED fitting tool and
model we look at the overall distribution of galaxy masses via the Galaxy
Stellar Mass Function (GSMF). We focus primarily on constructing a GSMF from
our fiducial Bagpipes results, and then demonstrate the effect of changing
this to an alternative Bagpipes or Prospector model.

We make a further cut to our galaxy sample here, removing the fields of El
Gordo, SMACS-0723 and Clio. The depths and available filters of these datasets
mean that they do not contribute significant volume to our GSMF estimates, but
increase the redshift uncertainties (Adams et al., 2023). This reduces the
number of galaxies used in the GSMF to 1092, and the total area used to 175
$\arcmin^{2}$.

The GSMF measures the abundance of galaxies of different masses at a given
redshift. A stellar mass function $\Phi(M,z)\Delta M$ is formally defined as
the number density of galaxies in a mass bin $\delta M$ at a given redshift
$z$. The evolution of the shape and normalisation of the stellar mass function
traces the global abundance of baryons across cosmic time, and hence
indirectly traces star formation. The integral of the stellar mass function
over mass gives the galaxy stellar mass density, which is the cumulative
formed stellar mass per unit volume at a given epoch.

We construct a galaxy stellar mass function (GSMF) from different mass
estimates in order to investigate possible evolutions of the GSMF at high
redshift. To derive the GSMF we use the 1/V${}_{\textrm{max}}$ methodology
(Rowan-Robinson & McCrea, 1968; Schmidt, 1968):

$\phi(M)d\log(M)=\frac{1}{\delta M}\sum_{i}\frac{1}{C_{i}V_{max,i}}$ (1)

with associated uncertainty

$d\phi(M)=\frac{1}{\delta
M}\sqrt{\left(\sum_{i}\frac{1}{C_{i}V_{max,i}}\right)^{2}}$ (2)

where $\delta M$ is the bin width in stellar mass, $C_{i}$ is the completeness
of the galaxy in bin $i$, and $V_{max,i}$ is the total observable volume of
the galaxy across all the fields. We use Equation 2 to calculate the
uncertainty in the bin, except in the case where there are the bin has very
low occupancy (N $\leq$ 4) where we instead directly calculate the Poisson
confidence interval with a more accurate estimator based on the $\chi^{2}$
distribution $I=[0.5\chi^{2}_{2N,a/2},0.5\chi^{2}_{2(N+1),1-a/2}]$ (Ulm, 1990;
Adams et al., 2023), which avoids uncertainties such as $1\pm 1$, which appear
infinite on a log-scale.

We iteratively shift the galaxy SED from the fiducial Bagpipes SED fitting in
small steps of $\Delta z=0.01$ before recalculating the bandpass-averaged
fluxes in the available NIRCam and HST filters for a given field. With these
fluxes we test if the galaxy would still be selected at every redshift step
given our selection criteria detailed in § 3.2. For the selection criteria
which are dependent on SNR requirements we base the detection strength on the
average depth for each field, as given in Table 1. This allows us to calculate
a $z_{\mathrm{max}}$ and $z_{\mathrm{min}}$ redshift for each galaxy in each
field, capped at the edges of each redshift bin. Accounting for a minimum
redshift is essential for accurately measuring the detectable volume given our
requirement for the shortest wavelength filter bandpass to fall blueward of
the Lyman break. We convert these maximum and minimum redshifts within each
bin to a volume by

$V_{\mathrm{max},i}=\sum_{\mathrm{fields}}\frac{4\pi}{3}\frac{A_{s}}{4\pi\textrm{
sr}}(d_{z_{\mathrm{max}}}^{3}-d_{z_{\mathrm{min}}}^{3})$ (3)

where $d_{z_{\mathrm{max/min}}}$ are the co-moving distances at the
maximum/minimum detectable redshifts (capped at bin edges). $A_{s}$ are the
survey areas, which are given in Table 1. For the purposes of constructing a
GSMF we exclude the Clio, El Gordo and SMACS-0723 fields due to their shallow
depth or lack of key broadband photometry (e.g. F115W).

The only field containing a significant lensing cluster that we include in the
GSMF is the MACS-0416 field. We exclude the cluster itself, and assume there
is no significant magnification for our galaxy candidates in the surrounding
NIRCam pointings. None of our candidates show evidence of strong lensing, and
we do not expect weak lensing by foreground neighbours to significantly affect
our sample.

We account for the posterior stellar mass uncertainties using a Monte Carlo
bootstrap methodology to bootstrap our GSMF. We draw stellar mass estimates
for each galaxy from the posterior stellar mass probability distribution
functions and compute 1000 independent realisations of the galaxy stellar mass
function from these posterior PDFs. We compute the 16th, 50th and 84th
percentiles of the distribution in order to quantify the uncertainty
introduced by our stellar mass estimates, and indirectly by the uncertainties
on the photometry. Typically these uncertainties are smaller than the Poisson
error, but can be large at the highest mass where individual objects near bin
edges can have a large impact. We require a bin to be occupied in more than
20% of realizations to avoid highly uncertain bins which have very low
occupancy.

###  5.1 Detection Completeness

We carry out completeness simulations on JWST data by inserting simulated
galaxies with an exponential light profile (Sérsic index of $n=1$) as galaxies
at high-$z$ are typically not concentrated (see Ferreira et al., 2022, 2023;
Ormerod et al., 2024) and with absolute UV magnitudes ranging from -16 to -24
in the detection image (inverse variance weighted stack of F277W, F356W and
F444W). We then run SExtractor on them using the same configuration as our
normal catalogue creation pipeline in order to measure the fraction recovered
as a function of apparent magnitude.

The magnitude range was swept in steps of 0.2, with 1000 sources inserted into
the image at a time in order to prevent introducing artificial overdensities.
We assume a size - luminosity relation of $r\propto L^{0.5}$ (Grazian et al.,
2012), with reference size set to $r=800$ pc at $M_{\mathrm{UV}}=-21$, the
maximum available size was on the order of 1000 pc. The inclination angle of
these sources was assumed to be uniformly random. We sweep a redshift range of
6.5 to 13.5 in steps of 0.5 to account for the changes in angular diameter
distance and obtain a 2D dependence of completeness on redshift and absolute
magnitude.

The above procedure was carried out on all fields considered in this paper.
Given that we perform our detection on a stacked image of the LW bandpasses,
our detection completeness is $>90$% across the redshift ranges and apparent
magnitudes of our high-$z$ sample, and the overall completeness is dominated
by the impact of our selection criteria.

###  5.2 Selection Completeness and Contamination

Figure 6: Top) Selection completeness as a function of redshift and stellar
mass for two example fields from our sample, based on our fiducial Bagpipes
SED fitting and full selection procedure. Completeness and contamination are
labelled as fractions, where $1=100\%$. Bottom) Total area-weighted
completeness and contamination for each of our galaxies as a function of
stellar mass, colored by redshift bin. We see higher completeness and lower
contamination at higher mass. The red shaded region in the contamination plot
shows our 50% contamination limit.

To estimate the contamination and completeness of our sample selection as a
function of stellar mass we use version 1.2 of the JADES Extragalactic Ultra-
deep Artificial Realization (JAGUAR) catalogue of synthetic galaxies released
by the JADES Collaboration (Williams et al., 2018). We choose this model for
our completeness simulations as the catalogues and spectra are readily
available, and Wilkins et al. (2023a) have shown that the JAGUAR model
reproduces the observed color evolution of galaxies at high redshift.

JAGUAR is a empirical model based on observational constraints on the mass and
luminosity functions at $z\leq 10$ from HST (Tomczak et al., 2014; Bouwens et
al., 2015; Stefanon et al., 2017). Mock spectra and photometry are generated
using BEAGLE (Chevallard & Charlot, 2016). We combine 5 different realisations
of the simulation in order to improve the statistics for rare high-$z$ and
high-mass sources. We filter the mock catalogues to the sources which are
detectable given the average depths of our observations. This ranges between
$8\times 10^{5}$ and $1.2\times 10^{6}$ realistic mock galaxies, which we run
through our full galaxy selection procedure.

We generate mock observations from these catalogues by estimating 1$\sigma$
errors in the measurements based on the average depths of our JWST reductions
in each field (see Table 1). We then perturb the catalogue photometry in flux
space within a Gaussian centered on the catalogue measurement and width equal
to the 1$\sigma$ error in that filter. We then run the galaxies through our
full selection procedure to measure how well our pipeline recovers the true
redshifts of the sources. This allows us to compare the robustness of our
selection criteria given the differing depths and filters in each set of
observations.

We test our ability to reconstruct the stellar masses of the JAGUAR galaxies
by fitting their perturbed photometry with our Bayesian SED-fitting using
Bagpipes. We estimate completeness and contamination for each redshift and
mass bin in our mass function by testing the recovery of simulated galaxies
into the correct bin. Completeness is defined as the number of galaxies which
our pipeline places in the correct redshift and mass bin, divided by the total
number of simulated galaxies within that bin. Contamination is defined as the
number of galaxies we place in the incorrect bin, divided by the total number
of galaxies in the bin. As an example, if there are 100 simulated galaxies in
a given redshift/mass bin, and our selection criteria selects 50, of which 10
are actually from a different bin, then the completeness would be 50% (50/100)
and the contamination would be 20% (10/50).

We see reasonable agreement between the true stellar mass and our inferred
stellar mass. Initially 50% of galaxies fall within the correct stellar mass
bin used in construction of the galaxy stellar mass function. In order to
account for the different choice of IMF and SFH parametrization between the
JAGUAR SEDs and our Bagpipes fitting, we derive a corrective scaling factor by
which we scale the JAGUAR masses in order to increase the agreement between
the mass estimates to 65%. As we bootstrap the GSMF using the stellar mass
PDFs, we find that 80% of galaxies will contribute to the correct mass bin.
The average offset between our recovered stellar mass and the simulated mass
is 0.02 dex, with a standard deviation of 0.28 dex. Only 2% of the time is the
deviation between the scattered and recovered masses greater than 0.75 dex,
which is the bin width used in the stellar mass functions.

Figure 6 shows the completeness and contamination for the CEERS and JADES
fields, which make up a significant fraction of our total volume.

To derive a single estimate of completeness and contamination for each galaxy,
we calculate the area-averaged completeness and contamination for each galaxy
in our sample, based on the fields in which the galaxy is found to be
observable in when we calculate V${}_{\textrm{max}}$. The bottom row of Figure
6 shows the total completeness, including both detection and selection, for
each galaxy as as function of stellar mass. Completeness is strongly dependent
on the stellar mass of the galaxy, and at masses below $\approx
10^{8}M_{\odot}$ we are $\leq 50\%$ complete in the majority of our area, and
dependent on the derived completeness corrections.

In order to reduce the possible effects of contamination on our derived mass
functions, we exclude objects where the estimated contamination fraction is
$\geq$50%. In the cases where no simulated galaxies fall within a given bin,
which occurs only for the highest mass bins where we expect high completeness,
we do not apply any correction factor. We also assume 100% completeness when
the total number of simulated galaxies in a mass-redshift bin is $\leq$5 due
to the high uncertainty in deriving correction factors with low statistics.

Multiple classes of interlopers with photometry similar to high-redshift
galaxies can contaminate estimates of the GSMF. In this section we detail our
methodology to remove common types of interlopers, including low-$z$ galaxies,
brown dwarfs, and AGN.

####  5.2.1 Low-z Galaxies

The misclassification of low-$z$ galaxies as high-$z$ through catastrophic
errors in photo-$z$ are commonly caused by misidentification of different
spectral features, e.g. confusion of the Lyman and Balmer breaks, or strong
rest-frame optical emission lines contributing to multiple wideband
photometric observations, giving the appearance of continuum emission (e.g.
the $z=16$ candidate in Donnan et al. (2022); Arrabal Haro et al. (2023),
which we do not select with as a robust high-$z$ galaxy with our EAZY-py SED
fitting and selection criteria).

Firstly, we compare our EAZY-py photo-$z$ estimates with spectroscopic
redshifts (spec-$z$’s) in Adams et al. (2023). In summary, we typically find
$>90\%$ accuracy (within 15%). With our EAZY-py photo-$z$’s, which we use as
an informative prior for Bagpipes and Prospector, we do not observe the
systematic photo-$z$ overestimation observed in other studies (Fujimoto et
al., 2023; Arrabal Haro et al., 2023).

We are able to fully test the rate of catastrophic photo-$z$ failure using our
JAGUAR completeness and contamination simulations. The rate of contamination
increases as we approach our 5$\sigma$ detection limits in the rest-frame UV,
and the vast majority of the contaminants have $<8\sigma$ detections in the
rest-frame UV. If we parametrize in terms of stellar mass, the highest rates
of contamination are seen at lower stellar mass where individual galaxies have
less impact on the GSMF. Candidates with SNR close to our selection limit are
also those most commonly rejected in our manual inspection of all candidates,
providing another layer of protection against possible interlopers.

We restrict the redshift range and fields used when constructing the GSMF to
exclude uncertain high-$z$ candidates. We do not include observations of El
Gordo, Clio or SMACS-0723 when constructing the GSMF because they are shallow
and show high rates of contamination in our JAGUAR simulations. Clio and
SMACS-0723 also lack observations in F115W, which Trussler et al. (2023) shows
lead to inaccurate photometric redshifts at $8<z<10$. Despite our sample
containing galaxies at $z\geq 14$, we do not attempt to measure the GSMF above
$z\geq$13.5 because we lack the appropriate medium band observations (or
spectroscopy) to remove interlopers such as the $z=16$ candidate in Donnan et
al. (2022). Below redshift 13.5 we have seen high accuracy in photo-$z$
estimates when compared to NIRSpec spectroscopy (e.g. Curtis-Lake et al.,
2023; Bunker et al., 2023a; Arrabal Haro et al., 2023; Wang et al., 2023a).

####  5.2.2 Brown Dwarfs

Y and T-type brown dwarfs within the Milky Way can masquerade as high-redshift
galaxies due to the appearance of a Lyman-break like dropout caused by strong
molecular absorption in the optical, coupled with bright near-infrared
emission. When mistaken for high-$z$ galaxies brown dwarfs can have a large
impact on the galaxy stellar mass function due to their high inferred stellar
masses.

We fit the photometry of all of our galaxy candidates with the synthetic
Sonora Bobcat brown dwarf templates (Marley et al., 2021) using a simple
$\chi^{2}$ minimization to scale the flux of the templates to each candidate.
We fit all the provided templates to each candidate, selecting the best-
fitting template for each galaxy. We remove candidates from our sample which
are both better fit by a brown dwarf template than a galaxy template, and show
a PSF-like morphology (50% encircled flux radius in F444W $<$ F444W PSF FWHM).
This removes 59 candidates in total (4% of the total sample). Figure 7 shows
an example of several brown dwarf candidates which would otherwise have been
selected as high-$z$ galaxy candidates.

Figure 7: Two example SEDs from our sample of 59 brown dwarf candidates (in
green) selected with the Sonora Bobcat models, shown along with their best-
fitting high-$z$ galaxy SEDs (in orange). The low redshift galaxy solution
(with an upper redshift limit of $z\leq 6$), which is used in our robust
galaxy selection criteria, is shown in blue. The redshift posteriors from
EAZY-py are shown for both the high and low redshift galaxy models, and are
overlaid with the statistics used within our selection criteria, indicating
the $\chi^{2}$ and integral of the primary PDF peak. The cutouts of both brown
dwarf candidates are shown below the best-fitting SEDs, demonstrating their
compact PSF-like morphology. The cutouts are 0$\farcs$9 across, as shown with
the scalebar and the white circle shows the extraction aperture used for SED-
fitting (0$\farcs$16 radius) The scalebar showing physical size in kpc is
calculated at the best-fitting galaxy photo-$z$ and is not applicable to the
brown dwarf solution. The coordinates for these brown dwarf candidates are
17h23m12.36s +65d49m38.8s and 17h21m53.00s +65d49m20.82s respectively.

None of the brown dwarf candidates selected in Hainline et al. (2023) are
contained within our robust galaxy sample in the Goods-South or CEERS fields.
As a check of our selection criteria we confirm that all of Hainline et al.
(2023)’s candidates located within the footprint of our datasets (GOODS-South,
CEERS) are selected as brown dwarfs by our brown dwarf selection criteria.

####  5.2.3 Active Galactic Nuclei

Given the apparent prevalence of high-$z$ AGN detected by multi-wavelength and
spectroscopic studies it is likely there is a population of galaxies
containing AGN within this sample. Our modelling does not fit AGN components
to the photometry, so if galaxies contain a significant AGN component it is
likely we will overestimate the stellar mass. In particular, obscured AGN with
red colors (e.g. Kocevski et al., 2023; Labbé et al., 2023) can mimic strong
Balmer breaks, leading to inference of an aged stellar population and high
stellar mass. Known as LRDs, a number of studies have found numerous
candidates up to $z\leq 9$ (see e.g. Greene et al., 2023; Furtak et al., 2023,
2023; Barro et al., 2023; Matthee et al., 2023; Labbe et al., 2023; Kokorev et
al., 2024). Pérez-González et al. (2024) uses MIRI observations of LRDs to
suggest that the characteristic photometry is driven mostly by stellar
emission, and that LRDs are mostly compact and highly obscured starburst
galaxies with young stellar ages. Other studies have struggled to distinguish
between AGN and stellar models (Barro et al., 2023) or favour an AGN model
(Noboriguchi et al., 2023). Further spectroscopic or MIRI observations can
break the degeneracy between AGN and old stellar populations. We attempt to
perform additional SED fitting using Prospector by including an AGN component,
which models AGN emission lines for the most massive galaxies in our sample
(M${}_{\star}\geq 10^{9.5}\ \textrm{M}_{\star}$). In almost all cases the AGN
fraction is less than 10%, and we see no reduction in stellar mass due to the
inclusion of AGN component. The addition of an AGN model does not
statistically improve the fit for these galaxies, although we can not strongly
exclude the presence of AGN in some of our candidates as the AGN model within
Prospector models only obscured AGN (Leja et al., 2018). Further fitting with
a wider range of AGN models, which is outside the scope of this paper, is
required to constrain the potential contribution of AGN emission to the
observed photometry. Given the uncertainty around the true nature of these
sources, we choose not to remove them from our sample. We discuss the impact
of LRDs on the GSMF in Appendix A.

We do however remove two galaxies from our sample which are spectroscopically
confirmed as AGN, or with significant AGN components within the photometry.
This includes the AGN within the CEERS field presented in Larson et al.
(2023).

###  5.3 Cosmic Variance

Cosmic variance is the field-to-field variance of the distribution of galaxies
at a given redshift due to galaxy clustering. Empirical measurements of cosmic
variance do not exist at $z\geq 7$. Estimates of cosmic variance from
simulations (e.g. _Bluetides_ , Bhowmick et al., 2020) are available, but the
accuracy of these estimates is difficult to measure given the limited
cosmological volume and dependence on the assumed cosmology and galaxy
formation physics. Bhowmick et al. (2020) argues that cosmic variance will be
a dominant source of uncertainty at these redshifts, as galaxies are predicted
to be highly clustered (Bhowmick et al., 2018).

We caution that the following methodology provides only an approximation of
the true cosmic variance as the underlying quantities, such as the galaxy
bias, are not well-quantified at $z\geq 7$ and are generally extrapolated from
wider area clustering results at lower redshift or taken from N-body
simulations (Moster et al., 2010). We estimate the uncertainty on the GSMF due
to cosmic variance following the prescription of Driver & Robotham (2010),
combining the cosmic variance from different surveys using the volume weighted
sum of squares from Moster et al. (2011). We treat nearby pointings as one
observation and sum their areas (e.g. the multiple pointings of CEERS or the
NEP, which are not widely separated enough on the sky to be considered
independent). We add the uncertainty in quadrature to Equation 2.

We compute a cosmic variance estimate for each galaxy individually based on
which of our fields they are found to be detectable in. Given that this work
combines multiple widely separated fields for the majority of our sample
cosmic variance is small compared to studies incorporating only a single
field. We see little difference in cosmic variance across our redshift bins.
In the best case, where a galaxy is detectable across every field, the average
cosmic variance across all redshift bins is found to be 18%. The highest
cosmic variance estimates are for faint galaxies which, given our selection
criteria, are only detectable in the JADES GOODS-South field, where we find a
cosmic variance of 42%.

###  5.4 Fitting the GSMF

The overall scientific consensus is that the stellar mass function at high-$z$
is well-described by the Schechter function (Schechter, 1976), given in
logarithmic form in Equation 4.

$\displaystyle\Phi(\mathrm{d}\log M)=\ln(10)\Phi^{\star}\mathrm{e}^{-10^{\log
M-\log M^{\star}}}\times$ (4) $\displaystyle\left(10^{\log M-\log
M^{\star}}\right)^{\alpha+1}\mathrm{~{}d}\log M$

The Schechter function, which consists of a power-law with an exponential
fall-off, is parameterised by 3 parameters: $\alpha$, M⋆, and $\Phi^{\star}$.
$\alpha$ controls the slope of the low-mass end of the SMF, while M⋆ gives the
turnover point at which the function turns from a power-law to an exponential
fall-off. $\Phi^{\star}$ controls the overall scaling of the SMF, and is the
only parameter that has been confidently shown to evolve with cosmic time
(Popesso et al., 2023). Stellar mass functions are sometimes split into star-
forming and quiescent, with a division in specific star formation rate (sSFR).
SMFs for quiescent galaxies are often parameterised by double Schechter
functions (Weigel et al., 2016).

Before JWST the mass function was constrained up to $z\approx 9$ using HST
WFC3/IR, Spitzer/IRAC and ground-based data (Visible and Infrared Survey
Telescope for Astronomy (VISTA), United Kingdom Infrared Telescope (UKIRT))
observations of the HUDF, CANDELS, Cluster Lensing and Supernova survey with
Hubble (CLASH) and Hubble Frontier Fields (HFF) surveys (Stark et al., 2013;
Bradač et al., 2014; Duncan et al., 2014; Bouwens et al., 2016; Song et al.,
2016; Stefanon et al., 2017; Bhatawdekar et al., 2019; Kikuchihara et al.,
2020). These studies generally agree that the mass function steepens out to
$z\approx 10$, with the low-mass slope, $\alpha$, $\rightarrow-2$. Duncan et
al. (2014); Bhatawdekar et al. (2019); Stefanon et al. (2021) find that the
stellar mass to halo mass ratio in galaxies (M⋆-Mh) shows no evolution over
$6\leq z\leq 10$ despite a 3 dex increase in overall stellar mass, suggesting
that stellar and halo mass grew together during reionization. Pre-JWST, the
accuracy of stellar mass functions at these redshifts is limited by the
uncertainty in stellar mass estimates from Spitzer data, where the broadband
photometry can be affected by strong nebular emission causing scatter in
stellar mass estimates, which is compounded by the limited number of
spectroscopically confirmed galaxies at these redshifts (Roberts-Borsani et
al., 2016; Strait et al., 2020). Previously, measurements of the mass function
above $z\geq 6$ have relied on UV-selected samples of galaxies, with mass
indirectly inferred through a calibrated L${}_{\textrm{UV}}$-M⋆ relationship
(Song et al., 2016; Harikane et al., 2016) rather than direct measurements of
the rest-frame optical wavelengths. Photometric mass-selected samples above
$z\geq$2 all have significant mass uncertainties at higher redshift (Retzlaff
et al., 2010; Laigle et al., 2016; Weaver et al., 2022),. Many studies combine
ground and space-based observations in order to probe the entire mass
function, although combining these datasets robustly is difficult due to
systematics between different surveys, including differing selection
functions, SED fitting techniques, detection bands, and survey depths.

We derive a mass function for 5 redshift bins: 6.5 $\leq z<$ 7.5, 7.5 $\leq
z<$ 8.5, 8.5 $\leq z<$ 9.5, 9.5 $\leq z<$ 11.5 and 11.5 $\leq z<$ 13.5, which
covers a time period of $\sim$500 Myr. We cover a mass range from
$7.75\leq\log_{10}M_{\star}/M_{\odot}\leq 11.5$ in steps of 0.75
$\log_{10}M_{\star}/M_{\odot}$. Our lower mass limit is driven by our
completeness simulations, and the bin width is chosen to ensure adequate
statistics in the majority of bins. Given the quantisation of our completeness
correction to each stellar mass/redshift bin, we ensure our mass bins are wide
enough to ensure our SED fitting places galaxies within the correct mass bin,
so that the completeness corrections are calculated correctly. by ensuring
that we can recover the masses of the simulated galaxies to within the correct
stellar mass bin. We fit the GSMF using Bayesian methods; specifically Markov-
Chain Monte Carlo via the emcee package (Foreman-Mackey et al., 2013).

For our highest redshift bin (11.5 $\leq z<$ 13.5), we do not attempt to fit
the mass function for a number of reasons; firstly given the low number of
available data point points a meaningful fit is not possible without fixing
multiple parameters, and secondly due to the increasing uncertainties in
stellar mass estimates and reliability of candidates at these redshifts.
Furthermore the stellar mass estimates are less reliable when there is less
rest-frame optical photometry available during SED fitting. Our JAGUAR
simulations also suggest that this bin has high contamination at lower mass,
with many galaxies below 108.5 M⊙ having $\geq$50% average contamination
across the fields they are detectable in.

We use wide uniform priors on $\Phi^{\star}$, M⋆ and $\alpha$, with $-8\leq\
\log\Phi^{\star}\leq-2,\ 8.0\leq\ \log M^{\star}\leq 12.5$ and
$-3.5\leq\alpha\leq-1.0$ in our MCMC fitting procedure and run until it has
converged before taking $50,000$ independent draws from the posterior. We
calculate the median posterior, the $\pm 1\sigma$ uncertainty from the
16th/84th percentiles as well as the maximum likelihood draw.

###  5.5 Mass Functions at 6.5 $<z<$ 13.5

Table 3 gives our tabulated GSMF for each redshift bin, as well as the
associated 68% (1$\sigma$) uncertainties. We also list the median number of
galaxies which contribute to each mass bin, as well as the average
completeness and contamination estimates. As discussed in § 5, bootstrapping
the GSMF means that the estimates for completeness, contamination, and
occupation shown are only averages as they will vary in each realization of
the GSMF.

Figure 8 shows our derived GSMF (red circles) and associated Schechter fit for
each redshift bin, with a comparison to relevant literature results. For the
lower redshift bins, which have the most comparisons in the literature, we
split our comparison into observational and simulated results to aid
readability. The shaded regions shown represent the 68% (1$\sigma$)
uncertainty combining the cosmic variance, Poisson error, and bootstrapping
via the mass PDFs. Mass bins with $\leq 50\%$ completeness are plotted without
a black border. We also show GSMF estimates without completeness corrections
applied with red diamonds. We convert the GSMF estimates from the literature
to a Kroupa (2001) where necessary, using the assumptions of Madau & Dickinson
(2014).

The uncertainty regions are likely underestimated; (Wang et al., 2023b) argue
the statistical uncertainties encoded in the posterior distributions used for
bootstrapping the GSMF do not represent the true uncertainty on the stellar
mass estimates due to underlying modelling assumptions, such as the SPS model,
IMF, nebular modelling and assumed SFH. We explore the impact of these
assumptions on the GSMF further in § 6.5.6, showing that alternative GSMF
estimates fall outside the uncertainty region.

Table 4 gives the maximum likelihood and median posterior estimates for our
Schechter function parameters for each redshift bin. Uncertainties correspond
to the 16\84th percentile of each parameter distribution. Figure 9 shows the
contours which correspond to the 68% (1$\sigma$) and 95% (2$\sigma$)
confidence levels.

Figure 9 shows the covariance of the Schechter parameters for each GSMF fit. A
range of parameters are able to fit our GSMF estimates, demonstrating the
highly covariant nature of the Schechter parameters. We note that the
confidence intervals for the turnover mass, M⋆, and the normalization,
$\phi^{\star}$, appear somewhat unconstrained at the high and low ends
respectively. The lack of constraint on M⋆ is a consequence of the wide
stellar mass bins used and consequential low number of overall mass bins.
Wherever we find M${}^{\star}>$11.5, the GSMF is fit by a pure power-law only,
and there is no reason to continually extend the prior outside the range where
we have available data to constrain M⋆.

Figure 10 shows the redshift evolution of these Schechter parameters, with a
comparison to those derived in the literature. It is difficult to measure any
significant redshift evolution, given the large uncertainty and covariance of
the results. We observe some evidence of evolution of the low-mass slope
$\alpha$, which steepens toward higher redshift in all but the highest
redshift bin.

In Figure 10 and Figure 9 we see that often our highest-likelihood Schechter
function does not coincide with the median posterior, despite the model being
fully-converged and a large number of posterior sample drawn.

Table 3: Tabulated GSMF from our fiducial Bagpipes SED-fitting. We give the average number of galaxies in each bin as well as the estimated average completeness and contamination estimates based on our JAGUAR simulations, along with the median redshift (in brackets) and stellar mass for all the objects in a given bin. See § 5.2 for the definition of completeness and contamination used in this work. This table is available for download at https://github.com/tHarvey303/EpochsIV. Redshift Bin | $\log_{{10}}(\frac{{M_{{\star}}}}{{M_{{\odot}}}})$ | Med($\log_{10}M_{\star}$) |  $\Phi\ (10^{-4}\textrm{dex}^{-1}$ Mpc-3) | Comp (%) | Cont (%) | Ngal
---|---|---|---|---|---|---
6.5 $\ <z\ \leq$ 7.5 | $8.125\pm 0.375$ | 8.07 | $49.65^{+10.69}_{-10.29}$ | 31 | 23 | 273
(6.94) | $8.875\pm 0.375$ | 8.75 | $7.37^{+1.79}_{-1.57}$ | 70 | 18 | 103
| $9.625\pm 0.375$ | 9.47 | $1.32^{+0.53}_{-0.42}$ | 79 | 10 | 18
| $10.375\pm 0.375$ | 10.40 | $0.64^{+0.46}_{-0.35}$ | 60 | 0 | 5
| $11.125\pm 0.375$ | 10.81 | $0.04^{+0.09}_{-0.04}$ | 100 | 0 | 1
7.5 $\ <z\ \leq$ 8.5 | $8.125\pm 0.375$ | 8.09 | $16.06^{+3.68}_{-3.48}$ | 35 | 25 | 109
(8.02) | $8.875\pm 0.375$ | 8.73 | $1.64^{+0.60}_{-0.45}$ | 75 | 17 | 30
| $9.625\pm 0.375$ | 9.45 | $0.17^{+0.17}_{-0.08}$ | 91 | 11 | 4
8.5 $\ <z\ \leq$ 9.5 | $8.125\pm 0.375$ | 8.02 | $8.16^{+2.30}_{-2.13}$ | 28 | 30 | 40
(8.86) | $8.875\pm 0.375$ | 8.71 | $0.50^{+0.35}_{-0.20}$ | 70 | 17 | 8
| $9.625\pm 0.375$ | 9.48 | $0.05^{+0.10}_{-0.05}$ | 93 | 18 | 1
9.5 $\ <z\ \leq$ 11.5 | $8.125\pm 0.375$ | 8.10 | $4.53^{+1.12}_{-1.00}$ | 26 | 18 | 42
(10.40) | $8.875\pm 0.375$ | 8.73 | $0.64^{+0.25}_{-0.19}$ | 65 | 18 | 16
| $9.625\pm 0.375$ | 9.42 | $0.03^{+0.10}_{-0.02}$ | 93 | 0 | 1
11.5 $\ <z\ \leq$ 13.5 | $8.125\pm 0.375$ | 8.09 | $2.13^{+0.99}_{-0.90}$ | 34 | 40 | 5
(11.94) | $8.875\pm 0.375$ | 8.70 | $0.22^{+0.22}_{-0.12}$ | 90 | 0 | 5
Figure 8: Galaxy Stellar Mass Functions (red markers) and the best fitting Schechter functions (red solid and dash-dot lines indicating the median and maximum likelihood draws from the fit posterior, and the 16th and 84th percentiles are shown by the red shaded region) for the EPOCHS v1 sample derived from the fiducial Bagpipes masses and photometric redshifts. We do not fit the $z\sim 12.5$ bin due to low statistics, uncertain stellar masses and the potential for high contamination. The bottom-right plot shows a comparison of the best-fitting Schechter functions for all redshift bins. Comparison to the following simulations are shown; _BLUETIDES_ (z$\leq$13, Feng et al. (2016); Wilkins et al. (2017)), _DRAGONS_ ($z\leq$7, Mutch et al. (2016)), _DREaM_ ($z\leq$10, Drakos et al. (2022)), _EAGLE_ ($z\leq$7, Furlong et al. (2015); Schaye et al. (2015)), _FIRE-2_ ($z\leq$10, Ma et al. (2018)), _FLARES_ ($z\leq$15, Lovell et al. (2021); Wilkins et al. (2023b)), _Illustris_ ($z\leq$10, Genel et al. (2014)), _Jaguar_ ($z\leq$10,Williams et al. (2018)), _DELPHI_ ($z\leq 20$, Mauerhofer & Dayal (2023)), _Santa Cruz SAM (GUREFT_ , ($z\leq$17, Yung et al. (2019, 2023)), _Universe Machine_ ($z\leq$10, Behroozi et al. (2019)) and the Feedback Free Burst model of Li et al. (2023) (FFB, $\epsilon_{\rm max}=0.2$, $5\leq z\leq 20$). We also show the SMF upper limit from Li et al. (2023), assuming a maximum star formation efficiency ($\epsilon_{\rm max}$) of unity. Comparisons to observational results from Song et al. (2016), Bhatawdekar et al. (2019) (disc-like galaxies), Kikuchihara et al. (2020), Stefanon et al. (2021), Furtak et al. (2021), Weaver et al. (2022), Navarro-Carrera et al. (2023), Gottumukkala et al. (2023) and Wang et al. (2024). The observational results have been converted to a Kroupa (2001) IMF where necessary. In a minority of cases there is $\Delta z\leq 0.5$ between the redshift of our SMF and literature comparisons.  Table 4: Best-fitting Schechter function parameters and uncertainties derived from fitting the derived GSMF. We give both the median posterior parameter values (with 1$\sigma$ uncertainties derived from the $16^{\mathrm{th}}-84^{\mathrm{th}}$ percentiles), as well as the values corresponding to the maximum likelihood draw from the posterior (given in brackets). This table is available for download at https://github.com/tHarvey303/EpochsIV. Redshift Bin | $\alpha$ | M⋆ | $\log_{10}\phi^{\star}$
---|---|---|---
$6.5<z\leq 7.5$ | $-1.94^{+0.10}_{-0.10}(-1.97)$ | $11.57^{+0.63}_{-0.85}(11.89)$ | $-5.98^{+0.87}_{-0.69}(-6.33)$
$7.5<z\leq 8.5$ | $-2.12^{+0.18}_{-0.17}(-2.07)$ | $10.82^{+1.11}_{-1.07}(10.05)$ | $-6.31^{+1.45}_{-1.36}(-5.24)$
$8.5<z\leq 9.5$ | $-2.26^{+0.31}_{-0.30}(-2.48)$ | $10.47^{+1.15}_{-1.18}(11.81)$ | $-6.63^{+1.78}_{-1.50}(-8.93)$
$9.5<z\leq 11.5$ | $-2.15^{+0.26}_{-0.22}(-2.10)$ | $10.51^{+1.18}_{-1.10}(9.48)$ | $-6.58^{+1.61}_{-1.48}(-5.16)$
Figure 9: Confidence intervals for our best-fitting Schechter parameters for
all fitted redshift bins. Overlaid are the locations of the median posterior
value (filled circles) and maximum likelihood draw (filled crosses). Filled
(shaded) regions show the 68% (95%) confidence levels. Figure 10: Evolution of
best-fitting Schechter function parameters with redshift, with the results of
this study shown by the red circles (median posterior) and red crosses
(maximum likelihood). Comparisons to Duncan et al. (2014), Grazian et al.
(2015), Caputi et al. (2015), Song et al. (2016), Bhatawdekar et al. (2019),
Kikuchihara et al. (2020), Stefanon et al. (2021), Furtak et al. (2021),
Weaver et al. (2022), and Navarro-Carrera et al. (2023) are shown. Markers
showing literature values without black outlines were fixed during fitting. M⋆
values have been adjusted to reflect a Kroupa (2001) IMF.

###  5.6 Stellar Mass Density

The stellar mass density measures the cumulative build-up of stellar mass
across cosmic time. This is given by the integration of the stellar mass
function $\phi(M,z)$ at a given epoch (redshift) between two limiting masses,
as shown in Equation 5. Care must be taken when comparing different estimates
of the stellar mass density (SMD), as the choice of IMF means mass estimates
must be converted; for example, to convert from a Chabrier et al. (2000) IMF
to a Salpeter (1955) IMF masses must be multiplied by a factor of 1.64 (Madau
& Dickinson, 2014).

$\rho_{\star}(M,z)=\int_{M_{\textrm{l}}}^{M_{\textrm{u}}}M\phi(M,z)dM$ (5)

The upper and lower integration limits $M_{\textrm{u}}$ and $M_{\textrm{l}}$
are normally taken to be 10${}^{13}M_{\odot}$ and 10${}^{8}M_{\odot}$ (see
e.g. Davidzon et al., 2017; Bhatawdekar et al., 2019; Weaver et al., 2022).

We perform the integration given in Equation 5, where $\phi(M,z)$ is the best-
fitting Schechter function fit for each redshift bin. We make use of the quad
function from SciPy (Virtanen et al., 2020), integrating between 108 M⊙ and
1013 M⊙. We compare the integration of both the maximum likelihood fit as well
as the median posterior fit. For comparison we also integrate the best-fitting
Schechter functions from the literature in each redshift bin, converting to a
Kroupa (2001) IMF if necessary.

Figure 11 shows our estimate for the stellar mass density at 6.5 $\leq z\leq
11.5$ from integration of the GSMF derived from our fiducial Bagpipes mass
functions. We show comparisons to observational SMD results at $z\geq 5.5$.
Estimates of the stellar mass density at $z\leq 5$ can be found in the
literature Madau & Dickinson (2014); Driver et al. (2018). We also tabulate
our SMD values in Table 5.

Table 5: Stellar mass density results from our fiducial Bagpipes SED fits calculated from the integral of the Schechter function for each redshift bin. Values with uncertainties come from the 16th, 50th and 84th percentiles of the posterior, and values in brackets are for the highest likelihood Schechter function. This table is available for download at https://github.com/tHarvey303/EpochsIV.  Redshift Bin |  $\rho_{\star}$ ($\log_{10}$ M⊙ Mpc-3)
---|---
$6.5<z\leq 7.5$ | 6.36${}^{+0.14}_{-0.17}$(6.42)
$7.5<z\leq 8.5$ | 5.52${}^{+0.14}_{-0.13}$(5.51)
$8.5<z\leq 9.5$ | 5.03${}^{+0.18}_{-0.18}$(5.02)
$9.5<z\leq 11.5$ | 4.93${}^{+0.18}_{-0.15}$(4.85)

We integrate all independent posterior draws in order to propagate our mass
function uncertainties into the stellar mass density. We show comparisons to
predictions and measurements of the cosmic star formation rate density (SFRD,
$\psi$) using Equation 6 from Madau & Dickinson (2014) in order to estimate
the inferred stellar mass density:

Figure 11: Stellar mass density as a function of redshift derived from the
integration of the Bagpipes fiducial GSMF. Markers with black borders show
comparisons to observational SMD results from Oesch et al. (2014); Duncan et
al. (2014); Grazian et al. (2015); Song et al. (2016); Bhatawdekar et al.
(2019); Kikuchihara et al. (2020); Weaver et al. (2022); Stefanon et al.
(2021). Markers with colored borders show the results of integrating the SFRD
predictions of JWST-era studies, which includes Donnan et al. (2022); Santini
et al. (2022); Bouwens et al. (2023); Finkelstein et al. (2023b); Willott et
al. (2023a); Harikane et al. (2023); Pérez-González et al. (2023). Where
necessary we convert literature results to a Kroupa (2001) IMF. Markers may be
shifted by up to $\Delta z=0.1$ for clarity. We also show the theoretical
predictions of the Madau & Haardt (2015) (constant star formation efficiency)
and Oesch et al. (2018) (DM halo evolution) models. The purple shaded area
shows the integral of the SFRD presented by Adams et al. (2023) for our
sample, which is consistent with the SMD derived from our Bagpipes SED
fitting. The red dotted errorbar at z$\sim$10.5 shows the SMD range we find
using different SED fitting tools and models, calculated from the Schechter
fits in Figure 13. We do not show the range of models at every redshift, but
it typically exceeds the uncertainty derived from the fit itself.

$\rho_{\star}(z)=(1-R)\int^{\infty}_{z}\psi\frac{dz}{H(z^{\prime})(1+z^{\prime})}$
(6)

The return fraction $R$ is estimated as 0.423 under the assumptions of Madau &
Dickinson (2014) for a Kroupa (2001) IMF (a well-mixed closed-box model with
constant IMF and metal yield and instantaneous recycling of metals). We
integrate the predictions of Madau & Dickinson (2014), Oesch et al. (2018) and
the UV luminosity density derived results of Adams et al. (2023), as well as
other JWST-era measurements. For Adams et al. (2023) we propagate
uncertainties in the inferred stellar mass density via a Monte Carlo
integration of their SFRD measurements and uncertainties. Star formation is
assumed to begin at $z=20$, but the results are relatively insensitive to the
exact formation redshift as long as it exceeds the redshift $z\geq 13.5$ limit
used in this work.

The inferred stellar mass density is heavily dependent on the best-fitting
Schechter function, and therefore indirectly dependent on the GSMF and stellar
mass estimates. For the alternative GSMF estimates in Figure 13, we find that
the use of the “continuity bursty” SFH, or stellar masses derived with our
Prospector fitting, can increase the stellar mass density at redshift 10.5 by
up to 0.75 dex, which is shown with a dashed errorbar. Similar variations are
possible in the other redshift bins, but for simplicity we do not show every
alternate GSMF and SMD derived for all redshift bins.

##  6 Discussion

We have presented our fiducial GSMF using Bagpipes at redshifts 6.5 $\leq
z\leq 13.5$, finding a steep low mass slope ($\alpha\lesssim-1.95$ and a high
stellar mass exponential cutoff (M${}^{\star}\gtrsim 10.5$) in all redshift
bins. We have also estimated the stellar mass density implied by our results,
finding an apparent flattening of the stellar mass density at $z>9$. In the
following section we will discuss our results in context of theory and
previous work. This includes examining the uncertainties in the GSMF derived
from our other SED fitting results, as well as considering how a modified top-
heavy IMF or contamination by hidden AGN would affect our results.

###  6.1 Massive Galaxies in the Early Universe

Numerous studies have reported an excess of galaxies with high inferred
stellar mass at high-redshift (Labbé et al., 2023; Endsley et al., 2023b;
Akins et al., 2023; Xiao et al., 2023). If galaxies of the inferred masses do
exist in the number densities implied by these studies, this could represent a
challenge to $\Lambda$CDM cosmology or our understanding of high-redshift
galaxy astrophysics, given the available timescale and available gas
reservoirs (Labbé et al., 2023; Lovell et al., 2023; Boylan-Kolchin, 2023).

In order to test whether our derived stellar mass and redshift estimates are
in tension with $\Lambda$CDM we use the Extreme-Value Statistics (EVS)
approach, presented in Lovell et al. (2023) and available
online666github.com/christopherlovell/evstats. We follow the methodology of
Lovell et al. (2023) throughout this section. Given the total available area
of our survey (187 arcmin2) we have computed the halo distribution PDF and
corresponding stellar mass PDF.

A universal baryon fraction of 0.16 is assumed based on cosmological results
(Planck Collaboration et al., 2016), and we compute the EVS limits for a range
of stellar fractions (fraction of baryons which form stars). In Figure 12 we
show the stellar fraction upper limit of unity with a dashed line, as well as
the more realistic predictions for stellar fraction, following a simple
truncated lognormal PDF with $\mu=e^{-2},\sigma=1$, with a dotted line and
shaded contours. Overlaid on Figure 12 we show the SED-fitting derived
redshifts and corrected stellar masses for our sample. We correct our derived
stellar masses for Eddington bias (Eddington, 1913) following $\ln M_{Edd}=\ln
M_{obs}+\frac{1}{2}\epsilon\sigma_{\ln M}^{2}$. $\ln M_{obs}$ is the stellar
mass estimate, $\sigma_{\ln M}$ is the uncertainty in the stellar mass taken
from the posterior PDF, and $\epsilon$ is the local slope of the halo mass
function. Jespersen et al. (2024) looks at the impact of cosmic variance on
the predictions of EVS, as this is not accounted for in the method of Lovell
et al. (2023) we use here. Given that our sample consists of multiple widely-
separated fields, we do not believe cosmic variance will play a significant
role on this EVS analysis.

No galaxies in our sample fall above the limit imposed by a stellar fraction
of unity, shown with a black line, which would require more stellar mass than
the available baryons to form stars. The majority of the EPOCHS v1 sample are
not in tension with the EVS limits and are shown with greyscale hexbins on
Figure 12, shaded according to number density within the bin. For galaxies
with redshift and stellar mass estimates which place them in possible tension
with the EVS limits assuming a realistic stellar fraction (more than
1.5$\times$ the mean value of the lognormal distribution, shown by the dashed-
dot line in Figure 12) we overlay the individual mass measurements as points
on Figure 12. For each of these galaxies we show the ‘fiducial’ Bagpipes
photo-$z$ and stellar mass, as well as the maximum and minimum stellar mass
estimate for that galaxy in our Bagpipes and Prospector fitting, including the
modified IMF results we discuss in § 4.5. We link individual mass estimates
for the same object with dotted lines. In the majority of cases, whilst the
maximum stellar mass estimate may suggest a possible tension, the minimum and
often the fiducial mass estimates are not in tension.

It is worth noting that all galaxies which show a possible tension at redshift
$z\leq 8$ in this figure are classed as ‘little red dots’ (LRDs), which are
discussed further in § 5.2.3 and potentially contain a significant AGN
component we do not account for. It is clear from Figure 12 that none of our
candidates are in tension with the predictions of $\Lambda$CDM, as seen by the
lack of sources above the solid black line, however a number of objects do
require very high stellar fractions at these redshifts. The galaxies which
require the highest star formation efficiencies seem to be found at $z\sim
7-8$, rather than the higher redshift probed in this study. This has also been
found for HST-dark galaxies observed as part of FRESCO, where the galaxies
with highest implied stellar mass densities are between $5\geq z\geq 6$ (Xiao
et al., 2023).

Further evidence for the compatibility of our results with standard
cosmological models can be seen in Figure 8, where we are below the SMF upper
limit calculated by Li et al. (2023) in all redshift bins, which is shown with
a black line. Li et al. (2023) present a bursty SFH model consisting of a
series of feedback-free bursts (on timescales of $\sim$10 Myr per burst) for
galaxies in halos above a given mass/redshift cutoff, resulting in higher star
formation efficiency, cosmic SFR density and stellar mass density above $z\geq
8$. The upper limit shown is computed assuming a maximum star formation
efficiency ($\epsilon_{\rm max}$) of unity. Our fiducial GSMF results at
$z\approx 7$ are close to the FFB predictions with $\epsilon_{\rm max}$ = 0.2,
but fall below this value at higher redshifts.

Wang et al. (2024) have shown that stellar masses of high mass galaxies are
typically overestimated by $\sim$0.4 dex when $>1\mu$m rest-frame emission is
not used in SED fitting, which requires MIRI observations in this redshift
regime. As MIRI observations are not used in this study, the degeneracies in
age-attenuation and the relative contributions of strong emission lines and
the stellar continuum observed by Wang et al. (2024) are not constrained,
possibly leading to an overestimation in stellar mass. Outshining and
stochastic star formation are however still likely to have an effect on
stellar masses derived from SED fitting, even when MIRI data is used
(Narayanan et al., 2023; Shen et al., 2023), and the derived discrepancy in
stellar mass without MIRI will likely depend on the assumed SFH, dust law and
parameter priors, as we have shown they can also systematically change stellar
mass estimates. Wang et al. (2024) also find that high mass galaxies at $6\leq
z\leq 8$ may require a higher star formation efficiency ($\epsilon\sim 0.3$)
than the local Universe, but they do not find any incompatibility with
standard cosmological models.

Figure 12: Photometric redshifts and stellar mass estimates for the EPOCHS v1
galaxy sample. We split the sample into two components, with the majority of
galaxies shown by the grayscale hexbins. For the most massive galaxies we show
individual datapoints in order to investigate possible tension with
$\Lambda$CDM. We use the Extreme-value Statistics (EVS) methodology of Lovell
et al. (2023) to calculate both the upper stellar mass limit assuming a
stellar fraction of unity (shown with a black line), as well as a more
realistic stellar fraction distribution parametrized by a truncated log-normal
distribution with $\mu=e^{-2}$ and $\sigma=1$ (median shown with a dash-dot
line, colored contours show 1-3$\sigma$ range). We assume a baryon fraction of
0.16 as predicted by Planck Collaboration et al. (2016) and use our full
survey area of 187 arcmin2 in the EVS calculation. For galaxies in the EPOCHS
v1 sample with a stellar mass estimate greater than 1.5$\times$ the limit for
a realistic stellar fraction (dash-dot line) we show the maximum, minimum and
fiducial stellar mass and redshift estimates for our Bagpipes and Prospector
fits, joined with a dotted black line. The shape and color of the points
references the SED-fitting tool, SFH and IMF used, as explained by the inset
legend. No galaxies exceed the maximum stellar mass predicted by $\Lambda$CDM,
although some galaxies do require a high stellar fraction, particularly the
’little red dots’ with high inferred masses at $z\sim 7$.

In Figure 3 we show examples of the photometry, best-fitting Bagpipes &
Prospector SEDs, and posterior redshift and stellar mass estimates for a few
of the most massive galaxies in the EPOCHS v1 sample at a range of redshifts.
For the galaxy labelled CEERSP1:7463, which is representative of the LRDs we
observe, the inferred star formation history for the non-parametric SED fits
suggest the bulk of stellar mass was formed $\sim$100-200 Myr ago, which may
put it and some of the other high-mass LRD galaxies in greater tension at
$\Lambda$CDM at earlier points in their star formation histories, if the
integrated SFH is considered. However we discuss in § 6.3, the masses and star
formation histories of the LRDs are highly uncertain, particularly when they
contribution of AGN is not considered.

###  6.2 Comparison of our sample to other studies

The study by Labbé et al. (2023) found a collection of unexpectedly massive,
high-redshift galaxies in the CEERS field. These objects typically have blue
colors in the short wavelength NIRCam filters, but extremely red colors at
longer wavelengths. We have cross-matched their 13 candidates to our catalogue
in order to compare our derived masses and redshifts. In the following text we
give both the ID in our catalogue followed by the ID given in Labbé et al.
(2023) for the sources we discuss here. 11 of their candidates appear within
our catalogue, of which 7 are selected by our selection criterion and are
among the most massive galaxies we find. The two sources we do not recover
appear to be blended with neighbouring sources in our SExtractor segmentation
maps. Of the four that we detect but do not select, 5346_CEERSP6 (11184) falls
within our image mask (close to the detector gap in the NIRCam SW
observations), and the other three fail our signal to noise requirement in the
rest-frame UV ($\geq$5 sigma). One of these, 2683_CEERSP3 (13050), has been
spectroscopically confirmed as a broadline AGN by Kocevski et al. (2023) at
$z=5.6$. We do not detect 37888 or 39575, which are among the lowest mass
candidates presented in Labbé et al. (2023). We do not consider 2 of the
selected galaxies in this study as their best-fitting redshift from EAZY-py is
marginally below our z$\geq$6.5 limit. For these galaxies, _8750_CEERSP3
(7274)_ and 2499_CEERSP1 (25666), we find lower redshift solutions than Labbé
et al. (2023), with $6.03_{-0.14}^{+0.36}$ and $6.45_{-0.18}^{+0.10}$ vs their
$7.77^{+0.05}_{-0.06}$ and $7.93^{+0.09}_{-0.16}$. For four of the remaining
five our photo-$z$ estimates agree well with Labbé et al. (2023), but for
1516_CEERSP2 (21834), we find a significantly higher photo-$z$ solution of
$10.4^{+1.3}_{-0.6}$ with EAZY-py. However our fiducial Bagpipes results,
which incorporate the EAZY-py photo-$z$ posterior as a prior, suggest a
photo-$z$ closer to $8.63^{+0.24}_{-0.32}$, which is in agreement with their
result of $8.54^{+0.32}_{-0.51}$. For these five galaxies, our stellar mass
estimates are on average a factor of 2.1$\times$ (0.32 dex) smaller than Labbé
et al. (2023), which is mostly due to the differing choice of IMF. Madau &
Dickinson (2014) suggests a scaling factor of $\approx 1.49\times$ to convert
a stellar mass estimate from Salpeter (1955) to Kroupa (2001), which does not
account for the full discrepancy between our mass estimates.

The most massive candidate presented by Labbé et al. (2023), 7463_CEERSP1
(38094), with a stellar mass of
$\log_{10}(M_{\star}/M_{\odot})=10.89^{+0.09}_{-0.08}$ is the second most
massive galaxy in our sample, with a fiducial Bagpipes mass of
$10.65^{+0.09}_{-0.10}$. Given that our sample contains $\>5\times$ the area
of the dataset used in Labbé et al. (2023), this may suggest that the CEERS
field is overdense. As a comparison we compute the GSMF implied by the results
of Labbé et al. (2023), assuming 100% completeness. The results of this are
shown in Figure 8. Their implied datapoint at both redshift 8 and 9 lies above
our best-fitting Schechter functions, due to a combination of a smaller cosmic
volume and higher stellar masses than our results. Wang et al. (2024) shows
that the stellar masses of high-z galaxies selected with the same color
criteria as Labbé et al. (2023) are reduced by $\sim 0.4$ dex when rest-frame
$>1\mu$m observations from MIRI are included in SED fitting, suggesting that
the stellar masses derived without MIRI data for the most massive galaxies in
this work and in Labbé et al. (2023) may be systematically overestimated.

Although our source detection is done at the longest available wavelengths,
our JWST sample is primarily selected in the rest-frame UV due to our Lyman-
break criteria. Galaxies with weak UV emission may not be selected or even
detected in our sample, whilst still containing significant stellar mass.
Highly obscured galaxies like Submillimeter Galaxies (SMGs), which have high
stellar masses, will not be included in our selection criteria, although they
are thought to have a low spatial density (e.g. Smail et al., 2023).

Whilst high-$z$ galaxies are predicted to be relatively compact (Roper et al.,
2022; Shen et al., 2024), the high resolution of JWST and relatively small
extraction apertures used also means we are more sensitive to compact, high-
surface brightness sources than extended, low-surface brightness sources. If a
population of extended, diffuse galaxies exists at these redshifts they are
likely to be excluded from this sample.

###  6.3 Contamination of the GSMF by hidden AGN

As discussed in § 5.2.3, obscured high-$z$ AGN can impact the high-mass end of
the GSMF at $z=7-8$ by masquerading as galaxies with evolved stellar
populations and therefore higher masses. Kokorev et al. (2024) presents
candidate LRDs in both the CEERS and GOODS-South fields, which covers the
NGDEEP and JADES-Deep observations. We compare their catalogue to our high-$z$
galaxy sample, finding 5(3) of their CEERS (GOODS-South) candidates in our
robust sample. They are all at $z\leq 8$ in our sample, with IDs 1080_CEERSP9,
2184_CEERSP1, 4842_CEERSP6, 7463_CEERSP1, 7520_CEERSP5, 6169_JADES-Deep-GS,
31738_JADES-Deep-GS and 36222_JADES-Deep-GS. It is worth noting that some of
these galaxies in CEERS were first found in Labbé et al. (2023) and are
discussed in § 6.2. We detect almost all of the $z\geq 6.5$ galaxies from
Kokorev et al. (2024) in the CEERS and GOODS-South fields (18 and 10
respectively), but for the galaxies that do not appear in our final sample we
either find a lower redshift solution (below our $z\geq 6.5$ cutoff), or we do
not robustly detect the Lyman break to 5$\sigma$, meaning they do not meet our
selection criteria. We see good agreement in redshift, with a maximum offset
of $\Delta z\approx 0.4$, and a mean offset of $\Delta z\approx 0.14$. These