# EPOCHS III: Unbiased UV continuum slopes at $6.5<z<13$ from combined PEARLS
GTO and public JWST NIRCam imaging

Duncan Austin Jodrell Bank Centre for Astrophysics, Alan Turing Building,
University of Manchester, Oxford Road, Manchester M13 9PL, UK  Christopher J.
Conselice Jodrell Bank Centre for Astrophysics, Alan Turing Building,
University of Manchester, Oxford Road, Manchester M13 9PL, UK  Nathan J. Adams
Jodrell Bank Centre for Astrophysics, Alan Turing Building, University of
Manchester, Oxford Road, Manchester M13 9PL, UK  Thomas Harvey Jodrell Bank
Centre for Astrophysics, Alan Turing Building, University of Manchester,
Oxford Road, Manchester M13 9PL, UK  Qiao Duan Jodrell Bank Centre for
Astrophysics, Alan Turing Building, University of Manchester, Oxford Road,
Manchester M13 9PL, UK  James Trussler Jodrell Bank Centre for Astrophysics,
Alan Turing Building, University of Manchester, Oxford Road, Manchester M13
9PL, UK  Qiong Li Jodrell Bank Centre for Astrophysics, Alan Turing Building,
University of Manchester, Oxford Road, Manchester M13 9PL, UK  Ignas
Juodžbalis Jodrell Bank Centre for Astrophysics, Alan Turing Building,
University of Manchester, Oxford Road, Manchester M13 9PL, UK  Kavli Institute
for Cosmology, University of Cambridge, Cambridge, Madingley Road, Cambridge,
CB3 0HA  Katherine Ormerod Jodrell Bank Centre for Astrophysics, Alan Turing
Building, 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  Lewi Westcott Jodrell Bank Centre for Astrophysics, Alan
Turing Building, University of Manchester, Oxford Road, Manchester M13 9PL, UK
Honor Harris Jodrell Bank Centre for Astrophysics, Alan Turing Building,
University of Manchester, Oxford Road, Manchester M13 9PL, UK  Stephen M.
Wilkins Astronomy Centre, Department of Physics and Astronomy, University of
Sussex, Brighton, BN1 9QH, 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  Joseph Caruana Department of
Physics, University of Malta, Msida MSD 2080, Malta  Institute of Space
Sciences & Astronomy, University of Malta, Msida MSD 2080, Malta  Dan Coe
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218,
USA  Association of Universities for Research in Astronomy (AURA) for the
European Space Agency (ESA), STScI, Baltimore, MD 21218, USA  Center for
Astrophysical Sciences, Department of Physics and Astronomy, The Johns Hopkins
University, 3400 N Charles St. Baltimore, MD 21218, USA  Seth H. Cohen School
of Earth and Space Exploration, Arizona State University, Tempe, AZ
85287-1404, USA  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
Jordan C. J. D’Silva 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  ARC Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia
Brenda Frye Department of Astronomy/Steward Observatory, University of
Arizona, 933 N Cherry Ave, Tucson, AZ, 85721-0009, USA  Lukas J. Furtak
Physics Department, Ben-Gurion University of the Negev, P.O. Box 653, Be’er-
Sheva 84105, 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,
Natural Science Building 102, Louisville KY 40292, USA  Rolf A. Jansen School
of Earth and Space Exploration, Arizona State University, Tempe, AZ
85287-1404, USA  Anton M. Koekemoer Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD 21218, USA  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, 34124 Trieste, Italy
Rafael Ortiz III School of Earth and Space Exploration, Arizona State
University, Tempe, AZ 85287-1404, USA  Nor Pirzkal Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Aaron Robotham
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  Russell E. Ryan, Jr. Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Jake Summers School of Earth and Space Exploration, Arizona State University,
Tempe, AZ 85287-1404, USA  Christopher N. A. Willmer Steward Observatory,
University of Arizona, 933 N Cherry Ave, Tucson, AZ, 85721-0009, USA  Rogier
A. Windhorst School of Earth and Space Exploration, Arizona State University,
Tempe, AZ 85287-1404, USA  Haojing Yan Department of Physics and Astronomy,
University of Missouri, Columbia, MO 65211, USA  Erik Zackrisson Observational
Astrophysics, Department of Physics and Astronomy, Uppsala University, Box
516, SE-751 20 Uppsala, Sweden  Swedish Collegium for Advanced Study,
Linneanum, Thunbergsvägen 2, SE-752 38 Uppsala, Sweden

###### Abstract

We present an analysis of rest-frame UV continuum slopes, $\beta$, using a
sample of 1011 galaxies at $6.5<z<13$ from the EPOCHS photometric sample
collated from the GTO PEARLS and public ERS/GTO/GO (JADES, CEERS, NGDEEP,
GLASS) _JWST_ NIRCam imaging across $178.9~{}\mathrm{arcmin}^{2}$ of unmasked
blank sky. We correct our UV slopes for the photometric error coupling bias
using $200,000$ power law SEDs for each $\beta=\\{-1,-1.5,-2,-2.5,-3\\}$ in
each field, finding biases as large as $\Delta\beta\simeq-0.55$ for the lowest
SNR galaxies in our sample. Additionally, we simulate the impact of rest-UV
line emission (including Ly$\alpha$) and damped Ly$\alpha$ systems on our
measured $\beta$, finding biases as large as $0.5-0.6$ for the most extreme
systems. We find a decreasing trend with redshift of $\beta=-1.51\pm
0.08-(0.097\pm 0.010)\times z$, with potential evidence for Pop. III stars or
top-heavy initial mass functions (IMFs) in a subsample of 68
$\beta+\sigma_{\beta}<-2.8$ galaxies. At $z\simeq 11.5$, we measure an
extremely blue $\beta(M_{\mathrm{UV}}=-19)=-2.73\pm 0.06$, deviating from
simulations, indicative of low-metallicity galaxies with non-zero Lyman
continuum escape fractions $f_{\mathrm{esc,LyC}}\gtrsim 0$ and minimal dust
content. The observed steepening of
$\mathrm{d}\beta/\mathrm{d}\log_{10}(M_{\star}/\mathrm{M}_{\odot})$ from
$0.22\pm 0.02$ at $z\simeq 7$ to $0.81\pm 0.13$ at $z\simeq 11.5$ implies that
dust produced in core-collapse supernovae (SNe) at early times may be ejected
via outflows from low mass galaxies. We also observe a flatter
$\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}=0.03\pm 0.02$ at $z\simeq 7$ and a
shallower $\mathrm{d}\beta/\mathrm{d}\log_{10}(M_{\star}/\mathrm{M}_{\odot})$
at $z<11$ than seen by _HST_ , unveiling a new population of low mass, faint,
galaxies reddened by dust produced in the stellar winds of asymptotic giant
branch (AGB) stars or carbon-rich Wolf-Rayet binaries.

##  1 Introduction

The first year of study with the James Webb Space Telescope (_JWST_ , Gardner
et al., 2023) has led to the discovery of a multitude of high-redshift galaxy
candidates in the first billion years of cosmic evolution (Castellano et al.,
2022; Finkelstein et al., 2022; Harikane et al., 2022; Naidu et al., 2022; Yan
et al., 2022; Adams et al., 2023a; Atek et al., 2023; Austin et al., 2023;
Bouwens et al., 2023a; Donnan et al., 2023; Leung et al., 2023) discovered in
deep near-infrared (NIR) imaging using the Near-InfraRed Camera (NIRCam, Rieke
et al., 2005, 2023) instrument. In addition, several studies including Tang et
al. (2023); Fujimoto et al. (2023); Nakajima et al. (2023); and Bunker et al.
(2023) have spectroscopically confirmed galaxies using JWST’s Near InfraRed
Spectrograph (NIRSpec) Micro Shutter Assembly (MSA, Ferruit et al., 2022;
Jakobsen et al., 2022; Rawle et al., 2022; Böker et al., 2023), with the most
distant sources found at $z\sim 12-13$ (e.g. Curtis-Lake et al., 2023; Wang et
al., 2023; Castellano et al., 2024).

In this early JWST data, measurements of the ultraviolet luminosity function
(UV LF) have revealed an overabundance in intrinsically UV bright sources at
$z>10$ (Adams et al., 2023b; Donnan et al., 2023; Finkelstein et al., 2023a;
Bouwens et al., 2023b; Pérez-González et al., 2023; Castellano et al., 2023;
McLeod et al., 2023; Leung et al., 2023; Finkelstein et al., 2023b; Donnan et
al., 2024), which may be explained by either the mapping of UV luminosity to
host halo mass (e.g. Mason et al., 2023), an increased star formation
efficiency (SFE, e.g. Inayoshi et al., 2022; Mirocha & Furlanetto, 2023), or
minimal dust obscuration at high redshift (Ferrara et al., 2023; Ziparo et
al., 2023). Further investigation into the rest-frame UV properties of these
galaxies is therefore of the utmost importance to help distinguish between
these plausible scenarios.

The rest-frame UV continuum slope, $\beta$, is commonly calculated via the
power law $f_{\lambda}\propto\lambda^{\beta}$ (Calzetti et al., 1994)
(henceforth C94), and is a key diagnostic for understanding the properties of
the stellar continuum, providing a tracer of massive O/B-type main sequence
stars and their surrounding HII nebular regions. It is primarily dependent on
the galactic dust content, meaning the UV dust attenuation ($A_{\mathrm{UV}}$)
can be estimated from $\beta$ via the empirical Meurer et al. (1999) relation,
$A_{\mathrm{UV}}=4.43+1.99\beta$. The $100~{}\mathrm{Myr}$ star formation
rates (SFRs) of star-forming galaxies (SFGs) can therefore be estimated from
dust-corrected intrinsic UV magnitudes, $M_{\mathrm{UV}}$, under an assumed
luminosity to SFR conversion ($\kappa_{\mathrm{UV}}$, see Madau & Dickinson,
2014). Additionally, $\beta$ has a more minor dependence on stellar ages and
metallicities (see e.g. Bouwens et al., 2012; Tacchella et al., 2022),
especially in bursty star formation history (SFH) models, evidence for which
has been observed in a spectroscopically confirmed $z=7.3$ galaxy by Looser et
al. (2023).

Since $\beta$ is closely related to the dust content in galaxies, we can use
it to study the global build up of galactic dust. It is expected that the
dominant dust production mechanism at early times is nucleation in the ejecta
of core-collapse supernovae (SNe, e.g. Todini & Ferrara, 2001; Bianchi et al.,
2009; Marassi et al., 2019), although dust destruction processes such as
sputtering, sublimation, and grain-grain collisions in the resulting reverse
shock are not well constrained even in the local Universe (e.g. Bocchio et
al., 2016; Micelotta et al., 2018; Kirchschlager et al., 2019, 2024). As well
as core-collapse supernovae (SNe), dust production also occurs in the stellar
winds of $0.8-8~{}\mathrm{M}_{\odot}$ asymptotic giant branch (AGB) stars
(e.g. Zhukovska et al., 2008; Gail et al., 2009; Höfner & Olofsson, 2018),
$12-30~{}\mathrm{M}_{\odot}$ red super giant (RSG) stars (Levesque et al.,
2006), and $\gtrsim 30~{}\mathrm{M}_{\odot}$ Wolf-Rayet (WR) stars that are
carbon-rich (WC stars) with an OB companion (e.g. Lau et al., 2020, 2022),
although the accompanying SN is expected to very quickly destroy dust produced
in both RSGs and WCs. Dust reprocessing is also expected in the interstellar
medium (ISM) via astration and depletion (e.g. Draine & Salpeter, 1979;
Draine, 2009, 2011), although these processes are not yet fully understood
(see Schneider & Maiolino (2023) for a thorough review). The dust chemistry
and grain size distribution impact the extinction curve shape (Salim &
Narayanan, 2020), which is often assumed to follow the Calzetti et al. (2000),
LMC (Gordon et al., 2003), or SMC (Pei, 1992) laws. Recent results regarding
the relation between IR excess
($\mathrm{IRX}=\log_{10}(L_{\mathrm{IR}}/L_{\mathrm{UV}})$) and $\beta$ (the
IRX-$\beta$ relation) from the Atacama Large Millimeter/submillimeter Array
(ALMA) Reionization Era Bright Emission Line Survey (REBELS, Bouwens et al.,
2022) data have, however, demonstrated the suitability of a “Calzetti-like”
attenuation curve at $z\simeq 7$ (Bowler et al., 2024).

Prior to the onset of JWST, much effort had been made to measure $\beta$ in
large, unbiased, photometric galaxy samples at $z\sim 4-10$ collated from deep
imaging by Hubble Space Telescope’s (HST) Advanced Camera for Surveys (ACS),
Near Infrared Camera and Multi-Object Spectrometer (NICMOS, e.g. Hathi et al.,
2008), and Wide Field Camera 3 (WFC3) Infra-Red (IR) instruments. The most
notable work at $z=7-10$ has been conducted with WFC3-IR, where $\beta$ is
measured either directly from F105W (Y-band), F125W (J-band), F140W (JH-band)
and F160W (H-band) colors (Bouwens et al., 2009, 2010; Finkelstein et al.,
2010; McLure et al., 2011; Dunlop et al., 2012, 2013; Bouwens et al., 2014;
Wilkins et al., 2016; Bhatawdekar & Conselice, 2021), or from the best-fitting
SED template (e.g. Finkelstein et al., 2012) in the 10 C94 filters designed to
omit UV absorption/emission features present in galaxy spectra. These studies
find relatively red $\beta\sim-1.9$ at $z\simeq 4$ (e.g. Hathi et al., 2013)
and moderately blue $\beta\sim-2.2/-2.4$ at $z\simeq 7$ in galaxies with
$M_{\mathrm{UV}}\simeq-19$, symbolic of low-metallicity, relatively dust-free
systems, corroborating predictions from simulations (Wilkins et al., 2012,
2013; Gonzalez-Perez et al., 2013).

More recently, the redder JWST/NIRCam filters have given access to the rest-
frame UV at $z\gtrsim 10$, allowing for both the first $\beta$ measurements at
these redshifts and improved $\beta$ measurements from vastly deeper rest-
frame NUV/optical coverage for galaxies at $z\gtrsim 5$. Photometric $\beta$
measurements have been made in abundance with JWST (Topping et al., 2022;
Cullen et al., 2023a; Nanayakkara et al., 2023; Morales et al., 2023), with
Topping et al. (2023) and Cullen et al. (2023b) observing, on average,
extremely blue UV slopes at $z\gtrsim 10$, implying little-to-no dust presence
at these epochs. In addition, JWST studies have found tentative evidence for
$\beta<-3$ in individual galaxies at $z\gtrsim 10$ (e.g. Atek et al., 2023;
Austin et al., 2023; Cullen et al., 2023a; Furtak et al., 2023), although
these may be a result of either photometric scatter or the known $\beta$ bias
from the increased selectability of faint blue galaxies compared to their
redder counterparts (see e.g. Rogers et al., 2013). These ultra-blue sources
may, however, provide evidence for exotic Pop. III stellar populations
(Zackrisson et al., 2011) and/or top-heavy/bottom-light IMFs (e.g. Rasmussen
Cueto et al., 2023; Steinhardt et al., 2023). As well as this, these blue
galaxies require high Lyman Continuum (LyC) escape fractions,
$f_{\mathrm{esc,LyC}}$, needed to adequately reduce the reddening of the
continuum by free-free, free-bound and two-photon nebular emission (Schaerer,
2002, 2003; Raiter et al., 2010a; Byler et al., 2017; Topping et al., 2022).
These $f_{\mathrm{esc,LyC}}$ are unfortunately impossible to measure directly
due to the opacity of the intergalactic medium (IGM) in the Epoch of
Reionization (EoR). Several indirect tracers have been suggested including
blue UV $\beta$ slopes (Zackrisson et al., 2013; Chisholm et al., 2022), low
[OIII]+H$\beta$ equivalent widths (EWs, Topping et al., 2022; Endsley et al.,
2023), small sizes (Mascia et al., 2023), strong Mg II-$\lambda\lambda
2796,2803$ Chisholm et al. (2020); Katz et al. (2022), etc, although it is
likely that a combination of these are required to accurately measure
$f_{\mathrm{esc,LyC}}$ (e.g. Choustikov et al., 2023). This is necessary to
estimate the relative importance of low mass galaxies, which are expected to
contribute a significant fraction of the total ionizing photon budget of
galaxies needed to reionize the neutral IGM by $z\simeq 6$ (see e.g. Bouwens
et al., 2015a; Planck Collaboration et al., 2016, 2020).

In this work we use the EPOCHS v1 sample photometric galaxy sample, presented
in EPOCHS-I (Conselice et al. in prep.), comprised from deep public and PEARLS
GTO (PI: R. Windhorst & H. Hammel, PIDs: 1176 & 2738, Windhorst et al., 2023)
JWST/NIRCam imaging. The content of this paper is as follows. We outline our
data sources and cataloguing procedure in § 2 before calculating UV properties
in § 3. $\beta$ biases are discussed in § 4 and $\beta$ scaling relations are
presented in § 5. We discuss the dust implications and possibility of exotic
Pop. III and top-heavy IMF scenarios in § 6 before concluding in § 7.

$\Lambda$CDM cosmological parameters of $\Omega_{m}=0.3$,
$\Omega_{\Lambda}=0.7$, $H_{0}=70\text{ kms}^{-1}\text{Mpc}^{-1}$, and AB
magnitudes (Oke, 1974; Oke & Gunn, 1983) are assumed throughout this analysis.

##  2 Data and cataloguing

###  2.1 PEARLS imaging

In this work, we utilize blank field JWST/NIRCam photometric imaging from the
Prime Extragalactic Areas for Reionization Science (PEARLS; PI: R. Windhorst &
H. Hammel, PIDs: 1176 & 2738, Windhorst et al., 2023) GTO survey. This
includes 3 NIRCam parallel fields surrounding the Hubble Frontier Field (HFF,
Lotz et al., 2017) MACS J0416.1-2403 (hereafter MACS-0416) lensing cluster at
$z=0.396$ and the El Gordo ($z=0.87$, Menanteau et al., 2012) cluster
parallel. In addition, we also include data from the 4 spokes of the North
Ecliptic Pole Time Domain Field (NEP-TDF, Jansen & Windhorst, 2018). As part
of the PEARLS survey design, these fields contain imaging in 4 short
wavelength (SW; F090W, F115W, F150W, F200W) and 4 long-wavelength (LW; F277W,
F356W, F410M, F444W) NIRCam filters. In addition, we also include bluer F606W
imaging from the HST ACS Wide Field Channel (WFC) in the NEP-TDF from the
GO-15278 (PI: R. Jansen) and GO-16252/16793 (PIs: R. Jansen & N. Grogin, see
O’Brien et al., 2024) HST programs. At the time of writing, the first 2 NEP-
TDF spokes, the first epoch of MACS-0416, and the El Gordo cluster are
publicly available from the Minkulski Archive for Space Telescopes (MAST)111An
overview of PEARLS data and papers, as well as access to both independently
reduced NIRCam imaging and cluster lens models are available at
https://sites.google.com/view/jwstpearls..

###  2.2 Public ERS and GO imaging

As part of the EPOCHS sample, we include public ERS and GO imaging from the
Grism Lens-Amplified Survey from Space (GLASS; PI: T. Treu, PID: 1324, Treu et
al., 2022), data release 1 of deep JWST Advanced Deep Extragalactic Survey
imaging in GOODS-South (JADES-Deep-GS; PI: D. Eisenstein, PID: 1180,
Eisenstein et al., 2023; Rieke et al., 2023), Cosmic Evolution Early Release
Science (CEERS; PI: S. Finkelstein, PID: 1345, Bagley et al., 2023a) and Next
Generation Deep Extragalactic Exploratory Public (NGDEEP; PIs: S. Finkelstein,
C. Papovich & N. Pirzkal, PID: 2079, Bagley et al., 2023b) surveys.
Photometric NIRCam imaging is available in F115W, F150W, F200W, F277W, F356W,
and F444W for every survey. Where available we include the F090W (GLASS and
JADES) wideband filter and the F410M (CEERS and JADES) and F335M (JADES)
medium band NIRCam filters.

In addition, we also include bluer data from HST ACS/WFC for CEERS and NGDEEP
to somewhat compensate for the lack of blue F090W NIRCam photometry. We use
ACS/WFC data in the F606W and F814W filters from the Cosmic Assembly Near-IR
Deep Extragalactic Legacy Survey (CANDELS; PIs: S. Faber & H. Ferguson, DOI:
10.17909/T94S3X, Grogin et al., 2011; Koekemoer et al., 2011) covering the
Extended Groth Strip (EGS, Groth et al., 1994) for CEERS, and deep F435W,
F606W, and F814W ACS/WFC data from v2.5 of the Hubble Legacy Fields (HLF,
Whitaker et al., 2019) covering the HUDF-Par2 for NGDEEP.

###  2.3 NIRCam data reduction pipeline

Our data reduction pipeline is similar to that of Ferreira et al. (2022);
Adams et al. (2023a); and Austin et al. (2023), and is identical to that used
in the rest of the EPOCHS series. A detailed account of this pipeline is
presented in Adams et al. (2023b) (henceforth EPOCHS-II), which is summarized
below.

We use version 1.8.2 of the official JWST data reduction pipeline (Bushouse et
al., 2022)222https://github.com/spacetelescope/jwst and Calibration Reference
Data System (CRDS) v1084, which includes improved in-flight LW flat fielding
which dramatically deepens these images compared to CRDS v0995. Between stage
1 and 2 we subtract templates of ‘wisps’ in F150W and F200W (SW artefacts
produced by reflection off the top secondary mirror support strut, which are
most prominent in the B4 detector) using the official STScI
templates.333https://stsci.app.box.com/s/1bymvf1lkrqbdn9rnkluzqk30e8o2bne We
do not remove ‘claws’ (artefacts from nearby bright foreground stars within
the susceptibility region) in the NIRCam imaging, but instead mask these post
reduction. After stage 2, we apply Chris Willott’s $1/f$ noise
correction.444See https://github.com/chriswillott/jwst Before stage 3 we
perform background subtraction on individual ‘cal.fits’ frames, which consists
of an initial flat background subtraction followed by a further 2D background
subtraction using photutils (Bradley et al., 2022). This skips the sky
subtraction step in stage 3 and allows for quicker background subtraction
assessment and fine-tuning. Post stage 3, we use the tweakreg part of the
DrizzlePac555https://github.com/spacetelescope/drizzlepac python package to
first align the F444W image onto the Gaia-DR3 (Data Release 3, Gaia
Collaboration et al., 2018; Vallenari et al., 2022) derived World Co-ordinate
System (WCS) before matching all remaining filters to this WCS. We finally
pixel-match to the F444W image with the use of astropy reproject (Hoffmann et
al., 2021)666https://reproject.readthedocs.io/en/stable/ to make drizzled
images with a pixel scale of $0\farcs 03$/pixel.

###  2.4 Catalogue creation and photo-$z$’s

To detect sources in each field, we perform forced photometry in apertures
detected in the inverse-variance weighted stacked F277W, F356W, and F444W
reduced images using SExtractor (Bertin & Arnouts, 1996) with the setup given
in EPOCHS-II to produce an initial photometric catalogue. This includes band
fluxes and magnitudes in both Kron (“FLUX_AUTO” and “MAG_AUTO”, Kron, 1980)
and 0$\farcs$32 diameter circular (“FLUX_APER” and “MAG_APER”) apertures as
well as an estimate of the half light radius (“FLUX_RADIUS”), among other
quantities. Once photometric fluxes have been measured for our JWST/NIRCam
images, we correct for the missing flux that is spread beyond our aperture by
the point spread function (PSF) using the simulated WebbPSF (Perrin et al.,
2012, 2014) PSFs. Our $0\farcs 32$ apertures contain approximately $70-80\%$
of the total point source flux to get the best balance of depth, number of
pixels used and minimize the aperture correction obtained from the potentially
uncertain WebbPSF models.

For the Hubble images with differing dimensions to our reduced JWST imaging,
we use the photutils (Bradley et al., 2022) python module to perform forced
photometry in the same $0\farcs 32$ diameter circular apertures from the
F277W, F356W, and F444W JWST/NIRCam stack. We also confirm that the fluxes we
obtain are consistent with those measured using sep (SExtractor for python;
Barbary, 2016). We PSF correct our Hubble ACS/WFC photometry using the PSF
models given in Bohlin (2016).

To calculate image depths in each band, we randomly place $0\farcs 32$
diameter apertures in empty regions defined by both the SExtractor
segmentation map and our band dependent image masks. These image masks are
manually created to mask out stellar diffraction spikes, a $50-100$ pixel
border near the shallower image edges, as well as any remaining image
artefacts, such as stray ‘snowballs’, ‘wisps’, or ‘claws’. Each aperture used
in these depth calculation is forced to not overlap any other aperture as well
as being at least $1\arcsec$ away from any source. The flux measurements of
the nearest 200 apertures to each detected source are used to estimate a
“local depth” and hence provide a photometric error defined as the Normalized
Median Absolute Deviation (NMAD) of measured fluxes. We choose to adopt this
flux error as opposed to the underestimated SExtractor output as this more
appropriately deals with the correlated noise in the photometric images. To
account for potential JWST/NIRCam zero point (ZP) uncertainties, we also
impose a $10\%$ minimum flux error in each band. Table 1 shows the average
$5\sigma$ depths obtained in each NIRCam band for the JWST surveys used in
this work. These are calculated as the NMAD of all apertures in empty image
regions as opposed to the average local depth to avoid aperture double
counting.

We next calculate photometric redshifts (photo-$z$’s) for our EPOCHS sample
using the EAZY-py (Brammer et al., 2008) spectral energy distribution (SED)
fitting code using the standard “tweak_fsps” templates produced using the
flexible stellar population synthesis (FSPS) package (Conroy & Gunn, 2010)
combined with set 1 and 4 from the bluer Larson et al. (2023) template set
(referred to henceforth as “fsps_larson”) with stronger emission line
treatment. The 3 pure stellar templates from Larson et al. (2023) (set 1) use
v2.2.1 of the Binary Populations and Spectral Synthesis (BPASS; Eldridge et
al., 2017; Stanway & Eldridge, 2018; Byrne et al., 2022) stellar population
synthesis (SPS) model for ages $\\{10^{6},10^{6.5},10^{7}\\}$ Gyr at a fixed
metallicity $Z_{\star}=0.05~{}\mathrm{Z}_{\odot}$ while adopting a standard
Chabrier (2003) IMF. In addition, template set 4 includes nebular continuum
and line emission calculated from CLOUDY v17 (Ferland et al., 2017) with
Ly$\alpha$ removed, ionization parameter $\log U=-2$,
$Z_{\mathrm{gas}}=Z_{\star}$, Lyman continuum escape fraction
$f_{\mathrm{esc,LyC}}=0$, $n_{\mathrm{H}}=300~{}\mathrm{cm}^{-2}$ gas cloud
hydrogen density, and spherical geometry. We run the EAZY-py photo-z fitting
twice with redshift both free ($0<z<25$) as well as a ‘low-redshift’ run,
where the redshift is limited to be $z<6$, in order to determine the
$\chi^{2}$ of the best-fitting low-z solution. This is required for our
selection criteria outlined in § 2.5 and has previously been performed for
high-$z$ galaxy identification in the past (e.g. Hainline et al., 2023a).

###  2.5 High-$z$ sample selection

One potential source of contamination in our EPOCHS sample, other than low-z
galaxy interlopers, comes from Milky Way Y- and T-type brown dwarfs. To
account for these, we fit the Sonora Bobcat brown dwarf templates (Marley et
al., 2021), to our candidate galaxies via least squares regression find the
best-fitting solution. In addition, we note that our sample may well contain
“little red dots” (LRDs), as found by (e.g. Matthee et al., 2023; Kokorev et
al., 2024) which are hypothesized to be active galactic nuclei (AGN; see Labbe
et al., 2023). The impact of contamination by these LRDs on the results of
this work is presented in Appendix A and an analysis of the impact on the
global stellar mass function (GSMF) is given in EPOCHS-IV (Harvey et al.,
2024).

To select a robust sample of high redshift galaxies from our photometric
catalogues, we use similar selection criteria as used in EPOCHS-I/II/IV,
outlined below:

  1. 1.

There must be at least 1 photometric band entirely bluewards of the Ly$\alpha$
break at $\lambda_{\mathrm{rest}}=1216~{}\mathrm{\AA}$. The limits used here
are defined by the upper and lower 50% transmission limits of the filters
taken from the Spanish Virtual Observatory (SVO; Rodrigo & Solano, 2020)
filter profile service. This constrains us to $z\gtrsim 6.5$ when using
JWST/NIRCam data alone.

  2. 2.

Any bands entirely bluewards of the Ly$\alpha$ break must be non-detected at
the $3\sigma$ confidence level, with bands straddling the break having no SNR
requirements.

  3. 3.

The two bands directly redwards (but not straddling) the Ly$\alpha$ break must
be detected at $>5\sigma$, and all other redwards NIRCam widebands must be
detected at $>2\sigma$.

  4. 4.

The integral of the output EAZY-py probability density function (PDF) must
have greater than $60\%$ contained within 10% of the maximum likelihood of the
EAZY-py redshift posterior.

  5. 5.

The best-fit EAZY-py (redshift free) solution must satisfy
$\chi^{2}_{\mathrm{red.}}<3$, and the difference between the ‘redshift free’
and ‘low redshift’ best-fitting solutions must satisfy $\Delta\chi^{2}>4$.

  6. 6.

The SExtractor half-light radius must be greater than or equal to $1.5$ pixels
in the LW F277W, F356W, and F444W JWST/NIRCam images to ensure the removal of
F200W drop-out (dithered) hot-pixels.

  7. 7.

If the SExtractor half-light radius is smaller than the PSF full width at half
maximum (FWHM) in F444W, then the difference between the best-fitting brown
dwarf solution and the ‘redshift free’ EAZY-py run must satisfy
$\Delta\chi^{2}_{\mathrm{red.}}>4$.

In addition to the above selection criteria, we also remove rogue candidates
from our sample by eye. This check was done independently by authors DA, TH,
QL, and NA. This splits our sample into “certain”, “uncertain” and “rejected”
candidates. In total we have 1214 unmasked galaxies and 59 brown dwarf
candidates in the full EPOCHS catalogue, which also includes the blank NIRCam
parallel fields of the SMACS-0723 and CLIO777The CLIO cluster is identified at
$z=0.42$ within the Galaxy and Mass Assembly (GAMA, Driver et al., 2011; Liske
et al., 2015) survey and imaged with NIRCam as part of the GTO PEARLS JWST
program. strong gravitational lensing clusters. Our “rejected” subsample
contains 49 sources removed completely by eye ($4\%$ of the full sample) which
are mainly sources containing obvious LW hot pixels which just pass our
criterion 6 above. This leaves 111 “uncertain” galaxy candidates ($9\%$ of the
full sample) and 1054 “certain” galaxy candidates. These “uncertain”
candidates mainly comprise sources that are contaminated by stray wisps,
claws, or diffuse foreground light. When limiting this to the fields and
redshift range used in this work ($6.5<z<13$), and removing all galaxies at
$6.5<z<7.5$ in fields without bluer HST ACS/WFC data (El Gordo, MACS-0416, and
GLASS) due to the extreme blue $\beta$ biases present (see § 3.1), we are left
with 1011 galaxies in our EPOCHS-III sample covering
$178.9~{}\mathrm{arcmin}^{2}$ of unmasked blank sky area (excluding CLIO and
SMACS-0723 included in the full EPOCHS sample). Full catalogues and an
associated README file will be distributed as part of EPOCHS-I.

Table 1: $5\sigma$ depths (in $0\farcs 32$ diameter apertures) and areas of the blank field surveys used in the EPOCHS-III sample. These differ from those presented in EPOCHS-I as we limit our redshift range to $z<13$ and also do not use the GTO PEARLS CLIO and ERS SMACS-0723 parallel fields in this analysis. This sample also differs slightly to EPOCHS-II/EPOCHS-IV due to our exclusion of galaxies at $6.5<z<7.5$ in El Gordo, MACS-0416, and GLASS due to the large $\beta$ biases associated with this redshift range. We note that NGDEEP is split into three regions; one has deep HST ACS/WFC coverage (4.03 arcmin2) and another has shallower ACS/WFC coverage (1.28 arcmin2), leaving 0.98 arcmin2 with JWST NIRCam data alone. These HST areas and depths for NGDEEP are collated from the work of Austin et al. (2023). Survey | HST ACS/WFC | JWST NIRCam | Area /
---|---|---|---
F606W | F814W | F090W | F115W | F150W | F200W | F277W | F335M | F356W | F410M | F444W | arcmin2
NEP-TDF | 28.74 | - | 28.50 | 28.50 | 28.50 | 28.65 | 29.15 | - | 29.30 | 28.55 | 28.95 | 57.32
MACS-0416 | - | - | 28.67 | 28.62 | 28.49 | 28.64 | 29.16 | - | 29.33 | 28.74 | 29.07 | 12.30
El Gordo | - | - | 28.23 | 28.25 | 28.18 | 28.43 | 28.96 | - | 29.02 | 28.45 | 28.83 | 3.90
CEERS P1-8,10 | 28.60 | 28.30 | - | 28.70 | 28.60 | 28.89 | 29.20 | - | 29.30 | 28.50 | 28.85 | 60.31
CEERS P9 | 28.31 | 28.32 | - | 29.02 | 28.55 | 28.78 | 29.20 | - | 29.22 | 28.50 | 29.12 | 6.08
NGDEEP | 29.20/30.30 | 28.80/30.95 | - | 29.78 | 29.52 | 29.48 | 30.28 | - | 30.22 | - | 30.22 | 6.29
JADES-Deep-GS | 29.07 | - | 29.58 | 29.78 | 29.68 | 29.72 | 30.21 | 29.58 | 30.17 | 29.64 | 29.99 | 22.98
GLASS | - | - | 29.14 | 29.11 | 28.86 | 29.03 | 29.55 | - | 29.61 | - | 29.84 | 9.76
|  |  |  |  |  |  |  |  |  |  |  |

###  2.6 Completeness and contamination

To estimate the completeness and contamination in our sample, we utilize 5
realizations of the JAdes extraGalactic Ultradeep Artificial Realizations
(JAGUAR; Williams et al., 2018) mock catalogue of star-forming galaxies
(SFGs). The mock photometry in the catalogue is generated from a family of
BEAGLE (Chevallard & Charlot, 2016) SEDs incorporating the Bruzual & Charlot
(2003) SPS and Gutkin et al. (2016) nebular emission models. The distribution
of redshifts and masses of these SEDs are selected to match the mass functions
of Tomczak et al. (2014), extrapolated to match the $z>4$ HST UV luminosity
functions (UV LFs) of Bouwens et al. (2015b, 2016); Calvi et al. (2016);
Stefanon et al. (2017); Oesch et al. (2013, 2018), with $M_{\mathrm{UV}}$
calculated from the $M_{\mathrm{UV}}-M_{\star}$ relations from 3D-HST (Skelton
et al., 2014; Momcheva et al., 2016). The UV $\beta$ slopes in the catalogue
are calculated from the $M_{\mathrm{UV}}-\beta$ relationship given in Bouwens
et al. (2009, 2015b).

For each survey used in this work, we scatter the mock JAGUAR photometry
within errors set by the $1\sigma$ measured depth with a minimum 10% flux
error in each filter before running through our selection procedure outlined
in § 2.5. The absolute magnitude, $M_{\mathrm{UV}}$, and UV continuum slope,
$\beta$, were calculated for the selected sample at the fixed EAZY redshifts
(as explained in § 3), with stellar masses calculated with Bagpipes for the
JAGUAR catalogue using the setup in Table 2 taking the same bandpass filters
and average depths from the CEERS survey.

To account for the SED modelling differences between BEAGLE and Bagpipes, we
calculate a median mass difference, $\langle\Delta
M_{\star}\rangle=\log_{10}(M_{\star,\mathrm{obs}}/M_{\star,\mathrm{int}})$,
between observed and intrinsic stellar masses in the CEERS JAGUAR catalogue
and apply these scaling factors to estimate observed stellar masses from the
JAGUAR catalogues simulating other EPOCHS-III fields. For our assumed Bagpipes
lognormal SFH, we find the differences between correctly identified high-$z$
SFGs, $\langle\Delta M_{\star,\mathrm{high-}z}\rangle=0.244$, and contaminant
Balmer break interlopers, $\langle\Delta
M_{\star,\mathrm{interlopers}}\rangle=0.527$. A closer look at the mass
differences for this sample using both the Bagpipes (Carnall et al., 2018) and
Prospector (Johnson et al., 2021) Bayesian SED fitting tools is presented in
EPOCHS-IV, although we note the large dependence of SFH assumption on the
stellar masses we measure.

For each survey, we calculate the contamination from the respective JAGUAR
catalogues following

$\mathrm{Cont}(\Theta_{\mathrm{obs}})=\frac{N_{\mathrm{selected,interlopers}}(\Theta_{\mathrm{obs}})}{N_{\mathrm{selected}}(\Theta_{\mathrm{obs}})}\mathrm{,}$
(1)

in bins of $\Theta=\\{(M_{\mathrm{UV}},\beta),(M_{\star},\beta)\\}$ in
intrinsic and observed galaxy property frames for completeness and
contamination respectively. In theory, the derived quantities in Equation 1
are also redshift dependent, however we decide not to bin the contamination by
redshift due to the relatively small size of our JAGUAR catalogues. The
contamination for our shallowest and deepest surveys (El Gordo and JADES-Deep-
GS, where a summary of survey depths is given in EPOCHS-I/II/IV and Table 1),
are shown in Figure 1.

Figure 1: Contamination from lower redshift SFGs in the El Gordo (upper
panels) and JADES-Deep-GS (lower panels) fields for the entire redshift range
used in this work ($7.5<z<13.0$ for NIRCam only fields and $6.5<z<13.0$
otherwise) as a function of $\Theta=(M_{\mathrm{UV}},\beta)$ (left panels) and
$\Theta=(\log_{10}M_{\star},\beta)$ (right panels) calculated from the JAGUAR
mock galaxy catalogue (Williams et al., 2018). The observed $\beta$ plotted
here is calculated using the photometric power law method. The total number of
selected low-z Balmer break interlopers and selected galaxies is shown in the
plot titles. As expected the most contamination occurs in the faintest, lowest
mass bins which becomes larger in the bluer bins due to the favouring
selection of redder galaxies at these intrinsic magnitudes/masses. JADES-Deep-
GS has less contamination in total which occurs in lower mass/fainter bins due
to its greater depth ($m_{\mathrm{F277W}}\simeq 30.2$) and increased number of
filters (HST/ACS_WFC/F606W and JWST/NIRCam/F335M) compared to El Gordo
($m_{\mathrm{F277W}}\simeq 29.0$).

We also use our JAGUAR catalogues to provide an estimate of the sample
completeness, calculated as

$\mathrm{Comp}(\Theta_{\mathrm{int}})=\frac{N_{\mathrm{selected}}(\Theta_{\mathrm{int}})}{N_{\mathrm{total}}(\Theta_{\mathrm{int}})}$
(2)

in the parameter spaces
$\Theta=\\{(z,M_{\mathrm{UV}},\beta),(z,M_{\star},\beta)\\}$, where we bin the
completeness in the redshift bins used in this work ($6.5<z<8.5$, $8.5<z<11$,
and $11<z<13$). We plot both $20\%$ completeness contours (red) and our JAGUAR
simulation limits (black) in Figure 11 and Figure 13. It is noteworthy that
although the effective sky area across our 5 JAGUAR realizations ($\sim
605~{}\mathrm{arcmin}^{2}$) far exceeds the sky coverage of our EPOCHS-III
galaxy sample, there still exists very few high mass galaxies in our highest
redshift bin. This is most-likely due to the rapid fall-off at the high-mass
end of the JAGUAR GSMF, which could be somewhat resolved by running more
JAGUAR realizations or using a simulation which more accurately match the high
mass end of updated JWST GSMFs. As well as this, the mock catalogues have a
limited range of $\beta$ covered almost entirety by our SED fitting template
set, making 20% completeness contours difficult to extend past mere faint, low
mass detection limits in the parameter spaces probed for this work.

##  3 Calculating UV properties

In this section we discuss the two methods of calculating $\beta$ from the
photometric data. The first method calculates $\beta$ via a power law fit
($f_{\lambda}\propto\lambda^{\beta}$) to the wideband photometric fluxes for
rest-frame wavelengths $1250<\lambda_{\mathrm{rest}}/\mathrm{\AA}<3000$. The
second fits the 10 C94 filters using the same power law to the Bayesian SED
template posterior for each source. These are discussed in § 3.1 and § 3.2
respectively. We compare these methods in § 3.3 and to spectroscopic $\beta$
in § 3.4. Since the Bayesian SED fitting method is systematically biased
against measuring the bluest $\beta$ due to the implicit prior on $\beta$, we
favour the photometric power law method in this work.

###  3.1 Calculating $\beta$ from photometric fluxes

Since the launch of JWST, we now have deep, high resolution imaging in the
Near InfraRed (NIR) redwards of $1.6~{}\mu\mathrm{m}$, making the direct
calculation of $\beta$ using multi-band photometry (as opposed to SED fitting)
the preferred method. Where previously we could not calculate $\beta$ in this
way at $z\gtrsim 8.5$ using only the reddest F160W HST/WFC3IR filter, the
vastly improved rest-frame UV coverage from the JWST/NIRCam F200W, F277W, and
F356W widebands now makes this possible at these high redshifts.

In this work, we fit a power law of the form
$f_{\lambda}\propto\lambda^{\beta}$ to the rest-frame UV photometry, as done
by e.g. Rogers et al. (2014); Bouwens et al. (2014) in the HST era and more
recently by e.g. Cullen et al. (2023a); Topping et al. (2022) using early JWST
ERS/ERO NIRCam data. We define which filters are used in this fitting
procedure as those which fall entirely within
$1250<\lambda_{\mathrm{rest}}~{}/~{}\mathrm{\AA}<3000$ as given by the same
50% flux limits as used in our selection criteria in § 2.5, where the rest
frame wavelengths are derived from the best-fitting EAZY-py redshift. For
filter $i$, this is calculated as

$f_{\lambda,i}\propto\frac{\int_{0}^{\infty}T_{i}(\lambda)\lambda^{1+\beta}\mathrm{d}\lambda}{\int_{0}^{\infty}T_{i}(\lambda)\lambda\mathrm{d}\lambda}$
(3)

in the photon-counting convention, where $T_{i}$ is the filter transmission
profile for filter $i$ taken from the SVO filter profile service. The filters
used are plotted as a function of redshift and rest wavelength coverage in the
top panel of Figure 2.

From these power law fits, we calculate $M_{\mathrm{UV}}$ from the bandpass
averaged flux in a $100~{}\mathrm{\AA}$ wide top-hat filter centred on
$\lambda_{rest}=1500~{}\mathrm{\AA}$. We correct these $M_{\mathrm{UV}}$ by a
factor
$-2.5\log_{10}(\mathrm{``FLUX\\_AUTO^{\prime\prime}/``FLUX\\_APER^{\prime\prime}})$
in the band with effective wavelength closest to
$\lambda_{\mathrm{obs}}=1500(1+z)~{}\mathrm{\AA}$ to account for extended
sources when required. Occasionally SExtractor over-corrects faint extended
sources, so in cases where “FLUX_AUTO”/“FLUX_APER”$>10$, we instead switch to
using the F444W band to apply the correction which avoids the issue.

In the lower panel of Figure 2 we plot the mean $\beta$ errors,
$\langle\sigma_{\beta}\rangle$, as a function of redshift for the entirely of
the EPOCHS-III sample (black) and in 3 magnitude bins
($M_{\mathrm{UV}}<-20.5$, $-20.5<M_{\mathrm{UV}}<-19.5$, and
$M_{\mathrm{UV}}>-19.5$). Here we outline two main reasons for the trend in
$\sigma_{\beta}$ across our sample: the number of JWST/NIRCam bands present in
the rest-frame UV and the SNR ratio of each galaxy in these aforementioned
filters. From Figure 2 we see that for the majority of our redshift range of
interest ($7.2<z<9.4$ and $9.5<z<12.3$) there are only two rest-UV widebands,
whereas this increases to three at $6.5<z<7.2$ and $12.3<z<13$. In the JADES-
Deep-GS data, the inclusion of the F335M medium band instead gives this
dataset three rest-UV filters at $10.8<z<12.3$ and four at $12.3<z<13$. It is
clear to see that the increased number of rest-UV bands reduces
$\sigma_{\beta}$ for our entire sample from $\sigma_{\beta}\simeq 0.4$ at
$z<7.2$ to $\sigma_{\beta}\simeq 0.55$ at $7.2<z<9.4$ and
$\sigma_{\beta}\simeq 0.47$ at $z>9.4$. We also note that objects with
brighter apparent UV magnitudes have reduced $\sigma_{\beta}$ due to their
smaller photometric errors. This is apparent in the lower panel of Figure 2,
where the average for the full sample with $\langle
M_{\mathrm{UV}}\rangle=-19.51$ lies between the $-20.5<M_{\mathrm{UV}}<-19.5$
and $M_{\mathrm{UV}}>-19.5$ bins.

Figure 2: Top: Rest-frame UV coverage of the JWST/NIRCam filters used in this
work as a function of both redshift and $\lambda_{\mathrm{rest}}$. The right-
hand side shows the rest-wavelength coverage of the 10 C94 filters in shaded
lime green which are used to avoid prominent rest frame UV nebular emission
lines shown as dash-dotted dark green lines in this same side-plot. Bottom:
Average power law measured $\beta$ errors, $\langle\sigma_{\beta}\rangle$ as a
function of both redshift and $M_{\mathrm{UV}}$. We plot our EPOCHS-III
sample, with $\langle M_{\mathrm{UV}}\rangle=-19.51$, in black as well as in 3
$M_{\mathrm{UV}}$ bins ($M_{\mathrm{UV}}<-20.5$ (purple),
$-20.5<M_{\mathrm{UV}}<-19.5$ (pink), and $M_{\mathrm{UV}}>-19.5$ (yellow)).
As expected, we observe a trend of decreasing $\sigma_{\beta}$ with increasing
intrinsic UV magnitude due to the increased apparent magnitude (and hence
decreased photometric errors) of the brighter sources in our sample.

###  3.2 Calculating $\beta$ via Bayesian SED fitting with Bagpipes

In addition to measuring $\beta$ with a power-law fit to the rest-frame UV
photometry, we run the Bagpipes Bayesian SED fitting code (Carnall et al.,
2018) to measure $\beta$ from the best-fit SED to all of the available
photometric fluxes rather than just those in the rest-frame UV. We use the
2016 version of the Bruzual & Charlot (2003) (BC03) SPS models assuming a
Kroupa (2001) IMF and Calzetti et al. (2000) dust attenuation law with
redshift fixed to that measured by EAZY-py. The Inoue et al. (2014) model for
InterGalactic Medium (IGM) attenuation is assumed and we do not model the
Ly$\alpha$ damping wing due to the uncertain nature of the patchy IGM during
reionization. We assume a lognormal SFH parametrized by the timescales for SF
onset ($t_{\mathrm{start}}$), SF peak ($t_{\mathrm{max.}}$), and lognormal
FWHM. The stellar mass ($M_{\star}$), metallicity ($Z$), V-band dust
attenuation ($A_{V}$), nebular ionization parameter ($U$), and Lyman continuum
escape fraction ($f_{\mathrm{esc,LyC}}$), are given wide priors in base 10
logarithmic space (equivalent to uniform priors on the logged parameter). We
fix the age of the nebular birth clouds to $t_{\mathrm{BC}}=10~{}\mathrm{Myr}$
which are also assumed to have no additional dust attenuation (as in the
default Bagpipes setup). A summary of the fixed and free parameters of our
Bagpipes fits, as well as the underlying assumptions made about the family of
SEDs we are using to perform the analysis, are given in Appendix B.

Table 2: Summary of fixed and free parameters used in our Bagpipes Bayesian SED fitting procedure. Parameter names, descriptions, and prior distributions/limits for the parameters defining the stellar, nebular and dust properties, as well as the lognormal star formation history are outlined. We note that we have recomputed the default Bagpipes BC03 SPS grids using CLOUDY v17.03 to extend the range of $\log U$ to $-1$. Parameter | Prior | Limits/Value | Description
---|---|---|---
Redshift
$z$ | Fixed | $z_{\mathrm{phot}}$ |  EAZY-py photo-z
Star formation history (lognormal)
$t_{\mathrm{start}}$ | $\log_{10}$ | $(1~{}\mathrm{Myr},t_{U})$ | Time of SF onset
$t_{\mathrm{peak}}$ | $\log_{10}$ |  $(10^{-3},15)$ Gyr | Time of peak SF
FWHM | $\log_{10}$ |  $(10^{-3},15)$ Gyr | SFH FWHM
Stellar properties
$\log_{10}(M_{\star}/\mathrm{M}_{\odot})$ | Uniform | $(5,12)$ | Stellar mass formed
$\log_{10}(Z_{\star})$ | Uniform | $(-6,1)$ | Stellar metallicity
$Z_{\mathrm{gas}}$ | Fixed | $Z_{\star}$ | Gas-phase metallicity
Nebular properties
$\log_{10}(U)$ | Uniform | $(-3,-1)$ | Ionization parameter
$f_{\mathrm{esc,LyC}}$ | $\log_{10}$ | $(0.001,1)$ | Lyman continuum
escape fraction
$t_{\mathrm{BC}}$ | Fixed |  $10$ Myr | Birth cloud lifetime
Dust properties
$A_{V}$ | $\log_{10}$ | $(10^{-4},10)$ | V-band attenuation
$\tau_{\mathrm{BC}}$ | Fixed | 0 | Birth cloud
optical depth
|  |  |

Several physical properties are calculated from these Bagpipes runs, including
$\beta$, $M_{\mathrm{UV}}$, and $M_{\star}$. $\beta$ is calculated by fitting
the same power law function as in § 3.1 to posterior spectrum in the 10 C94
top-hat filters and $M_{\mathrm{UV}}$ is calculated in the same way as in §
3.1. The posterior on stellar mass, $M_{\star}$, is an output from the base
Bagpipes code which is subsequently corrected for extended sources by
“FLUX_AUTO”/“FLUX_APER” in the F444W band where appropriate which traces the
rest-frame optical light.

We use the stellar masses calculated in our Bagpipes run when determining the
contamination likeliness of each galaxy in § 2.6 and when observing trends
between $\beta$ and $M_{\star}$ in § 5.3. Whilst we do not consider the
majority of the constrained Bagpipes parameters in this work, it is noteworthy
that the $f_{\mathrm{esc,LyC}}$ values in Figure 3 are illustrative only and
are not well constrained due to their degeneracy with stellar metallicities
and galaxy ages derived from the assumed lognormal SFH.

###  3.3 Comparison of photometric $\beta$ methods

In this section, we compare and contrast our two different photometric methods
to measure $\beta$ in order to determine which values to use in the analysis
performed in this work. A comparison of $\beta$ measurements (with the power
law $\beta$ corrected for the biases explored in § 4) is shown in the left-
hand panel of Figure 3.

The SED fitting method uses all available photometric data and produces more
precise values of $\beta$ for individual galaxies due to the mitigation of
photometric errors arising when fitting a power law to a small number of flux
measurements. In addition, using the SED method means that $\beta$ values are
calculated using the same underlying assumptions as the stellar masses,
leading to more consistent conclusions. In addition, the power law method is
more susceptible to photometric scatter (see § 4.4), rest-frame UV
contamination by nebular emission lines (see § 4.1), as well as the redshift
and filterset-dependent unequal coverage of the rest-frame UV (Figure 2).

Photometric SED fitting for $\beta$ is the favoured method of Finkelstein et
al. (2012); Bhatawdekar & Conselice (2021); Morales et al. (2023), however it
has previously been noted that this method is systematically biased by the
limited range of $\beta$ allowed by the choice of IMF and SPS model (e.g.
Rogers et al., 2013; Cullen et al., 2023a). This primarily impacts the range
of allowed blue $\beta$ slopes which biases the SED fitted $\beta$ redwards
(i.e. less negative values). Whilst our template set reaches as blue as
$\beta=-2.95$ for the youngest ($t\simeq 1$ Myr), most metal poor
($Z_{\star}=10^{-6}~{}Z_{\odot}$), dust-free, stellar dominated
($f_{\mathrm{esc,LyC}}=1$) galaxies, the likelihood of measuring blue $\beta$
slopes via SED fitting is impacted by the implicit $\beta$ prior (determined
as a complex function of our other input priors outlined in Table 2).

We plot this $\beta$ prior in the right-hand panel of Figure 3 and compare to
the $\beta$ posterior from our SED fitting method. It is clear that while the
bluest achievable $\beta$, shown as the black dash-dotted line, is
sufficiently blue, the probability of measuring $\beta\lesssim-2.7$ becomes
increasingly unlikely using the Bagpipes setup given in Table 2. We therefore
choose to adopt the power law photometric $\beta$ method in the rest of the
analysis done in this paper as to avoid this curtailing towards bluer $\beta$.
We test whether our $\beta$ prior is strongly dependent on SFH by re-running
our SED fitting procedure with the “continuity bursty” SFH parametrization
used in EPOCHS-IV and originally by Leja et al. (2019) and Tacchella et al.
(2022), finding minimal differences in the $\beta$ value at which this prior
falls off. Re-running instead using templates derived from BPASS as opposed to
BC03 pushes this soft limit bluewards by $\sim 0.1-0.2$ and the lowest
possible UV slope to $\beta=-3.1$; the value of these $\beta$ limits are
therefore impacted most strongly by choice of SPS model.

Figure 3: Left: Comparison of photometrically measured $\beta$ slopes for both
the SED fitting and bias corrected power law methods for our EPOCHS-III sample
colored by Bagpipes derived $f_{\mathrm{esc,LyC}}$.
$\beta_{\mathrm{PL,corr}}=\beta_{\mathrm{SED}}$ is shown in dashed orange,
which does not represent the best-fit to the data. Right: UV slope prior (red)
and posterior probability density (blue) for $\beta_{\mathrm{SED}}$. We see
that the posterior $\beta_{\mathrm{SED}}$ somewhat follows the prior in that
the extreme blue $\beta$ values are largely avoided due to the limited range
of parameter space explored by our Bagpipes template set. The black dot-dashed
line represents the bluest achievable $\beta_{\mathrm{SED}}=-2.95$.

###  3.4 Spectroscopic comparison

To assess the accuracy of our photometric $\beta$ measurements, we compare to
those derived from spectroscopic measurements using the low resolution $R\sim
100$ NIRSpec PRISM/CLEAR disperser-filter combination from the DAWN JWST
Archive (DJA).888https://dawn-cph.github.io/dja/ An outline of the data
reduction process using the public msaexp999https://github.com/gbrammer/msaexp
tool is given in Heintz et al. (2023a).

We cross-match our photometric sample from CEERS and JADES-Deep-GS with those
with robust spectroscopic redshifts in DJA finding a total of 55 matches, with
41 having robust NIRSpec redshifts (selected using “grade == 3”). This sample
consists of 15 galaxies from CEERS, 7 from the CEERS DDT (PI: P. Arrabal-Haro,
PID: 2750, Arrabal Haro et al., 2023a, b), 8 from JADES program 1180 (PI: D.
Eisenstein, PID: 1180, Eisenstein et al., 2023), and 11 from JADES program
1210 (PI: N. L$\ddot{\mathrm{u}}$tzgendorf, PID: 1210, Bunker et al., 2023). A
summary of the spectroscopically measured $\beta$ and previous literature
studies publishing these sources are given in Table 8.

Additionally, we search the DJA for medium ($R\sim 1,000$) and high resolution
($R\sim 2,700$) spectra, finding matches for 4 CEERS galaxies in the medium
resolution F100LP/G140M, F170LP/G235M, and F290LP/G395M filter-disperser pairs
(CEERSP3-559,2668,9946 and CEERSP6-7138). 4 JADES-Deep-GS galaxies have medium
resolution F070LP/G140M, F170LP/G235M and F290LP/G395M (IDs 12248, 15023,
15705, and 26579) grating spectra only and 6 more (IDs 14738, 15297, 18531,
18605, 19523, and 21391) have both medium resolution and additional high
resolution F270LP/G395H grating spectra.

We calculate spectroscopic $\beta$ slopes by fitting a power law in
$f_{\lambda}$ to the NIRSpec PRISM spectra both in the 10 C94 filters and in
the wavelength range $1250<\lambda_{\mathrm{{rest}}}~{}/~{}\mathrm{\AA}<3000$.
In addition, rest-frame UV signal-to-noise ratios (SNRs) are calculated over
this same wavelength range. We plot our spectroscopic $\beta$ measurements
against photometric $\beta$ calculated using both the power law (described in
§ 3.1) and SED fitting methods (see § 3.2) in Figure 4, including only those
galaxies with UV continuum SNR${}_{\mathrm{UV}}>3$.

Figure 4 also shows the comparison of spectroscopic and photometric redshifts
which, if significantly different, may provide a large impact on the measured
$\beta$. We categorize the photo-z error by a quantity $\epsilon$, where

$1+\epsilon=\frac{1+z_{p}}{1+z_{s}}=\frac{\lambda_{\mathrm{p}}}{\lambda_{\mathrm{s}}}$
(4)

also indicates a wavelength shift between photometric and spectroscopic
observed frame SED features at $\lambda_{\mathrm{p}}$ and
$\lambda_{\mathrm{s}}$ respectively. We find that 40/41 galaxies in this
spectroscopic sample fall within $\lvert\epsilon\rvert<0.1$, hence we conclude
that this is not a significant factor in the differences between $\beta$
measurements here.

We measure the inverse-variance weighted mean and standard error for our
different $\beta$ measurements, finding photometric results of
$\langle\beta_{\mathrm{phot,SED}}\rangle=-2.31\pm 0.00$ and
$\langle\beta_{\mathrm{phot,PL}}\rangle=-2.40\pm 0.02$. In addition, we also
correct the power law $\beta$ measurements for the $\beta$ bias explored in §
4 which reddens the inverse-variance weighted mean by $0.04$ to
$\langle\beta_{\mathrm{phot,PL,corr}}\rangle=-2.36\pm 0.02$. To assess whether
our photometric and spectroscopic $\beta$ measurements are statistically
likely to arise from the same underlying distribution (which they should since
they are the same galaxies), we perform a two-sided Kolmogorov-Smirnoff (KS)
test using the “scipy.stats.ks_2samp” (Hodges, 1958) python function. For our
power law photometric $\beta$ we measure $\mathrm{KS}=0.18$ with a
corresponding $\mathrm{p}=0.67$ after bias corrections compared to
$\mathrm{(KS,p)}=(0.29,0.11)$ for the SED fitting photometric $\beta$,
suggesting that the SED fitting results likely contain larger systematic
errors.

Since the SED fitting method uses the same wavelength coverage as the
spectroscopic results, we may naively expect these results to match closely,
however they remain systematically red (with $6.8\sigma$ significance) due to
the limited parameter space covered by the SED template set. This is portrayed
by the Bagpipes $\beta$ prior in the right-hand panel of Figure 3. Even though
the photometric power law fits for $\beta$ do not have complete coverage over
$1250<\lambda_{\mathrm{rest}}~{}/~{}\mathrm{\AA}<3000$ (with redshift and
filter-set dependence shown in Figure 2), they disagree at just the
$2.0\sigma$ significance level with the spectroscopic results. In addition, we
note that the difference could be compounded by additional $\beta$ biases
caused by strong rest-frame UV emission lines, Ly$\alpha$ emission, or damped
Ly$\alpha$ (DLA) systems that we do not correct for in this work (see § 4).

In the rest of this work, we favour the power law photometric $\beta$
calculation method due to the fewer systematic biases associated with
measuring $\beta$ in this way, albeit at the expense of greater photometric
$\sigma_{\beta}$ errors.

Figure 4: Comparison of redshifts and $\beta$ slopes for our $6.5<z<13$ EPOCHS
photometric sample cross-matched with $>3\sigma$ UV continuum detected
PRISM/CLEAR NIRSpec spectroscopic results from the DAWN JWST Archive (DJA).
Top: Power law measured $\beta$ that are both corrected (grey circles) and
uncorrected (hollow red circles) for the $\beta$ bias in § 4 against
spectroscopic $\beta$ measured over
$1250<\lambda_{\mathrm{rest}}~{}/~{}\mathrm{AA}<3000$. The blue square shows
the inverse-variance weighted mean for our bias corrected results. Bottom: SED
measured $\beta$ from our Bagpipes fitting against spectroscopic $\beta$
measured in the 10 C94 filters. Inverse-variance weighted mean values are
given in the upper left and a photometric vs spectroscopic redshift comparison
is given in the lower right. $\beta_{\mathrm{phot}}=\beta_{\mathrm{spec}}$ is
shown as a thick dashed line. Two-sided KS-test results and corresponding
p-values are shown in the lower left of both panels, where power law $\beta$
in the upper panel has been bias corrected.

##  4 UV continuum slope biases

###  4.1 The impact of strong rest-frame UV emission lines

It is well known that modelling the rest-frame UV as a power law is not
strictly correct due to the presence of the $2175~{}\mathrm{\AA}$ dust bump
(Stecher & Donn, 1965; Hensley & Draine, 2023) from carbonaceous dust grains
(see Witstok et al. (2023a) for more details), as well as rest-frame UV
emission lines and nebular continuum emission in SFGs. In this section, we
simulate the impact of rest-UV emission lines on $\beta$. We produce $10,000$
mock power law SEDs evenly spaced in redshift between $6.5<z<13$ at fixed
intrinsic $\beta_{\mathrm{int}}=-2.5$, with photometric errors calculated
assuming the galaxy has $m_{\mathrm{UV}}=26$ with each filter having a
$5\sigma$ depth of $m_{\mathrm{AB}}=30$. Doppler broadened rest-frame UV
emission lines with fixed Doppler parameter, $b=150~{}\mathrm{kms}^{-1}$, and
rest-frame equivalent widths (EWs) $1-25$ Å were individually added to the
SEDs before calculating bandpass-averaged fluxes in the standard 8 PEARLS
JWST/NIRCam bands. From these photometric fluxes, we then measured $\beta$
using our preferred power law method at the fixed input redshift with the bias
in $\beta$ measured as
$\Delta\beta_{\mathrm{line}}=\beta_{\mathrm{line}}-\beta_{\mathrm{int}}$. We
do this for the rest-UV emission lines C IV $\lambda 1549$, He II $\lambda
1640$, O III] $\lambda 1665$, and C III] $\lambda 1909$ which contaminate the
F115W, F150W, and F200W SW NIRCam wideband filters. The resulting bias as a
function of input redshift is shown in Figure 5 for
$\mathrm{EW}_{\mathrm{rest}}=10$ Å. We note that this simulation is
independent of $\beta$, $m_{\mathrm{UV}}$, and Doppler $b$, and that the
biases scale linearly with rest-frame EW such that

$\Delta\beta(\mathrm{EW}_{\mathrm{rest}})=\mathrm{EW}_{\mathrm{rest}}\frac{\Delta\beta(\mathrm{EW}_{\mathrm{rest}}=10~{}\mathrm{\AA})}{10}\mathrm{.}$
(5)

Previous observations of SFGs from the ground have found rest-UV emission line
equivalent widths (EWs) as high as 27 Å for the semi-forbidden C III] $\lambda
1909$ (Rigby et al., 2015; Stark et al., 2017) and $\sim 8-10$ Å for O III]
$\lambda 1665$ and He II $\lambda 1640$ (Nanayakkara et al., 2019; Saxena et
al., 2020) at $z\sim 4$. Strong He II $\lambda 1640$ emitters are often
associated with AGN when C IV $\lambda 1549$ emission is also present (e.g.
Feltre et al., 2016), although it has been shown that BPASS SPS models can
explain the observed ratios through star-formation alone (Xiao et al., 2018).
In addition, Jaskot & Ravindranath (2016) find that the BPASS models show that
C III] $\lambda 1909$ is the strongest line bluewards of $2700$ Å, which
becomes more prominent in young, metal-poor systems with high ionization
parameters.

More recently, several authors have reported strong rest-frame UV emission
lines at $z>6$ in NIRSpec spectroscopic data. At $z\simeq 8-8.5$, Tang et al.
(2023) find C III] $\lambda 1909$ EWs as large as $12-16~{}\mathrm{\AA}$, with
Hsiao et al. (2023) and Fujimoto et al. (2023) finding highly ionized bubbles
at $z=8.5-13$ and $z=10.17$ in the Ultradeep NIRSpec and NIRCam Observations
before the Epoch of Reionization survey (UNCOVER; PIs: I. Labbe & R. Bezanson;
PID: 2561; Bezanson et al., 2022) and the MACS-0647 lensing cluster. Recent
JWST observations have also detected strong C IV $\lambda 1549$ due to the
presence of hard ionizing fields at early times, with Topping et al. (2024)
observing rest-frame C IV $\lambda 1549$ EW of $34~{}\mathrm{\AA}$ in a
gravitationally lensed $z=6.11$ object behind the RXCJ2248 lensing cluster
(see also Mainali et al., 2017; Schmidt et al., 2017). Castellano et al.
(2024) also find the strongest EoR C IV $\lambda 1549$ EW of $45.8\pm
1.2~{}\mathrm{\AA}$ in the spectroscopically confirmed GHZ2/GLASS-z12 object
at $z=12.34$; this source was initially detected in GLASS NIRCam photometry by
Naidu et al. (2022) and Castellano et al. (2022) and is included in our
EPOCHS-III sample (see EPOCHS-I, Figs. 2/3) and provides evidence for UV line
emission bias in our sample. Additionally, Charbonnel et al. (2023); Cameron
et al. (2023a); Senchyna et al. (2023) also observe strong Nitrogen abundance
in the well-studied GN-z11 (Oesch et al., 2016; Bunker et al., 2023).

Figure 5: $\beta$ bias from rest-frame UV line emission. Shown are results
assuming rest-frame EWs of $10~{}\mathrm{\AA}$ for C IV $\lambda 1549$, He II
$\lambda 1640$, O III] $\lambda 1665$, and C III] $\lambda 1909$ applied to a
set of $10,000$ power law SEDs with $\beta_{\mathrm{int}}=-2.5$ evenly spaced
across $6.5<z<13$. These biases scale linearly with emission line EW and are
true under the assumption that $\epsilon=0$, or that the emission line does
not impact the photometrically derived redshift.

While we cannot determine rest-frame UV line EWs in our sample using broadband
photometry alone, these studies give us an idea of the expected $\beta$ biases
and ionization conditions in reionization era galaxies. The rest-UV line EWs
and intrinsic $\beta$ slopes measured in the C94 filters for the blue Larson
et al. (2023) templates used in this work are given in Table 3. Due to the
strong rest-frame UV emission lines in these templates, there is a redshift
and filter-set dependent bias between the power law fit to the photometry and
from the SED in the C94 filters. Even when accounting for the emission lines,
we find that using the C94 filters biases $\beta$ red by up to 0.1 compared to
our power law method. We attribute this to the increased wavelength coverage
of the C94 filters towards the nebular continuum two-photon turnover in the
youngest Larson et al. (2023) templates quantified by
$\Delta\beta_{\mathrm{neb}}=\beta_{\mathrm{C94}}(f_{\mathrm{esc,LyC}}=0)-\beta_{\mathrm{C94}}(f_{\mathrm{esc,LyC}}=1)$.

$\log_{10}$(Template Age / yr) | 6.0 | 6.5 | 7.0
---|---|---|---
EWrest(Ly$\alpha$) / Å | 0.0 | 0.0 | 0.0
EWrest(C IV $\lambda 1549$) / Å | 5.6 | 2.0 | 0.5
EWrest(He II $\lambda 1640$) / Å | 0.3 | 1.0 | 1.1
EWrest(O III] $\lambda 1665$) / Å | 5.2 | 3.6 | 1.0
EWrest(C III] $\lambda 1909$) / Å | 28.7 | 21.5 | 7.3
$\beta_{\mathrm{C94}}(f_{\mathrm{esc,LyC}}=0)$ | -2.36 | -2.32 | -2.52
$\beta_{\mathrm{C94}}(f_{\mathrm{esc,LyC}}=1)$ | -3.11 | -2.82 | -2.77
$\Delta\beta_{\mathrm{neb}}$ | 0.75 | 0.50 | 0.25
|  |  |
Table 3: Rest-frame EWs of the most prominent UV emission lines from the blue
Larson et al. (2023) templates used in our SED fitting procedure (set 4) for
ages ranging $1-10$ Myr. There is no Ly$\alpha$ present in these templates and
C III] $\lambda 1909$ is consistently the strongest of these lines which
appears especially potent at the youngest ages. Intrinsic $\beta$ slopes
calculated in the 10 C94 filters for $f_{\mathrm{esc,LyC}}=\\{0,1\\}$ (i.e.
including and excluding CLOUDY nebular emission) are also given. For a
decreasing galaxy age, the reddening of the underlying continuum by nebular
emission ($\Delta\beta_{\mathrm{neb}}$) increases in addition to the stellar
continuum becoming bluer.

###  4.2 Proximate damped Ly$\alpha$ systems and the Ly$\alpha$ damping wing

When traversing the IGM in the EoR, Ly$\alpha$ photons from the source galaxy
are absorbed by neutral hydrogen along the line-of-sight to the observer.
This, combined with attenuation from the Ly$\alpha$ forest bluewards of 1216
Å, produces the characteristic red-skewed Ly$\alpha$ damping wing. The most
widely used prescription of the Ly$\alpha$ damping wing is presented in
Miralda-Escudé (1998). This model assumes that the IGM has a constant (volume-
averaged) neutral hydrogen fraction, $x_{\mathrm{HI}}$, between the source
redshift, $z_{\mathrm{gal}}$, and that the end of reionization occurs at
$z_{\mathrm{Re,end}}$, and neglects the patchy nature of reionization. The
Ly$\alpha$ optical depth through the neutral IGM,
$\tau_{\mathrm{IGM,Ly\alpha}}(\lambda_{\mathrm{obs}},z_{\mathrm{gal}},x_{\mathrm{HI}},R_{\mathrm{b}},\tau_{\mathrm{GP}},z_{\mathrm{Re,end}})$101010It
is noteworthy that $z_{\mathrm{Re,end}}$ has a minor impact on calculated
values of $\tau_{\mathrm{IGM,Ly\alpha}}$ in the Miralda-Escudé (1998) model.,
is also given as a function of $R_{\mathrm{b}}$, the radius of the surrounding
ionized HII bubble, and $\tau_{\mathrm{GP}}$, the cosmology-dependent Gunn &
Peterson (1965) optical depth. Using NIRSpec data from the first year of
_JWST_ operation, large ionized bubbles with $R_{\mathrm{b}}\gtrsim 1$ cMpc
have been inferred from these Ly$\alpha$ damping wings at high-redshift
(Whitler et al., 2023; Witstok et al., 2023b; Umeda et al., 2023), with a
large implied contribution from faint galaxies with $M_{\mathrm{UV}}>-16$ in
large-scale galaxy overdensities (Lu et al., 2024). Recent comparisons to
theorized damping wing profiles assuming a more realistic patchy reionization
process (e.g. Keating et al., 2023) suggest that the presence of small amounts
of residual HI as low as $x_{\mathrm{HI}}\sim 10^{-5}$ (McQuinn, 2016) within
these ionized bubbles and strong Ly$\alpha$ emission (seen observationally by
Saxena et al., 2023a) may instead be responsible for these large
$R_{\mathrm{b}}$ values. Due to the uncertainties in modelling this damping
wing, we choose not to quantify any associated $\beta$ biases. We do note,
however, that this is likely to bias our $\beta$ slightly redwards due to the
softening of the Ly$\alpha$ break leading to minor redshift overestimations
when SED fitting using templates that do not include this dampening.

In addition, DLA systems with high column densities of neutral hydrogen
$N_{\mathrm{HI}}>2\times 10^{20}~{}\mathrm{cm}^{-2}$ (Lanzetta, 2000), arising
from dense gas clouds in the vicinity of starburst galaxies, have been
observed at high-redshift both in quasar (Totani et al., 2006) and galaxy
spectra (Heintz et al., 2023a). These systems act to soften the Lyman-break
and increase photometric Lyman-Balmer break degeneracy, decreasing photo-z
accuracy. The UV $\beta$ slopes of these systems are often biased red due to
the reduction in flux associated with the first band redwards of
$1250~{}\mathrm{\AA}$. In addition, photo-z’s are often overestimated when SED
fitting codes confuse the reduction in flux in the first band redwards of the
Lyman-$\alpha$ break with the break itself. At certain filter-set dependent
redshifts, $\Delta\beta_{\mathrm{DLA}}$ becomes smaller with increasing
$N_{\mathrm{HI}}$ when this band is no longer included in the rest-frame UV
range used to calculate $\beta$. We do not attempt to identify DLAs in our
selection procedure, but instead simulate the impact on simple power law SED
models.

We produce $10,000$ mock power law SEDs with $\beta=-2.5$ and
$m_{\mathrm{UV}}=26$ are produced with redshifts evenly spaced in the interval
$6.5<z<13$ assuming Inoue et al. (2014) IGM attenuation post-reionization.
These SEDs are then propagated through the ISM with fixed Ly$\alpha$ velocity
offset, $\Delta v_{\mathrm{Ly}\alpha}=0$, Doppler $b=150~{}\mathrm{kms}^{-1}$,
and a range of neutral hydrogen column densities
$\log_{10}(N_{\mathrm{HI}}/\mathrm{cm}^{-2})=\\{21,21.5,22,22.5,23,23.5\\}$.
The DLA optical depth is calculated as

$\tau_{\mathrm{DLA}}(\lambda_{\mathrm{rest}})=CaN_{\mathrm{HI}}H(a,x(\lambda_{\mathrm{rest}}+\Delta
v_{\mathrm{Ly\alpha}},b))\mathrm{,}$ (6)

where $C$ and $a$ are the Ly$\alpha$ photon absorption constant and damping
parameters, and $H(a,x(\lambda,b))$ is the Voijt-Hjerting approximation to the
Voigt profile (Tepper-García, 2006) which is dependent on the Doppler
parameter, ${b=\sqrt{2k_{B}T_{\mathrm{HI}}/m_{p}}}$, via

$x(\lambda,b)=\frac{c(\lambda-\lambda_{\mathrm{Ly\alpha}})}{b\lambda_{\mathrm{Ly\alpha}}}\mathrm{.}$
(7)

This Voigt profile is insensitive to the choice of $b$ for DLAs with large
$N_{\mathrm{HI}}$ where the Lorentzian wings from the naturally broadened
Ly$\alpha$ dominate.

We then calculate bandpass-averaged fluxes from the standard 8-band PEARLS
JWST/NIRCam filterset and re-run through our EAZY-py photo-z fitting procedure
before re-calculating $\beta$ using the photometric power law method at the
derived photo-z to calculate the bias
$\Delta\beta_{\mathrm{DLA}}=\beta_{\mathrm{DLA}}-\beta_{\mathrm{int}}$. It is
noteworthy that in this process we include only photometric bands at
$\lambda_{\mathrm{rest}}<3000$ Å to avoid the non-power law nature of the
Balmer jump at $3646$ Å, Balmer break at $\sim 4000$ Å, and strong rest-
optical nebular emission lines including [OIII]/H$\beta$.

We plot $\Delta\beta_{\mathrm{DLA}}$ as a function of $N_{\mathrm{HI}}$, in
Figure 6, finding $\Delta\beta_{\mathrm{DLA}}=0.5$ even in the most extreme
DLA scenarios with $N_{\mathrm{HI}}=10^{23.5}~{}\mathrm{cm}^{-2}$ , meaning
that our sample may well be biased red by DLAs in some specific cases. This
upper $N_{\mathrm{HI}}$ limit is set by the required HI column density should
the Cameron et al. (2023b) objects be contrived DLAs as opposed to nebular
dominated systems. As of now, the highest DLA column densities observed in the
early Universe ($z>10$) are approximately
$\log_{10}(N_{\mathrm{HI}}~{}/~{}\mathrm{cm}^{-2})=22-22.5$ (Heintz et al.,
2023a; Umeda et al., 2023; D’Eugenio et al., 2023), and large spectroscopic
studies have recently been performed to assess the abundance of these DLA
systems (Heintz et al., 2024). Much is still to be determined about these
systems, however, including whether large neutral gas reservoirs can form in
the presence of strong winds, producing larger Ly$\alpha$ velocity offsets
$\Delta v_{\mathrm{Ly\alpha}}\gtrsim+400~{}\mathrm{kms}^{-1}$, as well as how
common DLAs are at high-redshift where gas masses are expected to be larger.
We conclude that it is likely that some of our sources are biased red, however
the extent of this in our photometric sample is largely unknown.

Figure 6: Left: Power law measured $\beta$ bias as a function of input
redshift and DLA HI column density, $N_{\mathrm{HI}}$. Biases as large as
$0.5$ redwards are observed in DLAs with the largest column densities,
$N_{\mathrm{HI}}=10^{23.5}~{}\mathrm{cm}^{-2}$, at redshifts
$z\sim\\{6.5,9,12\\}$. Right: DLA $\beta$ bias as a function of photo-z error,
$\epsilon$, as defined in Equation 4, with red values showing overestimated
photo-z, reaching a maximum of $15\%$. The black region towards the highest
redshifts shows the region that would fail the observed $z<13$ cut due to the
photo-z overestimation.

###  4.3 The impact of Ly$\alpha$ emission

Although Lyman Alpha Emitters (LAEs) are known to exist at low-redshift,
during the EoR the observed Ly$\alpha$ rest-frame EWs are expected to reduce
due to absorption by neutral Hydrogen along the line-of-sight (i.e. the
Lyman-$\alpha$ damping wing). Although the exact nature of this damping wing,
and indeed the size of ionized bubbles in the early Universe, is largely
unknown, several studies have observed high EW Ly$\alpha$ systems at $z>6$ in
recent spectroscopic JWST data (e.g. Saxena et al., 2023a, b; Nakane et al.,
2023; Witstok et al., 2023b; Chen et al., 2024; Napolitano et al., 2024). In
this section, we study the impact of these LAEs on measured $\beta$ slopes in
JWST/NIRCam wideband studies. Even though both the number density (Haiman,
2002; Malhotra & Rhoads, 2006) and EW evolution (Tang et al., 2024) of these
EoR LAEs has already been studied, we do not quantify the average bias in our
sample as the selection efficiency of LAEs by our specific EPOCHS criteria in
§ 2.5 is not well known.

The impact of Ly$\alpha$ emission on photometrically measured $\beta$ slopes
has in fact already been studied by Rogers et al. (2013) who find a blue bias
that increases with EW. Their work analyses the bias in $\beta$ measured from
a power law fit to HST/WFC3IR F105W (Y105), F125W (J125), F140W (J140), and
F160W (H160) photometry at $6.5<z<7.5$, finding
$-0.8\lesssim\Delta\beta_{\mathrm{Ly\alpha}}\lesssim-0.5$ at $z\simeq 7$ for
rest-frame EWs of $50-100~{}\mathrm{\AA}$. In this section, we analyse the
bias in $\beta$ measured using JWST/NIRCam widebands which may be
systematically different from the $\beta$ measured using Y105J125J140H160
HST/WFC3IR photometry.

We calculate the impact of Ly$\alpha$ emission on the measured $\beta$ slopes
using a similar method to that used to calculate the DLA bias in § 4.2.
$10,000$ mock power law SED templates are created with fixed $\beta=-2.5$,
$m_{\mathrm{UV}}=26$ at redshifts $6.5<z<13$ for each rest-frame Ly$\alpha$
$\mathrm{EW}=\\{5,10,20,30,40,50,75,100,150,200,300\\}$ Å including Ly$\alpha$
forest IGM attenuation following the Inoue et al. (2014) prescription, with
the upper EW limit set by the largest observable values in Saxena et al.
(2023a). Bandpass-averaged fluxes are generated for the standard 8 JWST/NIRCam
filters used in PEARLS, with errors produced assuming a $5\sigma$ depth of
$m_{\mathrm{AB}}=30$, and the photometry is run through EAZY-py to determine
photo-z’s. $\beta$ is then calculated by the power law method using the rest-
frame UV filters determined by the EAZY-py photo-z and compared to the
intrinsic $\beta=-2.5$ to determine $\Delta\beta_{\mathrm{Ly\alpha}}$, with
results shown in Figure 7.

In order to account for the increased flux in the first redwards band of the
Ly$\alpha$ break due to this Ly$\alpha$ emission, SED fitting procedures with
underestimated Ly$\alpha$ emission typically underestimate photo-z’s (right
panel of Figure 7). This introduces the band containing Ly$\alpha$ into the
rest-frame UV and hence biases $\beta$ blue. As can be seen from Figure 7,
there are two redshift ranges of interest (left and central panels) where
Ly$\alpha$ emission may produce significant
$\Delta\beta_{\mathrm{Ly\alpha}}\simeq-0.6$. At $z=7.3$, the photo-z
underestimation induced by the Ly$\alpha$ emission causes the F090W filter to
be incorrectly identified as the dropout filter, with the highly-elevated
F115W flux density now strongly biasing $\beta$ blue. The scenario is similar
at $z=9.8$ instead with the F115W and F150W filters, although the effect is a
lot weaker than at $z=7.3$ due to the gap between the F115W and F150W filters
not present between F090W and F115W resulting in weaker Ly$\alpha$ throughput
and hence little impact on $\Delta\beta$. We note that this estimate of bias
is likely an upper limit due to the additional redshift constraints provided
by photometric data at $\lambda_{\mathrm{rest}}>3000$ Å.

Figure 7: Left & centre: $\beta$ bias as a function of rest-frame Ly$\alpha$
EW at input redshift $z\simeq 7.3$ (left panel) and $z\simeq 9.8$ (centre
panel). It is worth noting the colorbar scaling in these two panels, where the
bias at $z\simeq 7.3$ is far greater than at $z\simeq 9.8$ even before
accounting for the reduction in observed EW due to absorption in the neutral
IGM along the line-of-sight. Right: Photo-z error as a function of input
redshift and Ly$\alpha$ EW. We see that the photo-z is underestimated by as
much as $\epsilon=-0.15$ (i.e. 15%) in the case of strong Ly$\alpha$ emission
at $z\simeq\\{8.5,12.5\\}$.

The bias methodology mentioned above falls short of a complete implementation
of Ly$\alpha$ due to our neglect of Doppler broadening due to HI velocity
dispersion, $\sigma_{\mathrm{Ly\alpha}}$, galactic dynamics/outflows shifting
the major Ly$\alpha$ peak by $\Delta v_{\mathrm{Ly\alpha}}$ or doubly or
triply peaked emission. The profiles of LAEs are explored in more detail using
simulations (Blaizot et al., 2023), however this is not expected to have a
major impact on our bias simulation results. We also do not consider the
impact of the Ly$\alpha$ damping wing, the shape of which is not well known
due to the unknown volume averaged neutral fraction of HI as a function of
redshift, $\langle x_{\mathrm{HI}}(z)\rangle$, and the patchy reionization
process. A further discussion of the Ly$\alpha$ damping wing is presented in §
4.2.

###  4.4 Photometric error coupling bias

The most widely reported $\beta$ bias is the so-called “photometric error
coupling bias” (termed by Bouwens et al., 2012) which biases faint objects
approaching our detection limit blue. This bias impacts high redshift galaxy
samples which are commonly selected based on SNR requirements either side of
the Ly$\alpha$ break at $1216~{}\mathrm{\AA}$ or a specific Ly$\alpha$ break
strength, both of which favour the up-scattering of photometric measurements
at luminosities just below the survey detection limit. In contrast, galaxies
where the first photometric band at
$\lambda_{\mathrm{rest}}>1216~{}\mathrm{\AA}$ scatters down have a reduced
Ly$\alpha$ break strength and a more prominent Balmer break photo-z solution
and hence are less likely to be selected in our sample.

This bias has been analysed extensively for $\beta$ measurements using
HST/WFC3IR F125W, F140W and F160W fluxes in the 2009 Hubble Ultra Deep Field
(HUDF09; PI: G. Illingworth, HST PID: 11563, see Bouwens et al., 2011), HUDF12
(PI: R. Ellis, HST PID: 12498, Ellis et al., 2013; Koekemoer et al., 2013) and
CANDELS (Koekemoer et al., 2011; Grogin et al., 2011) campaigns at $z\sim 6-8$
(e.g. Finkelstein et al., 2012; Bouwens et al., 2012; Dunlop et al., 2012;
Rogers et al., 2013; Dunlop et al., 2013; Bouwens et al., 2014), as well as in
the HFF MACS-0416 lensing cluster (Bhatawdekar & Conselice, 2021). More
recently, this has been done for large photometrically selected JWST samples
using both power law (Cullen et al., 2023a, b) and SED fitting (Morales et
al., 2023) measurements of $\beta$.

To analyse this photometric error coupling bias we produce $200,000$ power law
SEDs for each intrinsic $\beta=\\{-1,-1.5,-2,-2.5,-3\\}$, which appropriately
match the UV slopes in our observed data, for each field used in this work.
Each SED is next randomly assigned an intrinsic $m_{\mathrm{UV}}$ in the
interval $26<m_{\mathrm{UV}}<30$ and bandpass-averaged fluxes are measured
with the appropriate filterset for the field. Photometric errors are
calculated exactly as for the real data, assuming each galaxy has local depth
given by the average $5\sigma$ depths in the field with a 10% minimum error
floor to reflect the same NIRCam ZP tolerance allowed in our real photometry.
We next scatter our photometric flux measurements within their respective
Gaussian errors and re-run through the same high-$z$ EAZY-py photo-z
calculation and selection as in § 2.4. As before, we re-measure $\beta$ using
the power law method and calculate $\Delta\beta$ for selected objects only.
Results of this simulation for the NEP-TDF field are shown as a function of
input $m_{\mathrm{UV}}$ in Figure 8 and redshift in Figure 9.

In Figure 8, we see that bluer objects are inherently more selectable with the
EPOCHS selection criteria; approximately 57% of mock galaxies with intrinsic
$\beta=-3$ are selected compared to 37% with $\beta=-1$, meaning redder faint
galaxies are more likely biased blue compared to their intrinsically bluer
counterparts in the NEP at a fixed input $m_{\mathrm{UV}}$. The reasoning for
this is two-fold, namely that faint red galaxies more regularly fail the
$5\sigma$ SNR cut in the first band entirely at
$\lambda_{\mathrm{rest}}>1216~{}\mathrm{\AA}$ (criteria 3 in § 2.5) and are
more likely confused with Balmer break galaxies at lower redshift (criteria 4
and 5 in § 2.5) due to their shallower Ly$\alpha$ breaks. This trend is
observed across all of the fields in this study.

In Figure 9, we observe both a red bias at $z\simeq 7.6,10.5$ and a blue bias
at $z\simeq 8.6,11.8$ caused mainly by photo-z inaccuracy. At these redshifts,
the Ly$\alpha$ break passes between the F090W/F115W and F115W/F150W bands,
leading to incorrect photo-z measurements (the second of these is also noted
in Cullen et al., 2023b). We note that these redshift-dependent effects can be
mitigated considerably by the inclusion of deep HST/WFC3IR and medium band
JWST/NIRCam photometry.

Figure 8: $\beta$ bias as a function of $m_{\mathrm{UV}}$ for the NEP-TDF
field calculated from $200,000$ power law SEDs with intrinsic
$\beta=\\{-1,-1.5,-2,-2.5,-3\\}$ shown in green, light blue, yellow, dark
blue, and red respectively. The number of selected SEDs is given in the legend
for each intrinsic $\beta$, with bluer $\beta$ (at fixed $m_{\mathrm{UV}}$)
more efficiently selected. The normalized number density as a function of
$m_{\mathrm{UV}}$ is plotted in the lower left inset. Since the redder SEDs
are less efficiently selected, they exhibit the greatest $\beta$ bias, which
can reach as large as $-0.85$ for $\beta=-1$, $m_{\mathrm{UV}}=29$ in the NEP.
Figure 9: Top: $\Delta\beta$ for galaxies with intrinsic $\beta=-3$ in the
NEP-TDF as a function of input redshift. The mock galaxies are colored by the
photometric redshift error, where blue(red) points have
underestimated(overestimated) photo-z’s respectively. The overplotted thick
black line shows the average bias in the selected sample, a zoomed in version
of which is shown in the lower right. Bottom: $\beta$ bias as a function of
$m_{\mathrm{UV}}$ in the NEP-TDF from $200,000$ power law SEDs with results
for intrinsic $\beta=\\{-1,-1.5,-2,-2.5,-3\\}$ shown in green, light blue,
yellow, dark blue, and red respectively. Median $\Delta\beta$ values,
$\langle\Delta\beta\rangle$, are shown as horizontal dashed lines which become
bluer with reddening intrinsic $\beta$.

###  4.5 Summary of bias correction factors

In this section, we have outlined potential photometric $\beta$ biases from UV
line emission, LAEs, and DLAs using power law SEDs and the 8 wideband PEARLS
NIRCam filters. These features, however, are not detectable in our EPOCHS-III
photometric data and we therefore do not correct for these biases in this
work. The $\beta$ correction factors used in this work are estimated based on
the redshift and $m_{\mathrm{UV}}$ dependence of the photometric error
coupling bias only. This is calculated for each individual galaxy following

$\Delta\beta_{\mathrm{tot}}=\Delta\beta_{z,i}+\Delta\beta_{m_{\mathrm{UV}},i}-\langle\Delta\beta\rangle_{i}\mathrm{.}$
(8)

Here $\Delta\beta_{\mathrm{tot}}$ is the total $\beta$ bias,
$\Delta\beta_{z,i}$ and $\Delta\beta_{m_{\mathrm{UV}},i}$ are the photometric
error coupling biases in terms of redshift and $m_{\mathrm{UV}}$ in field $i$,
and $\langle\Delta\beta\rangle_{i}$ is the median bias included to avoid
double counting. For our EPOCHS-III sample, we measure an average bias
measurement of $\langle\Delta\beta\rangle=-0.06^{+0.06}_{-0.10}$. While this
measurement seems like a minor correction, we note that at the extremes we
measure galaxies with $\Delta\beta=-0.55$ and $\Delta\beta=0.24$. It is also
clear from Figure 8 that the bias is greater for redder galaxies that are
fainter and have lower masses due to the preferential up-scattering of red
compared to blue galaxies at a given $m_{\mathrm{UV}}$. This means the bias
corrections made in this work will mostly act to redden already red galaxies,
increasing the slope of our measured $\beta-M_{\mathrm{UV}}$ and decreasing
the slope of our $\beta-M_{\star}$ relations in § 5.2 and § 5.3 respectively.

The use of power law SEDs, however, does not describe the full picture from
the real Universe. Additional constraints in the rest-frame optical (i.e.
nebular emission lines, the $\lambda_{\mathrm{rest}}=3646~{}\mathrm{\AA}$
Balmer jump, or the Balmer break at $\lambda_{\mathrm{rest}}\simeq
4000~{}\mathrm{\AA}$) will increase the accuracy of photo-z measurements. We
conclude, therefore, that the bias corrections as a result of photo-z
uncertainties are likely overestimated in all fields due to this necessary
simplification. In addition, we note that this technique neglects possible
degeneracies that may exist between $\Delta\beta_{z,i}$ and
$\Delta\beta_{m_{\mathrm{UV}},i}$, and that an even more computationally
expensive procedure is required to adequately measure these using a more
statistically significant sample size.

##  5 Results

###  5.1 Redshift evolution

We now determine whether the decreasing trend of $\beta$ with redshift found
in HST data (e.g. Finkelstein et al., 2012; Bouwens et al., 2014; Bhatawdekar
& Conselice, 2021) extends to $z>10$ in JWST data. We plot our bias corrected
power law $\beta-z$ relation for the EPOCHS-III sample in Figure 10, with
circular beige points showing bootstrapped median data points at $z\simeq 7$,
$z\simeq 9$ and $z\simeq 11.5$. For each of our $10,000$ bootstraps we 1)
randomly scatter each galaxy within their redshift and $\beta$ PDFs, 2) bin
the galaxies in redshift and randomly select (with replacement) a number of
galaxies chosen from a Gaussian centred on the number of galaxies in the bin
with width given by the Poisson error on this, 3) calculate median of each
bin, 4) calculate median and $16^{\mathrm{th}}-84^{\mathrm{th}}$ percentiles
of the bootstrapped median values in each bin. We also follow this method to
calculate bootstrapped median data points in § 5.2 and § 5.3. Our bootstrapped
median data points, as well as the corresponding median $\langle
M_{\mathrm{UV}}\rangle$ values are given in Table 4.

Table 4: Bootstrapped binned results of median $z$ and bias corrected power law $\beta$ for our EPOCHS-III sample as plotted in Figure 10. We show median $M_{\mathrm{UV}}$ values in each redshift bin with the errors representing the $16^{\mathrm{th}}-84^{\mathrm{th}}$ percentiles of the $M_{\mathrm{UV}}$ distribution. It is clear that at higher redshift our EPOCHS-III selection criteria selects on average intrinsically brighter sources. $z$ bin | $\langle z\rangle$ | $\langle\beta\rangle$ | $\langle M_{\mathrm{UV}}\rangle$
---|---|---|---
$6.5<z<8.5$ | $6.98^{+0.03}_{-0.04}$ | $-2.36^{+0.03}_{-0.03}$ | $-19.27^{+0.79}_{-0.75}$
$8.5<z<11$ | $8.97^{+0.04}_{-0.03}$ | $-2.50^{+0.11}_{-0.10}$ | $-19.54^{+0.51}_{-0.74}$
$11<z<13$ | $11.64^{+0.19}_{-0.10}$ | $-2.69^{+0.15}_{-0.15}$ | $-19.79^{+0.53}_{-0.56}$
|  |  |

We find that the individual data points in Figure 10 are best fit by the power
law

$\beta=-1.51\pm 0.08-(0.097\pm 0.010)\times z\mathrm{,}$ (9)

implying a decreasing $\beta$ at earlier cosmic times. The negative slope of
$\beta-z$ is steeper than the JWST photometric results of Topping et al.
(2023) and shallower than those of Cullen et al. (2023b), who measure
$\mathrm{d}\beta/\mathrm{d}z=-0.030^{+0.024}_{-0.029}$ and
$\mathrm{d}\beta/\mathrm{d}z=-0.28^{+0.04}_{-0.04}$ respectively at similarly
high redshifts, although it remains consistent with the
$\mathrm{d}\beta/\mathrm{d}z=-0.10\pm 0.06$ derived at $z<6$ by Bouwens et al.
(2014). These trends are dependent on both the method used to measure $\beta$
and the average $M_{\mathrm{UV}}$ of the galaxy sample. The differences may be
explained by the dependence of $\mathrm{d}\beta/\mathrm{d}z$ on the average
intrinsic UV brightness of the galaxy sample, as noted by Topping et al.
(2023), with brighter $\langle M_{\mathrm{UV}}\rangle$ samples showing a
steeper evolution. Our $\langle M_{\mathrm{UV}}\rangle=-19.35$ is brighter
than that of Topping et al. (2023), who measure $-18.61<\langle
M_{\mathrm{UV}}\rangle<-18.16$ in their $0.78-1.15~{}\mu\mathrm{m}$ dropout
samples from deep JADES data in GOODS-South. At $z\simeq\\{9.5,10.5,11.5\\}$,
Cullen et al. (2023b) measure $\langle M_{\mathrm{UV}}\rangle=\\{-18.9\pm
0.8,-19.5\pm 0.7,-19.1\pm 0.5\\}$ in their combined sample, although their
$\langle M_{\mathrm{UV}}\rangle=-21.2\pm 1.6$ at $7.5<z<9$ is significantly
brighter due to their use of much wider ground based COSMOS/UltraVISTA data
which drives their measured $\mathrm{d}\beta/\mathrm{d}z$ steeper as a result
of their redder $\langle\beta\rangle$ at $z\simeq 8-8.5$. The stacked
spectroscopic results of Roberts-Borsani et al. (2024) at $z>5$ yield
$\mathrm{d}\beta/\mathrm{d}z=-0.06\pm 0.01$, deviating from our results by
$\simeq 2.6\sigma$

We test the $M_{\mathrm{UV}}$ dependence of our derived
$\mathrm{d}\beta/\mathrm{d}z$ by splitting our sample into bright
($M_{\mathrm{UV}}<-19.5$) and faint ($M_{\mathrm{UV}}>-19.5$) sub-samples
before repeating our $\beta-z$ power law fitting procedure. In our faint
($\langle M_{\mathrm{UV}}\rangle=-18.94$) and bright ($\langle
M_{\mathrm{UV}}\rangle=-19.92$) sub-samples we measure
$\mathrm{d}\beta/\mathrm{d}z=-0.108\pm 0.015$ and
$\mathrm{d}\beta/\mathrm{d}z=-0.092\pm 0.013$ respectively. It is clear that
there exists no clear evolution in $\beta-z$ with $M_{\mathrm{UV}}$ in our
sample which we attribute to the flatter
$\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ found in this work compared to
Cullen et al. (2023b), Topping et al. (2023) and Bouwens et al. (2014). This
highlights the susceptibility of conclusions related to the evolution of
sample averaged $\beta$ to minor differences in SFG selection criteria between
different studies. A further discussion regarding trends with
$M_{\mathrm{UV}}$ is presented in § 5.2.

Our derived $\beta-z$ slope implies that, over cosmic time, there is an
increasing metal enrichment and dust content within galaxies. We compare our
results to the representative UV continuum slope of dwarf galaxies enriched by
Pop. III stars given as $\langle\beta\rangle=-2.51\pm 0.07$ by Jaacks et al.
(2018). Our Equation 9 suggests that this transition occurs at a redshift
$z=10.3$, although we note that this is later than suggested by the $\beta-z$
fits of Cullen et al. (2023b) and Topping et al. (2023) who measure $z=11.1$
and $z=14.3$ respectively.

Figure 10: Bias corrected power law $\beta$ evolution as a function of $z$.
Black background points show measurements for individual galaxies and beige
circular points show our bootstrapped median points. The
$16^{\mathrm{th}}-84^{\mathrm{th}}$ percentiles of the posterior power law fit
to the data, as given in Equation 9, is shown in blue. We compare our results
to the observational HST studies of Bouwens et al. (2014) (blue stars),
Wilkins et al. (2016) (lime green squares), Bhatawdekar & Conselice (2021)
(orange squares) and more recent JWST work by Nanayakkara et al. (2023) (red
squares), Cullen et al. (2023b) (purple triangles), Topping et al. (2023)
(dark green triangles), and Roberts-Borsani et al. (2024) (yellow triangles).

###  5.2 Correlations with $M_{\text{UV}}$

We next analyse the trends of observed $\beta$ with $M_{\mathrm{UV}}$, with
$\beta$ calculated using both methods explained in § 3. We correct our power
law $\beta$ results using the average correction factors outlined in § 4.4
based on interpolated $m_{\mathrm{UV}}$, $\beta$, and $z$. Results using bias
corrected power law $\beta$ measurements are shown in Figure 11, where we
split the galaxies into our three redshift bins ($6.5<z<8.5$, $8.5<z<11$, and
$11<z<13$). We calculate bootstrapped median $\langle M_{\mathrm{UV}}\rangle$
and $\langle\beta\rangle$ for our sample (in 2, 3, and 4 magnitude bins
respectively for the $6.5<z<8.5$, $8.5<z<11$, and $11<z<13$ redshift bins)
following the same procedure introduced in § 5.1, with $\langle\beta\rangle$
and $\langle M_{\mathrm{UV}}\rangle$ for each bin given in Table 5.
Additionally, we quantify the amplitude and slope of these relations by
fitting a power law of the form

$\beta=\frac{\mathrm{d}\beta}{\mathrm{d}M_{\mathrm{UV}}}(M_{\mathrm{UV}}+19)+\beta(M_{\mathrm{UV}}=-19)$
(10)

to each individual data point in our three redshift bins to allow direct
comparison with recent JWST studies on $\beta$ (e.g. Topping et al., 2023;
Cullen et al., 2023b). The fitted amplitude at $M_{\mathrm{UV}}=-19$,
$\beta(M_{\mathrm{UV}}=-19)$, and slope,
$\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$, for fits using SED fitting
$\beta$ as well as both bias corrected and uncorrected power law $\beta$ are
given in Table 7 and plotted in relation to other observation results in
Figure 12.

We compare our results with binned observational studies from both HST
(Finkelstein et al., 2012; Dunlop et al., 2012; Bouwens et al., 2012; Dunlop
et al., 2013; Bouwens et al., 2014; Bhatawdekar & Conselice, 2021) and JWST
(Topping et al., 2023; Cullen et al., 2023b; Roberts-Borsani et al., 2024) in
the left hand panels of Figure 11. In addition, we also plot individual
galaxies from both Wilkins et al. (2016) and Morales et al. (2023) as starred
points. In addition to observational comparisons, we also compare to a wide
range of simulated results from THESAN (Kannan et al., 2022, $z=\\{7,9\\}$),
SC SAM GUREFT (Yung et al., 2023, 2024, $z=\\{7,9,11\\}$), DELPHI (Mauerhofer
& Dayal, 2023, $z=\\{7.1,9.1,12.2\\}$), FLARES (Lovell et al., 2021; Vijayan
et al., 2021; Wilkins et al., 2023a, $z=\\{7,9,11\\}$), DREaM (Drakos et al.,
2022, $z=\\{6.5,8.5,11\\}$), and SIMBA-EoR (Davé et al., 2019; Wu et al.,
2020, $z=9$). It is challenging to identify the reasons for the differences
between these simulations due to the differing initial assumptions made and
techniques used. For the most part, however, there is no large discrepancy
between the simulations at $6.5<z<8.5$, although they begin to diverge at
$11<z<13$. This is most likely due to the fact that these simulations are
often calibrated to observational results, more of which are available from
HST at $z\simeq 7$ than at $z\simeq 12$. When comparing to these studies, our
results show a far flatter $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ in our
lowest redshift bin at $6.5<z<8.5$ as well as much bluer $\beta$ slopes at
$11<z<13$ than the simulations (also seen in Cullen et al. (2023b)).

####  5.2.1 Flat $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ at $6.5<z<8.5$

The flatter observed $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ at
$6.5<z<8.5$ could be due to a number of factors. In order for this slope to be
flattened compared to other results, we could have either discovered a new
population of faint red galaxies and/or missed a large population of faint
blue galaxies compared to other studies. In our bias simulations in § 4, we
have found that our SED fitting templates favor the selection of blue galaxies
over red ones, making it unlikely that we have missed a large sample of blue
galaxies unless they have $\beta<-3.1$ (the bluest Larson et al. (2023)
template), which is challenging to reproduce with standard SPS models and
IMFs. Our sample covers less sky area than many of these other HST studies and
even the JWST work of Cullen et al. (2023b) which contains bright $z>7.5$
objects from COSMOS/UltraVISTA (see Donnan et al., 2023) meaning we almost
certainly miss the rarer bright red galaxies at $6.5<z<8.5$ in this work
compared to others. As well as this, we note that we could have either over-
corrected for the $\beta$ bias, which can reduce our
$\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ by $\sim 0.04$ towards the other
observational results or we could be impacted by sample contamination (this is
explored in more detail in § 5.4).

####  5.2.2 Blue $\beta(M_{\mathrm{UV}}=-19)$ at $11<z<13$

In our highest redshift bin, we observe a bluer $\beta$ than predicted by the
SC SAM GUREFT, DELPHI, FLARES, and DREaM simulations. Previous JWST
observations by Topping et al. (2023) and Cullen et al. (2023b) find
$\beta(M_{\mathrm{UV}}=-19)=\\{-2.42\pm 0.13,-2.63\pm 0.09\\}$ at $z\simeq 12$
and $11<z<12$ respectively, whereas we observe bluer values of
$\beta(M_{\mathrm{UV}}=-19)=-2.73\pm 0.06$. We note that this may well be
exacerbated by sample contamination, which we assess in § 5.4, although for
the rest of this subsection we explore the physical interpretation of this
ultra-blue average $\beta$ measurement.

We use the FLARES simulations to test the impact of changing
$f_{\mathrm{esc,LyC}}$ and $A_{\mathrm{UV}}$ of both
$\beta(M_{\mathrm{UV}}=-19)$ and $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$
in our highest redshift $11<z<13$ bin. The results from FLARES that are
plotted in Figure 11 assume $f_{\mathrm{esc,LyC}}=0$ with the dust attenuation
switched on. We find that turning this dust attenuation off both changes
$\beta(M_{\mathrm{UV}}=-19)$ by $\simeq-0.2$ to
$\beta(M_{\mathrm{UV}}=-19)\simeq-2.53$ and flattens the slope from
$\simeq-0.2$ to $\simeq 0.0$. The FLARES
$\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ predictions are approximate only
as their $\beta-M_{\mathrm{UV}}$ does not strictly follow a power law. Since
our measured $11<z<13$ $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}\simeq-0.2$
(which is fixed to the slope in the $8.5<z<11$ bin) is similar to the dusty
FLARES simulations, we conclude that turning off the dust law completely does
not adequately solve the problem. If we additionally switch off the galactic
nebular emission (i.e. by setting $f_{\mathrm{esc,LyC}}=1$), we obtain a
significantly bluer $\beta(M_{\mathrm{UV}}=-19)\simeq-2.78$ and shallower
$\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}\simeq-0.05$. We therefore conclude
that a combination of a reduction in dust attenuation and increase in Lyman
continuum escape fractions is expected towards higher redshifts. If the v2.2.1
BPASS SPS models (Stanway & Eldridge, 2018) and Chabrier (2003) IMF assumed in
FLARES are indeed correct, then a galaxy at $M_{\mathrm{UV}}=-19$ is expected
to be completely dust free with $f_{\mathrm{esc,LyC}}=1$ in order to match our
observations.

The FLARES multi-component dust model calculates the UV optical depth by
including contributions from both the ISM, $\tau_{\mathrm{ISM}}$ as well as
additional attenuation from young stars in birth clouds (BC) $<10$ Myr old,
$\tau_{\mathrm{BC}}$ (Charlot & Fall, 2000). These V-band optical depths are
calibrated by the dust-to-metal (DTM; see equation 15 from Vijayan et al.,
2019) and integrated line-of-sight metal column densities (for the ISM optical
depth), and spatially resolved stellar metallicities (for
$\tau_{\mathrm{BC}}$). In addition, the ISM normalizing factor
$\kappa_{\mathrm{ISM}}$ is chosen to match the $z=5$ Bouwens et al. (2015b) UV
LF, and the BC normalizing factor $\kappa_{\mathrm{BC}}$ is chosen to match
the $z=5$ $\beta$ observations from Bouwens et al. (2012, 2014) and the $z=8$
[OIII] $\lambda\lambda 4959,5007$+H$\beta$ EW distribution from De Barros et
al. (2019). Since we expect there to be some contribution from dust in order
to retain the $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ slope at these high-
redshifts, there may exist either a lag between dust and metal production
and/or a change in grain shape, size, and composition (impacting
$\kappa_{\mathrm{ISM}}/\kappa_{\mathrm{BC}}$). A more in depth discussion of
the dust implications of this work is presented in § 6.2.

Table 5: Bootstrapped binned results of median $M_{\mathrm{UV}}$ and bias corrected power law $\beta$ for our EPOCHS-III sample as plotted in Figure 11. We show results at $6.5<z<8.5$ in four $M_{\mathrm{UV}}$ bins, as well as in three and two bins at the higher redshifts $8.5<z<11$ and $11<z<13$ respectively. $M_{\mathrm{UV}}$ bin | $\langle M_{\mathrm{UV}}\rangle$ | $\langle\beta\rangle$
---|---|---
$6.5<z<8.5$
$M_{\mathrm{UV}}<-20.5$ | $-20.81^{+0.06}_{-0.07}$ | $-2.27^{+0.15}_{-0.15}$
$-20.5<M_{\mathrm{UV}}<-19.5$ | $-19.82^{+0.03}_{-0.04}$ | $-2.38^{+0.07}_{-0.07}$
$-19.5<M_{\mathrm{UV}}<-18.5$ | $-19.07^{+0.03}_{-0.03}$ | $-2.27^{+0.05}_{-0.06}$
$M_{\mathrm{UV}}>-18.5$ | $-18.07^{+0.05}_{-0.05}$ | $-2.21^{+0.09}_{-0.09}$
$8.5<z<11.0$
$M_{\mathrm{UV}}<-20.0$ | $-20.36^{+0.08}_{-0.10}$ | $-2.20^{+0.21}_{-0.20}$
$-20.0<M_{\mathrm{UV}}<-19.0$ | $-19.49^{+0.05}_{-0.05}$ | $-2.43^{+0.13}_{-0.16}$
$M_{\mathrm{UV}}>-19.0$ | $-18.59^{+0.15}_{-0.13}$ | $-2.59^{+0.28}_{-0.26}$
$11.0<z<13.0$
$M_{\mathrm{UV}}<-19.5$ | $-20.05^{+0.11}_{-0.10}$ | $-2.61^{+0.20}_{-0.21}$
$M_{\mathrm{UV}}>-19.5$ | $-19.18^{+0.16}_{-0.12}$ | $-2.55^{+0.24}_{-0.25}$
|  |
Figure 11: $\beta-M_{\mathrm{UV}}$ split by redshift for $6.5<z<8.5$ (upper),
$8.5<z<11.0$ (middle) and $11.0<z<13.0$ (lower). The shaded blue line shows
the $16^{\mathrm{th}}-84^{\mathrm{th}}$ percentile of the fit to the power law
$\beta$ bias-corrected EPOCHS-III sample (black background points), with the
slope at $11.0<z<13.0$ bin fixed to that at $8.5<z<11.0$. Comparisons to
observations (Finkelstein et al., 2012; Dunlop et al., 2012; Bouwens et al.,
2012; Dunlop et al., 2013; Bouwens et al., 2014; Wilkins et al., 2016;
Bhatawdekar & Conselice, 2021; Topping et al., 2023; Morales et al., 2023;
Cullen et al., 2023b; Roberts-Borsani et al., 2024) and simulations (THESAN
(Kannan et al., 2022, $z=\\{7,9\\}$), SC SAM GUREFT (Yung et al., 2023, 2024,
$z=\\{7,9,11\\}$), DELPHI (Mauerhofer & Dayal, 2023, $z=\\{7.1,9.1,12.2\\}$),
FLARES (Lovell et al., 2021; Vijayan et al., 2021; Wilkins et al., 2023a,
$z=\\{7,9,12\\}$), DREaM (Drakos et al., 2022, $z=\\{6.5,8.5,11\\}$), and
SIMBA-EoR (Davé et al., 2019; Wu et al., 2020, $z=9$)) are plotted in the left
and right panels respectively. The black lines show limits of the JAGUAR
catalogue which mimics the CEERS observations, and the red lines show the 20%
completeness contours. Figure 12: Amplitude (upper panel) and power law slope
(lower panel) of the $\beta-M_{\mathrm{UV}}$ relations in Figure 11 showing
the redshift evolution of $\beta(M_{\mathrm{UV}}=-19)$ and
$\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ respectively, as defined in
Equation 10. $\beta$ results for both the bias corrected/uncorrected power law
(pink/turquoise diamonds) and SED fitting (beige circles) methods are shown.
Observational results collated from the literature from HST (Bouwens et al.,
2012, 2014; Bhatawdekar & Conselice, 2021) and JWST (Cullen et al., 2023b;
Topping et al., 2023) are also plotted.

###  5.3 Comparisons with stellar mass

The relationship between the UV spectral slope and stellar mass of galaxies
has been studied in deep blank field surveys since the emergence of HST/WFC3IR
data at $1.0<\lambda_{\mathrm{obs}}/\mu\mathrm{m}<1.6$ in the early 2010s.
However, with the launch of JWST we now have access to rest-frame optical data
up to $z\sim 7-8$ allowing for better constraints on stellar masses at high
redshift, which presents an ideal opportunity to investigate this relation
further.

We plot our bias corrected power law $\beta$ values against our Bagpipes
derived $\log_{10}M_{\star}$ in Figure 13, fitting a power law function of the
form

$\beta=\beta_{0}+\frac{\mathrm{d}\beta}{\mathrm{d}\log_{10}(M_{\star}/\mathrm{M}_{\odot})}\log_{10}(M_{\star}/\mathrm{M}_{\odot})$
(11)

where $\beta_{0}=\beta(\log_{10}(M_{\star}/\mathrm{M}_{\odot})=0)$ is a
physically meaningless normalization factor. 20% completeness contours derived
from 5 realizations of the mock JAGUAR catalogue (see § 2.6 for details) are
shown as red lines, with the limits of the JAGUAR catalogue shown in black.
Observational (Finkelstein et al., 2012; Bhatawdekar & Conselice, 2021) and
simulated results from DELPHI (Mauerhofer & Dayal, 2023), SC SAM GUREFT (Yung
et al., 2023, 2024), DREaM (Drakos et al., 2022), and FLARES (Lovell et al.,
2021; Vijayan et al., 2021; Wilkins et al., 2023a) are shown in the left hand
panel. Bootstrapped median $\langle\beta\rangle$ and
$\langle\log_{10}M_{\star}\rangle$ values in stellar mass bins are given in
Table 6 and
$\mathrm{d}\beta/\mathrm{d}\log_{10}(M_{\star}/\mathrm{M}_{\odot})$ for bias
corrected/uncorrected power law $\beta_{\mathrm{PL}}$ and Bagpipes
$\beta_{\mathrm{SED}}$ are tabulated in Table 7 and plotted in comparison to
observational results from Finkelstein et al. (2012) and Bhatawdekar &
Conselice (2021) in Figure 14.

In our lowest redshift bin in Figure 13, we see that the bootstrapped median
points do not match the best-fit power law to the individual galaxies
particularly well. While the choice of bin size plays a role here, we note
that a power law is probably not the best fitting function to use in the
future, although we use it here to provide direct comparisons to previous HST
studies. In an attempt to quantify the threshold mass, $M_{\mathrm{thresh}}$,
that our bootstrapped median data points galaxies become “dust-free”, we fit a
piecewise power law function that is flat below $M_{\mathrm{thresh}}$ and a
power law at $M_{\star}>M_{\mathrm{thresh}}$ although without much success due
to the large intrinsic scatter in the data. We note also that our results at
$11<z<13$ suggest less dust attenuation on average than predicted by the
simulations, which was explored in more detail in § 5.2.2.

The two most important takeaways from Figure 13, however, are that the slope
at $6.5<z<8.5$ and $8.5<z<11$ is shallower than observed by both Finkelstein
et al. (2012) and Bhatawdekar & Conselice (2021) due to the presence of a
large sample of low mass ($M_{\star}\sim 10^{7.5}~{}\mathrm{M}_{\odot}$) red
($-2\lesssim\beta\lesssim-1$) galaxies, and that the slope at $11<z<13$ is
much steeper than at lower redshift due to the non-detection of these low mass
red objects.

####  5.3.1 Shallow $\mathrm{d}\beta/\mathrm{d}M_{\star}$ at $6.5<z<8.5$

While we do not go as deep as Bhatawdekar & Conselice (2021), who utilize the
strong gravitational lens of the MACS-0416 cluster to probe down to
$M_{\mathrm{UV}}\lesssim-13.5$ at $z=6$, we notice that our selection criteria
detect faint ($M_{\mathrm{UV}}\simeq-17.5$), low mass red objects more
efficiently than their bluer counterparts. The reasoning behind this is
somewhat unclear and depends on the difference in selection criteria used,
although we note that these redder $\beta$ could potentially become accessible
due to the increased rest-frame UV coverage up to $3000~{}\mathrm{\AA}$
provided by the F200W and F277W JWST/NIRCam wideband filters compared to
HST/WFC3IR. The extent to which we are able to select these low mass red
objects with JWST is apparent when comparing the 20% completeness curves
plotted in Figure 13 to those by Finkelstein et al. (2012) and Bhatawdekar &
Conselice (2021) (see their Fig. 7 and Fig. 6 respectively). At $z\simeq 7$,
these pre-JWST completeness curves do not cover detections of galaxies redder
than $\beta\simeq-2$ at $\log_{10}(M_{\star}/\mathrm{M}_{\odot})=7.5$, whereas
our JAGUAR completeness curves extend to $\beta\simeq-1.6$ at $6.5<z<8.5$.

####  5.3.2 Steepening of $\mathrm{d}\beta/\mathrm{d}M_{\star}$ with
increasing redshift

From Figure 14 we can see that our
$\mathrm{d}\beta/\mathrm{d}\log(M_{\star}/\mathrm{M}_{\odot})$ increases
significantly in our highest redshift bin from $0.31\pm 0.06$ at $8.5<z<11$ to
$0.81\pm 0.13$ at $11<z<13$ due to the non-detection of low mass red galaxies
that are seen in the lower redshift bins. At $11<z<13$, the most massive
galaxies at $M_{\star}=10^{9.5}~{}\mathrm{M}_{\odot}$ are very red (with the
reddest having $\beta=-1.3$), meaning that dust formation channels must
already exist at these redshifts. Due to the young ages of these galaxies,
Type II SNe are the most promising candidates as the major dust production
mechanism. Since such red galaxies are not observed at lower masses, it is
likely that the dust produced by these Type II SNe is lost in outflows instead
of being retained by the large gravitational potential wells of the most
massive galaxies (Finkelstein et al., 2012).

While our JAGUAR 20% completeness contours suggest that this parameter space
should be detectable, we note that there are few galaxies in our JAGUAR
simulation at $11<z<13$. This is complemented by the fact that the JAGUAR mock
SEDs may not adequately reproduce the real Universe at high redshift, meaning
that our completeness contours may not accurately represent the real
completeness limits in this redshift bin. Galaxies with $\beta\gtrsim-2$, such
as our sample of NIRCam selected red sources, lie outside of the HST 20%
completeness contours of Finkelstein et al. (2012) and Bhatawdekar & Conselice
(2021), meaning that these studies are likely completeness limited at these
redshifts. In our highest redshift bin, we therefore cannot conclude that the
trend is not induced by sample incompleteness since previous studies seem to
have been.

Figure 13: Redshift evolution of $\beta-M_{\star}$ for our EPOCHS-III sample, where the power law measured $\beta$ has been bias corrected. Bootstrapped median points are shown as beige circles, with individual galaxies in gray and the $16^{\mathrm{th}}-84^{\mathrm{th}}$ percentiles of the power law fit shown in blue. 20% completeness contours from JAGUAR are shown as red lines, with the JAGUAR limit shown in black. We compare to observations (Finkelstein et al., 2012; Bhatawdekar & Conselice, 2021; Morales et al., 2023; Roberts-Borsani et al., 2024) and simulations (SC SAM GUREFT (Yung et al., 2023, 2024, $z=\\{7,9,11\\}$), DELPHI (Mauerhofer & Dayal, 2023, $z=\\{7.1,9.1,12.2\\}$), FLARES (Lovell et al., 2021; Vijayan et al., 2021; Wilkins et al., 2023a, $z=\\{7,9,12\\}$), DREaM (Drakos et al., 2022, $z=\\{6.5,8.5,11\\}$) and SIMBA-EoR (Davé et al., 2019; Wu et al., 2020, $z=9$)) in the left and right panels respectively. The black lines show limits of the JAGUAR catalogue which mimics the CEERS observations, and the red lines show the 20% completeness contours. Table 6: Bootstrapped binned results of median $\log_{10}(M_{\star}/\mathrm{M}_{\odot})$ and bias corrected power law $\beta$ for our EPOCHS-III sample as plotted in Figure 13. We show results at $6.5<z<8.5$ in five $\log_{10}(M_{\star}/\mathrm{M}_{\odot})$ bins, as well as in three and two bins at the higher redshifts $8.5<z<11$ and $11<z<13$ respectively. The positive $\mathrm{d}\beta/\mathrm{d}M_{\star}$ slope is evident from these values in all three redshift bins. $\log_{10}(M_{\star}/\mathrm{M}_{\odot})$ bin | $\langle\log_{10}(M_{\star}/\mathrm{M}_{\odot})\rangle$ | $\langle\beta\rangle$
---|---|---
$6.5<z<8.5$
$\log_{10}(M_{\star}/\mathrm{M}_{\odot})<7.5$ | $7.26^{+0.02}_{-0.02}$ | $-2.48^{+0.07}_{-0.07}$
$7.5<\log_{10}(M_{\star}/\mathrm{M}_{\odot})<8.25$ | $7.82^{+0.03}_{-0.03}$ | $-2.40^{+0.06}_{-0.06}$
$8.25<\log_{10}(M_{\star}/\mathrm{M}_{\odot})<9.0$ | $8.66^{+0.03}_{-0.03}$ | $-2.31^{+0.07}_{-0.07}$
$9.0<\log_{10}(M_{\star}/\mathrm{M}_{\odot})<9.75$ | $9.25^{+0.05}_{-0.05}$ | $-2.10^{+0.11}_{-0.11}$
$\log_{10}(M_{\star}/\mathrm{M}_{\odot})>9.75$ | $10.04^{+0.28}_{-0.15}$ | $-1.38^{+0.27}_{-0.36}$
$8.5<z<11.0$
$\log_{10}(M_{\star}/\mathrm{M}_{\odot})<8.0$ | $7.57^{+0.05}_{-0.05}$ | $-2.59^{+0.16}_{-0.16}$
$8.0<\log_{10}(M_{\star}/\mathrm{M}_{\odot})<9.0$ | $8.47^{+0.11}_{-0.10}$ | $-2.54^{+0.18}_{-0.18}$
$\log_{10}(M_{\star}/\mathrm{M}_{\odot})>9.0$ | $9.28^{+0.15}_{-0.11}$ | $-2.02^{+0.23}_{-0.25}$
$11.0<z<13.0$
$\log_{10}(M_{\star}/\mathrm{M}_{\odot})<8.5$ | $8.00^{+0.13}_{-0.19}$ | $-2.81^{+0.18}_{-0.21}$
$\log_{10}(M_{\star}/\mathrm{M}_{\odot})>8.5$ | $8.87^{+0.09}_{-0.09}$ | $-2.55^{+0.22}_{-0.23}$
|  |
Figure 14: Evolution of the slope of the $\beta-M_{\star}$ relation, $\mathrm{d}\beta/\mathrm{d}\log(M_{\star}/\mathrm{M}_{\odot})$, for both our bias corrected (pink diamonds) and uncorrected (turquoise diamonds) power law $\beta$ measurements. We compare our results to those given in Finkelstein et al. (2012) (lime green) and Bhatawdekar & Conselice (2021) (orange), which use the SED fitting method to measure $\beta$. Consistently shallower results are measured at $z\simeq 7$ and $z\simeq 9$, and the first measurement at $z\simeq 11.5$ is shown. We caution that we may be completeness dominated in the highest redshift bin. Table 7: Median redshifts, amplitudes and slopes for our power law fits to the $\beta-M_{\mathrm{UV}}$ and $\beta-M_{\star}$ relations in our 3 redshift bins ranging $6.5<z<13$. We show results using both the power law ($\beta_{\mathrm{PL}}$) and SED fitting ($\beta_{\mathrm{SED}}$) methods to measure $\beta$. In the upper panel results for the bias corrected power law $\beta$ are shown in brackets, and have the same error as the uncorrected values. In our highest redshift bin ($11<z<13$) we fix the slope of our $\beta-M_{\mathrm{UV}}$ to that of our $8.5<z<11$ bin (shown by †). z bin | $N_{\mathrm{gals}}$ | $\langle z\rangle$ | $\beta(M_{\mathrm{UV}}=-19)$ | $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ | $\mathrm{d}\beta/\mathrm{d}\log(M_{\star}/\mathrm{M}_{\odot})$
---|---|---|---|---|---
$\beta_{\mathrm{PL}}$ ($\beta_{\mathrm{PL,corr}}$)
$6.5<z<8.5$ | $823$ | $6.98^{+1.01}_{-0.39}$ | $-2.28(-2.19)\pm 0.01$ | $-0.01(0.03)\pm 0.02$ | $0.22(0.22)\pm 0.02$
$8.5<z<11$ | $138$ | $8.97^{+1.56}_{-0.38}$ | $-2.59(-2.49)\pm 0.06$ | $-0.19(-0.17)\pm 0.06$ | $0.29(0.31)\pm 0.06$
$11<z<13$ | $50$ | $11.64^{+0.49}_{-0.43}$ | $-2.83(-2.73)\pm 0.06$ | $-0.19(-0.17)^{\dagger}$ | $0.70(0.81)\pm 0.13$
$\beta_{\mathrm{SED}}$
$6.5<z<8.5$ | $823$ | $6.98^{+1.01}_{-0.39}$ | $-2.40\pm 0.01$ | $-0.03\pm 0.02$ | $0.09\pm 0.01$
$8.5<z<11$ | $138$ | $8.97^{+1.56}_{-0.38}$ | $-2.48\pm 0.02$ | $-0.05\pm 0.03$ | $0.09\pm 0.04$
$11<z<13$ | $50$ | $11.64^{+0.49}_{-0.43}$ | $-2.53\pm 0.02$ | $-0.05^{\dagger}$ | $0.18\pm 0.07$
|  |  |  |  |

###  5.4 The impact of sample contamination

As in all galaxy samples, there will be some contaminant objects in our
sample; these most likely arise as lower redshift galaxy interlopers due to
Balmer-Lyman break degeneracy. In this work, we predict the number of these
interlopers we expect to find using the JAGUAR simulation (Williams et al.,
2018), as explained in more detail in § 2.6. We assign a contamination
percentage likelihood to each galaxy dependent on its origin survey and
position in the ($M_{\mathrm{UV}},\beta$) and ($M_{\star},\beta$) parameter
spaces. Adding these percentages together, we expect our sample to have
$\\{$40-90/823, 6-15/138, 2-3/50$\\}$ contaminants in the respective
$6.5<z<8.5$, $8.5<z<11$ and $11<z<13$ redshift bins (which are likely to be
the bluest in our sample) corresponding to 5-10%. This, of course, assumes
that the distribution of galaxy colours in JAGUAR matches the real Universe,
and may not accurately represent the real contamination in our sample. A
thorough comparison of the colours of simulations, including JAGUAR, with
NIRCam wide band photometric data from the CEERS survey is presented in
Wilkins et al. (2023b). Consequently, our simulation work identifies which
regions of parameter space may be subject to contamination, but the precise
contamination level may not be reliable. The reconciliation of this potential
systematic is deemed beyond the scope of this work.

We exclude objects from our relationship fits that lie in the regions of
parameter space found to have 10%, 25%, and 50% contamination likelihood
according to our JAGUAR based simulations and re-fit our $\beta-
M_{\mathrm{UV}}$ and $\beta-M_{\star}$ scaling relations to observe the
difference that this has on our results. We perform the power law fitting of
Equation 10 in the same fashion as in Figure 10, except this time we weight
each galaxy, $i$, by a factor
$w_{i}=1-\mathrm{Cont}_{i}(M_{\mathrm{UV}},\beta)$, where
$\mathrm{Cont}_{i}(M_{\mathrm{UV}},\beta)$ is calculated for each galaxy as
explained in § 2.6.

The impact of this on our $\beta-M_{\mathrm{UV}}$ fits is plotted in Figure
15. We find that the largest impact on the power law slope of our $\beta-
M_{\star}$ relation occurs at $6.5<z<8.5$, where we observe a decrease of
$\simeq 0.043$ from $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}=0.026\pm 0.017$
to $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}=-0.017\pm 0.017$ when removing
all galaxies with contamination likelihood $>10\%$. Negligible evolution in
$\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ is observed at $z>8.5$. We also
observe a significant reddening of $\beta(M_{\mathrm{UV}}=-19)$ by $\simeq
0.07$ ($6.5<z<8.5$), $\simeq 0.10$ ($8.5<z<11$), and $\simeq 0.16$ ($11<z<13$)
when removing contaminant galaxies with contamination likelihood $>10\%$. This
somewhat dampens our discussion on blue $\beta$ slopes should this high level
of contamination be accurate, although we note that even with this reddening
we still observe $\beta(M_{\mathrm{UV}}=-19)=-2.57\pm 0.06$ in our $11<z<13$
bin which is still significantly bluer than the most up-to-date simulated
results. Due to the small dependence on $\beta$ when binning the contamination
in $\Theta=(M_{\star},\beta)$, we find very little impact on the slope of
$\beta-M_{\star}$ at any redshift.

Figure 15: The impact of JAGUAR estimated contamination on the amplitude (at
$M_{\mathrm{UV}}=-19$; top panel) and slope (bottom panel) on our power law
$\beta$ bias corrected $\beta-M_{\mathrm{UV}}$ relation measured in our 3
redshift bins. Shown on the right are the results when no contaminant objects
are removed and moving leftwards we remove objects with
$\mathrm{Cont}(M_{\mathrm{UV}},\beta)<\\{0.5,0.25,0.1\\}$. Since we fix the
$11<z<13$ slope, the trend of $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$ with
contamination removal percentage in the lower panel matches the $8.5<z<11$
bin.

##  6 Discussion

###  6.1 An abundance of faint, low mass, red galaxies at $6.5<z<11$

In the two lowest redshift bins of Figure 13 we observe an abundance of
$\log_{10}(M_{\star}/\mathrm{M}_{\odot})<8.0$, $\beta>-2$ galaxies which are
missed by studies of Finkelstein et al. (2012) and Bhatawdekar & Conselice
(2021). These objects are also observed with $M_{\mathrm{UV}}\gtrsim-19$ at
$6.5<z<8.5$, which flatten the $\mathrm{d}\beta/\mathrm{d}M_{\mathrm{UV}}$.
These could be a result of our failure to accurately trace the UV continuum;
we may be biased red by high HI column density DLAs (already observed by JWST
(see Heintz et al., 2023a)) or strong resonant Mg II$-\lambda\lambda
2796,2803$ line emission (henceforth Mg II), which may help trace
$f_{\mathrm{esc,LyC}}$ (Chisholm et al., 2020; Katz et al., 2022).

In addition, Schaerer (2002) note that nebular dominated galaxies with hard
ionizing fields (large $\log U$ values) can be biased red should they exist in
the real Universe. Candidate nebular dominated galaxies, such as those
presented in Cameron et al. (2023b), would also provide evidence for a long
hypothesized top-heavy IMF in the early Universe should they exist. We use the
synthesizer111111https://github.com/flaresimulations/synthesizer (Vijayan et
al., 2021) python code to measure the expected $\beta$ slope of a pure nebular
dominated galaxy using a template with assumed $\log U=-2$, BPASS v2.2 SPS
model, and Chabrier (2003) IMF. We find $\beta_{\mathrm{neb,C94}}\simeq-1.0$
in the C94 filters (which have been seen to be biased red due to the increased
wavelength coverage approaching two-photon continuum emission turnover) and
$\beta_{\mathrm{neb,2-window}}\simeq-1.4$ in two rest wavelength windows
covering $\lambda_{\mathrm{rest}}=\\{1250-1750,2250-2750\\}~{}\mathrm{\AA}$.
This two-window $\beta$ more accurately reflects our power law measurement
technique, and therefore any nebular dominated galaxies, should they exist in
our sample, would appear with $\beta\simeq-1.4$.

We do not attempt to quantify the number density of DLAs and nebular dominated
galaxies in our EPOCHS-III sample as the defining features of these SEDs are
not traceable with the JWST/NIRCam wide-bands used in these early blank-field
surveys. Increased medium band coverage utilizing JWST/NIRCam F210M, F250M,
F300M, and F335M probing
$2000<\lambda_{\mathrm{rest}}~{}/~{}\mathrm{\AA}<3000$ at $6.5<z<8.5$ can be
used to both estimate Mg II EWs and detect the Balmer jump at
$3646~{}\mathrm{\AA}$ in nebular dominated galaxies. The detection of these
features in unbiased medium-band photometric surveys, and subsequent modelling
of the expected number densities of these systems, will allow us to
distinguish whether these red galaxies are indeed as dusty as their red
$\beta$ suggests.

###  6.2 Dust implications

From the results of our $\beta-M_{\star}$ relations, we propose possible
scenarios for the average build up of galactic dust. We start in our highest
redshift bin at $11<z<13$ where we observe a steep
$\mathrm{d}\beta/\mathrm{d}\log_{10}(M_{\star}/\mathrm{M}_{\odot})=0.81\pm
0.13$. The most massive galaxies in this redshift bin have $M_{\star}\simeq
10^{9.5}~{}\mathrm{M}_{\odot}$ and $\beta\simeq-1.5$, meaning galactic dust
must have been formed at these epochs. Due to the relatively young ages of
galaxies at these redshifts, we attribute this to dust formed in Type II SNe
on the smallest $\sim 10$ Myr timescales (see also Finkelstein et al. (2012)).
Since these SNe must also occur in the lowest mass systems, the lack of
observed red, $M_{\star}<10^{8}~{}\mathrm{M}_{\odot}$ galaxies means that this
dust is not retained, most likely via SNe feedback induced outflows which
remove dust from the small gravitational potential wells of the low mass host
galaxies. As well as this, low mass galaxies may well be prefentially bursty,
meaning radiation pressure from O/B-type main sequence stars may contribute
here.

Following the initial phase of Type II SNe dust production domination at
$z>11$, we observe an average reddening in galaxies at all stellar masses.
This reddening is more prominent for the lower mass galaxies which flattens
the slope of the $\beta-M_{\star}$ relation. One plausible explanation for
this reddening is by dust production on slightly longer timescales than Type
II SNe from short-lived carbon-rich Wolf-Rayet (WC) stars, which have been
seen to produce copious amounts of dust
($10^{-10}-10^{-6}~{}\mathrm{M}_{\odot}\mathrm{yr}^{-1}$ per WC) in the local
Universe (Lau et al., 2020). This production mechanism has been somewhat
overlooked in the past due to the requirement of an O/B main-sequence
companion star for efficient dust production, however with newer BPASS SPS
models this can now be studied more thoroughly (see Lau et al., 2020).

While WR stars are quite rare in typical local IMFs, more top-heavy IMFs which
may be naively expected in the early Universe (proposed by e.g. Rasmussen
Cueto et al., 2023; Steinhardt et al., 2023) may produce enough WC stars (if
sufficiently carbon-enriched) to fully account for the reddening observed
between $z\simeq 12$ and $z\simeq 9$ should dust destruction processes not be
prevalent in the reverse shock of the resulting SNe. Quantitatively, the IMF
required for the two-photon nebular continuum to dominate in GS-NDG-9422 from
Cameron et al. (2023b) would produce a WR star for every $\sim
140~{}\mathrm{M}_{\odot}$ of normal stellar population compared to every $\sim
1300~{}\mathrm{M}_{\odot}$ in a typical local IMF. Evidence for dust produced
via WC binaries, for instance from the
$\lambda_{\mathrm{rest}}=2175~{}\mathrm{\AA}$ carbonaceous signature found at
$z=6.71$ by Witstok et al. (2023a), would likely point towards a more top-
heavy IMF.

By $z=6.5$ the Universe is approximately $825$ Myr old, with $500$ Myr of time
having passed since the upper redshift limit of this study at $z=13$. Should
all galaxies in our $6.5<z<8.5$ bin be $\sim 500$ Myr old, we expect the
stellar winds of asymptotic giant branch (AGB) stars to dominate the
production of dust. Although these stars are less massive
($0.5-8~{}\mathrm{M}_{\odot}$), longer lived ($\gtrsim 100$ Myr timescales)
and produce less dust than WR stars, they are also more common. In addition,
ISM dust grain growth may also be prevalent in this stage of galaxy evolution,
although little is known regarding the size, shape, and chemical composition
of dust grains at $z=6.5$ beyond theoretical predictions by, for example,
Hensley & Draine (2023).

As well as the average reddening, we also observe an increasing scatter in
$\beta-M_{\star}$ with decreasing redshift. Even though at lower redshift,
there is more diversity in ages, metallicities, and environment at a given
mass, Bouwens et al. (2012) show (in their Fig. 18) that the scatter in dust
attenuation provides the largest impact on the scatter in $\beta$. We
conclude, therefore, that the increase in scatter between $z=13$ and $z=6.5$
is due to the variety of dust production methods available in these galaxies
(from Type II SNe, ISM dust grain growth, or WC/RSG/AGB stars), which are on
average older, and may lead to a diversity of dust attenuation laws diverging
from the standard Calzetti et al. (2000) or SMC/LMC laws adopted in the local
Universe.

###  6.3 A robust sample of blue $\beta<-2.8$ galaxies

Our fits to $\beta-M_{\mathrm{UV}}$ in Figure 12 show that on average the
galaxy population is uniformly blue at $z\gtrsim 11.0$, most likely due to the
lack of dust at these high redshifts. Much work has been done in the past to
investigate galaxies with ultra-blue $\beta$ slopes with minimal
$A_{\mathrm{UV}}$ dust extinction, including the possibility of Pop. III stars
and an increasingly top-heavy IMF (e.g. Schaerer, 2002, 2003; Raiter et al.,
2010b; Zackrisson et al., 2011).

We produce a tiered sample of 68 ultra-blue galaxies with
$\beta+\sigma_{\beta}<-2.8$, 35 of which have
$\mathrm{Cont}(M_{\mathrm{UV}},\beta)<0.2$ with 16 galaxies from NGDEEP and
JADES-Deep-GS having $\mathrm{Cont}(M_{\mathrm{UV}},\beta)<0.1$. We plot this
tiered sample in Figure 16 with transparency defined by
$1-\mathrm{Cont}(M_{\mathrm{UV}},\beta)$.

Figure 16: Sub-sample of 68 ultra-blue $\beta+\sigma_{\beta}<-2.8$ candidates
with 33 and 16 having less than 20% and 10% contamination likelihood plotted
as coloured points with more solid colours (less transparency). The rest of
the EPOCHS-III sample within our redshift range of interest is plotted in