# Unveiling the main sequence of galaxies at $z\geq 5$ with the James Webb
Space Telescope: predictions from simulations

Jordan C. J. D’Silva,1 Claudia D. P. Lagos,1,2,3 Luke J. M. Davies,1
Christopher C. Lovell,4 Aswin P. Vijayan3,5
1International Centre for Radio Astronomy Research (ICRAR), M468, University
of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia
2ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D).
3Cosmic Dawn Center (DAWN).
4Centre for Astrophysics Research, School of Physics, Astronomy & Mathematics,
University of Hertfordshire, Hatfield AL10 9AB, UK.
5DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs.
Lyngby, Denmark E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We use two independent, galaxy formation simulations, Flares, a cosmological
hydrodynamical simulation, and Shark, a semi-analytic model, to explore how
well the James Webb Space Telescope (JWST) will be able to uncover the
existence and parameters of the star-forming main sequence (SFS) at $z=5\to
10$, i.e. shape, scatter, normalisation. Using two independent simulations
allows us to isolate predictions (e.g., stellar mass, star formation rate,
SFR, luminosity functions) that are robust to or highly dependent on the
implementation of the physics of galaxy formation. Both simulations predict
that JWST can observe $\geq 70-90$% (for Shark and Flares respectively) of
galaxies up to ${z\sim 10}$ (down to stellar masses of $\approx 10^{8.3}\,\rm
M_{\odot}$ and SFRs of $\approx 10^{0.5}\,\rm M_{\odot}\,yr^{-1}$) in modest
integration times and given current proposed survey areas (e.g. the Web COSMOS
$0.6\,\rm deg^{2}$) to accurately constrain the parameters of the SFS.
Although both simulations predict qualitatively similar distributions of
stellar mass and SFR, there are important quantitative differences, such as
the abundance of massive, star-forming galaxies, with Flares predicting a
higher abundance than Shark; the early onset of quenching as a result of black
hole growth in Flares (at $z\approx 8$), not seen in Shark until much lower
redshifts; and the implementation of synthetic photometry, with Flares
predicting more JWST-detected galaxies ($\sim 90\%$) than Shark ($\sim 70\%$)
at $z=10$. JWST observations will distinguish between these models, leading to
a significant improvement upon our understanding of the formation of the very
first galaxies.

###### keywords:

cosmology:theory – galaxies:star formation – galaxies:high redshift –
infrared:JWST

††pubyear: 2022††pagerange: Unveiling the main sequence of galaxies at $z\geq
5$ with the James Webb Space Telescope: predictions from simulations–A

##  1 Introduction

A ubiquitous feature of galaxies is the tight coupling between star formation
rates (SFR) and stellar masses that has been observed to exist out to high
redshifts with very little scatter over the breadth of the $\mathrm{SFR-
M_{\star}}$ plane (Brinchmann et al., 2004; Salim et al., 2007; Noeske et al.,
2007b, a; Whitaker et al., 2012a); the so called star forming main sequence
(SFS). To first order $\mathrm{SFR\propto M_{\star}^{\alpha}}$ where
$\mathrm{\alpha}$ is a constant; though, differing opinions exist in the
literature as to whether the SFS is best described as a single component power
law model (e.g. Wuyts et al. 2011; Whitaker et al. 2012b; Speagle et al. 2014;
Pearson et al. 2018) or a two component power law model (e.g. Lee et al. 2015;
Leslie et al. 2020; Thorne et al. 2021). Two component functions model the low
and high stellar mass end of the SFS separately to account for distinct slopes
in these mass regimes, which accordingly trace different populations of
galaxies (Whitaker et al., 2012a). The use of a double power law function to
describe the SFS is motivated by a turn over in the SFS that has been shown to
occur for massive galaxies on the $\rm{SFR-M_{\star}}$ plane that experience a
downturn toward lower SFRs away from a linear trend (Abramson et al., 2014;
Cook et al., 2020). The physical mechanism behind the SFS and its low scatter
is believed to be that galaxies relax toward the locus of the SFS by the
internal and self-consistent regulation of star formation. The inflow rate of
gas that is needed for star formation is in equilibrium with a combination of
the rate at which stars are being formed and the rate at which gas is diverted
from the galaxy via outflows driven by feedback from stars, supernovae and
active galactic nuclei (e.g. Lilly et al. 2013; Tacchella et al. 2016). This
equilibrium needs to be quickly regained when galaxies are thrown out of it to
explain the tightness of the SFS. Indeed, this has been demonstrated to be the
case in simulations (e.g. Schaye et al. 2010; Lagos et al. 2011; Lagos et al.
2014; Matthee & Schaye 2019).

Stellar masses and SFRs have been previously determined out to $z\sim 8$
(Bouwens et al., 2015; González et al., 2011; Katsianis et al., 2017) with
instruments like the Advance Camera for Surveys on board the Hubble Space
Telescope (HST; Beckwith et al. 2006b; Koekemoer et al. 2007; Grogin et al.
2011). The trouble with these results is that they are generally short
wavelength surveys that can be complicated by the effect of dust attenuation
(Draine, 2003). The results of these surveys are sensitive to the methods that
are used to correct for dust attenuation. This is particularly problematic as
it is unclear whether dust obscured galaxy populations have been robustly
accounted for, which have shown to become dominant at $z\sim 1-2$ (Casey et
al., 2014). At higher redshifts, it is unclear whether dust-obscured galaxies
can make a significant contribution to the cosmic SFR density or not (e.g.
Casey et al. 2018.) Sub-millimetre (sub-mm) instruments, like the Atacama
Large Millimetre Array (ALMA), have been able to constrain the sub-mm emission
of hundreds of galaxies (Fudamoto et al., 2021) at $z>1$, and as sub-mm light
is the result of the dust re-radiating absorbed short wavelength photons,
these measurements can be used to constrain the SFRs of these galaxies. The
problem with these observations is two fold. First, ALMA has a small field-of-
view, and hence it could potentially be missing a number of intermediate-to-
bright sub-mm bright galaxies, if their frequency is rare, leading to
systematic uncertainties in the estimated cosmic SFR density of more than a
factor of $2$ at $z\gtrsim 4$ (Casey et al., 2018). Second, it is difficult to
determine stellar masses from sub-mm observations (Michałowski et al., 2014).
While studies have attempted to confirm the existence of the SFS at $z\gtrsim
5$, the stellar mass range is too small ($<1$ dex) to obtain meaningful
measurements of the slope and scatter of the SFS (Pearson et al., 2018; Leslie
et al., 2020; Thorne et al., 2021), and the UV rest-frame wavelengths used to
measure stellar masses can be extremely affected by dust attenuation. This is
about to dramatically change, thanks to the recently launched James Webb Space
Telescope (JWST; Gardner et al., 2006), which is observing in the infrared
($\mathrm{\lambda\sim 0.5\to 30.0\mu m}$) with exquisite resolution and
sensitivity, a wavelength range that is less affected by dust than those
traced by the HST, even at high redshift. The Near Infrared Camera (NIRCam) on
board the JWST is an order of magnitude more sensitive and covers a greater
area per pointing than the already existing infrared capabilities of the Wide
Field Camera 3 on board the HST (Gardner et al., 2006). A key goal of the JWST
is to provide a complete census of galaxies out to $z\sim 10$.

A key outcome of new observations is to confirm or challenge the current
understanding of astrophysical processes. It is thus important to produce
tailored predictions for what the JWST will be able to uncover as it continues
to observe the distant Universe. For this, we can turn our attention to
physically motivated galaxy formation models. Galaxy formation simulations
have a rich history (e.g., Somerville & Davé, 2015), and have proven to be a
vital tool for interpreting and predicting physical observations of the
Universe, especially those from large multi-wavelength galaxy surveys (Cole et
al., 2000; Baugh, 2006; Benson, 2010; Vogelsberger et al., 2014; Schaye et
al., 2015; Croton et al., 2016; Lagos et al., 2018). In order to make
meaningful predictions for the JWST, these simulations must include a
description of galaxy formation in combination with a description of the
spectral output of those simulated galaxies to directly test possible biases
affecting the observations.

In this work we use state-of-the-art simulations to address whether the JWST
will be able to predict the existence, and properties, of the SFS from $z=5\to
10$. Recent studies have explored JWST predictions (Vogelsberger et al., 2020;
Shen et al., 2021; Curtis-Lake et al., 2021; Wilkins et al., 2022a, b, c),
mostly focusing on luminosity functions (LF) and expected numbers of galaxies.
In addition, these generally present predictions within an individual model,
without exploring the dependence of their predictions on the details of galaxy
formation modelling. In this work instead, we make use of two independent
galaxy formation models, using drastically different techniques, the
cosmological hydrodynamical simulation suite Flares (Lovell et al., 2021;
Vijayan et al., 2021), and the semi-analytic model of galaxy formation Shark
(Lagos et al., 2018), with the aim of isolating predictions that appear robust
to the details of the models, and those that are highly model dependent.

Broadly, we wish to predict what kinds of galaxies the JWST will observe and
their distributions of stellar mass and SFRs of those galaxies from $z=5$ to
$z=10$. The main question is whether those observations would be enough to
establish the existence of a SFS out to $z=10$ and the parameters describing
it. We choose to make predictions on the galaxies that will be observed in the
$\mathrm{\sim 2\mu m}$, F200W filter with NIRCam on the JWST. As the F200W is
the most sensitive imaging filter it will be an important component in the
construction of high redshift sources (Rigby et al., 2022).

This paper is laid out as follows. In Section 2 we describe the main
differences between Flares and Shark, and we make note of the cosmological
parameters and initial mass functions (IMFs) that are used in the two
simulations. In Section 3 we calculate predictions of stellar mass functions
(SMF) and SFRs. In Section 4 we calculate predictions on the parameters of the
SFS and the stellar mass dependent scatter in the SFS. In Section 5 we predict
the cosmic stellar mass and star formation history. In Section 6 we discuss
our results, and in Section 7 we present our final conclusions and main
points. Unless explicitly specified, we have absorbed cosmological
dependencies in our results. Furthermore, all magnitudes are quoted in AB,
unless stated otherwise.

##  2 Simulations

We introduce the two galaxy formation simulations used in this work: Flares, a
hydrodynamical zoom simulation (Section 2.1), and Shark, a semi-analytic model
of galaxy formation (Section 2.2). We explore differences between the dust and
stellar population synthesis (SPS) models between the two simulations in
Section 2.3.

###  2.1 Flares

Flares (First Light And Epoch of Reionisation) was introduced in Lovell et al.
(2021) and Vijayan et al. (2021). Below we summarise the baryon model ($\S$
2.1.1) and the way the spectral energy distribution (SED) of galaxies is
computed in Flares ($\S$ 2.1.2).

####  2.1.1 Modelling and resimulation method

Flares is a series of cosmological hydrodynamic zoom simulations built on the
EAGLE galaxy formation model (Schaye et al., 2015; Crain et al., 2015) and
therefore utilises the modelling techniques and sub-grid recipes used in
EAGLE. EAGLE is a suite of smoothed-particle-hydrodynamics simulations run on
the P-GADGET3 N-body Tree-PM code that was last described in Springel et al.
(2005) and uses the ANARCHY (see Schaller et al. (2015) for details) code to
solve the coupled equations of hydrodynamics. EAGLE includes a swathe of sub-
grid recipes including radiative cooling and photo-heating, star formation,
stellar evolution and chemical enrichment, black hole growth and feedback from
active galactic nuclei (AGN) and massive stars (see Schaye et al. (2015);
Crain et al. (2015) for further details on the sub-grid models). EAGLE was
tuned to the $z\sim 0$ galaxy SMF, stellar mass-black hole relation and galaxy
sizes, but has also been shown to agree well with observations not explicitly
used in the tuning of the sub-grid model free parameters. Relevant for this
work, Furlong et al. (2015) showed that the SMF and SFS were reasonably well
reproduced in EAGLE up to $z=4$ and $z=3$, respectively; Katsianis et al.
(2017) showed that the SFR function was reasonably well reproduced up to
$z=4$; and both Katsianis et al. (2019) and Davies et al. (2019) showed that
EAGLE reproduced well the SFS’s scatter at $z=0$.

Unfortunately, the dynamic range of stellar masses and star formation rates in
the EAGLE suit is insufficient to probe the most massive and star forming
galaxies that exist in the most overdense environments. Being a series of zoom
simulations, Flares is able to sample these rare overdense environments from a
larger dark matter-only simulation thereby providing a statistically complete
set of environments. Flares concentrates on redshift snapshots
$z=[5,6,7,8,9,10]$ and is therefore perfectly suited to provide predictions on
what future surveys will be able to conclude about high-redshift galaxies and
the epoch of reionisation. Previous studies have sampled the rarest cluster
environments. Barnes et al. (2017); Bahé et al. (2017a) for example sampled
and resimulated 30 massive galaxy clusters, also with the EAGLE model, with
$z=0$ halo masses $\mathrm{10^{14}\leq M_{200}/M_{\odot}\leq 10^{15.4}}$.
Attempting to use such a sample to parametrise the universal properties of
high redshift galaxies is however troublesome because of the bias toward only
massive cluster galaxies that is not representative of the diverse range of
galaxy environments present in the Universe. A prediction of the properties of
all galaxies in the early Universe thus requires an unbiased set of simulated
galaxies existing in a range of environments and possessing a range of masses.
Flares has generated a set of galaxies that are representative of the entire
Universe by sampling from a range of 40 overdensities and applying statistical
weights to each region (Lovell et al., 2021). The zoom regions are selected
from the $\mathrm{3.2^{3}Gpc^{3}}$ dark matter-only EAGLE box, which is also
used in Barnes et al. (2017). 40 spherical zoom regions with radius
$\mathrm{14\mathit{h}^{-1}cMpc}$ are resimulated down to $z=4.67$ at the same
resolution as the fiducial EAGLE simulation and with identical physics and
parameters. The benefit of this is a sample representative of the early
universe with sufficient number statistics on the most massive, clustered
galaxies, but unbiased to those.

The range of overdensities of the large dark-matter only box is partitioned
into 50 bins of equal width in $\log_{10}(1+\delta)$, and the $i^{\rm th}$ bin
is assigned a statistical weight, $w_{\rm true,i}$, according to the
proportion of resimulated overdensities that exist in the bin against the
total quantity of overdensities in the box, such that $\sum_{\rm i}{w_{\rm
true,i}}=1$. Each of the resimulated regions are also distributed over the bin
size, and are assigned weights, $w_{\rm i,j}$, according to the proportion of
resimulated overdensities in the bin against the total quantity of the 40
resimulated overdensities, once again so that $\sum_{i}{w_{\rm ij}}=1$. The
weight per bin associated with the resimulated regions is then $w_{\rm
sample,i}=\sum_{\rm j}{w_{\rm ij}}$. Each density bin is then weighted by
$r_{\rm i}=w_{\rm true,i}/w_{\rm sample,i}$. In order for each resimulation
region, $\mathrm{j}$, to be statistically representative of the cosmic
distribution, it is weighted by $f_{j}=\sum_{i}r_{i}w_{i,j}$ with
$\sum_{j}f_{j}=1$. The weights do not change with redshift since the relative
ordering of the overdensities is expected to be invariant through $z=10\to 5$
as even the greatest overdensities would be a result of only mildly non-linear
evolution.

Structures are found in Flares first from a Friends-Of-Friends finder (Davis
et al., 1985), and bound substructure is further identified with the SubFind
algorithm (Springel et al., 2001).

The dark matter particle mass resolution of the large
$\mathrm{3.2^{3}Gpc^{3}}$ dark matter-only box is $\mathrm{m_{DM}=8.01\times
10^{10}M_{\odot}}$ and the gravitational softening length is
$\mathrm{59\,ckpc}$ (note this is the same simulation used in Barnes et al.
2017; Bahé et al. 2017b to select the regions to resimulate). The gas particle
mass and dark matter particle mass in the fiducial EAGLE simulation, and
therefore Flares, is $\mathrm{m_{gas}=1.8\times 10^{6}M_{\odot}}$ and
$\mathrm{m_{DM}=9.7\times 10^{6}M_{\odot}}$ respectively; and the softening
length is $\mathrm{\epsilon=2.66\,ckpc}$. EAGLE, and by extension Flares, use
the Planck Collaboration et al. (2014) $\Lambda$ cold dark matter
($\Lambda$CDM) cosmology with $\mathrm{\mathit{H_{0}}=\mathit{h}\times
100Mpc\,kms^{-1}}$, $\mathrm{\mathit{h}=0.6777}$; $\mathrm{\Omega_{m}=0.307}$
and $\mathrm{\Omega_{\Lambda}=0.693}$. Flares adopts a universal Chabrier
initial mass function (Chabrier, 2003).

####  2.1.2 Dust attenuation calculation in Flares

The light emitted by stars is attenuated by a two-phase dust medium, the birth
clouds, which attenuate young stars, and the diffuse dust in the interstellar
medium, which attenuates both young and older stars. Dust attenuation can be
quantified by an examination of the optical depth in the V-band (550nm),
$\mathrm{\tau_{V}}$, and seeing how this is affected by the dust. The emission
from each simulated galaxy is computed using version 2.2.1 of the Binary
Population and Spectral Synthesis (BPASS) stellar population synthesis (SPS)
code (Eldridge et al., 2017; Stanway & Eldridge, 2018). Stellar clusters form
below sub-kpc scales meaning that dust attenuation from the birth clouds must
be calculated using a sub-grid recipe, as:

$\tau_{\rm BC,V}=\begin{cases}\kappa_{\rm
BC}\times(Z_{\star}/0.01)&\mathrm{t\leq 10^{7}yr},\\\
0&\mathrm{t>10^{7}yr},\end{cases}$ (1)

where $\kappa_{BC}$ is a normalisation that encodes information about the
physical properties of the dust grains and the dust-to-metal ratio of the
birth clouds, and $Z_{\star}$ is the metallicity of the star particle. The
optical depth is sensitive to the age of the stellar populations and the
piecewise function is obtained from Charlot & Fall (2000) who find that birth
clouds disperse by $\mathrm{10^{7}}$ years, beyond which there are few
attenuating particles around the stellar cluster and hence the optical depth
of the birth cloud tends to zero.

In the case of the diffuse dust component, the metal content of the ISM is
used as a proxy of dust content (Vijayan et al., 2021). The metal content of
the interstellar medium is directly determined from the SPH particles.

The equation to calculate the V-band optical depth is

$\tau_{\rm ISM,V}=\mathrm{DTM}\,\kappa_{\mathrm{ISM}}\,\Sigma(x,y)$ (2)

where $\mathrm{DTM}$ is the dust-to-metal mass ratio of the galaxy,
$\kappa_{\mathrm{ISM}}$ is a normalisation that encodes information about the
physical properties of the dust grains and $\Sigma(x,y)$ is the metal column
density integrated over a line of sight. The $\mathrm{DTM}$ ratio adopted in
Flares is a fitting function Vijayan et al. (2019, Equation 15) that connects
the dust-to-metal mass ratio to the age and metallicity of the galaxy. This
function fits the relation between the three latter quantities in the semi-
analytic model of galaxy formation L-galaxies, and specifically the version
introduced in Vijayan et al. (2019) which includes a model for dust formation,
growth and destruction. $\Sigma(x,y)$ is evaluated by integrating the density
field of SPH particles and linking this to the metallicity and mass of each
particle along the line of sight. The viewing angle is fixed such that the
line of sight is along the z-axis of the simulation. This prescription yields
a sufficient calculation of interstellar medium attenuation without having to
simulate the distinct geometry and properties of the dust.

Attenuation from both physical components, described with Equation 1 and
Equation 2, is combined with a general wavelength dependence as

$\tau_{\rm\lambda}=(\tau_{\mathrm{ISM}}+\tau_{\mathrm{BC}})\times(\lambda/550\mathrm{nm)}^{-1}.$
(3)

The form of Equation 3 is reminiscent of the Charlot & Fall (2000) attenuation
curve. This is the attenuation curve for each star particle along the line of
sight. Lower metallicities at high redshift motivates the use of an
attenuation curve that is similar to the Small Magellanic Cloud; Vijayan et
al. (2021) remark that the slope of the attenuation curve in Flares is flatter
in the UV region when compared to the same curve of the Small Magellanic
Cloud, but not as flat as the attenuation curve of Calzetti et al. (2000).
Furthermore, the clumpiness of the ISM can affect the resultant attenuation
curve. The dust model is calibrated using the UV light function, UV continuum
slope and the [OIII]+H$\beta$ EW distribution (for further details see the
appendix of Vijayan et al., 2022). At all redshifts, $\kappa_{\mathrm{BC}}=1$
and $\kappa_{\mathrm{ISM}}=0.0795$. Flares does not include a redshift
evolution of the dust grain sizes, masses or composition.

###  2.2 Shark

The second simulation used in this work is the semi-analytic model of galaxy
formation Shark, first introduced in Lagos et al. (2018).

####  2.2.1 N-body skeleton and semi-analytic method

The cosmic structure in Shark is provided by the N-body dark matter only
simulation Synthetic UniveRses For Surveys: SURFS (Elahi et al., 2018). SURFS
was designed to provide theoretical test beds to compliment ongoing and
upcoming galaxy surveys. SURFS was run using a memory-lean version of the
GADGET-2 code (Springel, 2005), producing a variety of box sizes between
$\mathrm{40\mathit{h}^{-1}cMpc}$ and $\mathrm{210\mathit{h}^{-1}cMpc}$, each
including between $\mathrm{512^{3}}$ to $\mathrm{1563^{3}}$ dark matter
particles. Different simulation boxes also cover a variety of particle masses
and softening lengths. The SURFS box that underpins the Shark simulation used
in this work is the L210N1563 simulation that has box size
$\mathrm{L_{box}=210\mathit{h}^{-1}cMpc}$, number of dark matter particles
$\mathrm{N_{p}=1563^{3}}$, dark matter particle mass
$\mathrm{m_{DM}=2.21\times 10^{8}\mathit{h}^{-1}M_{\odot}}$ and softening
length $\mathrm{\epsilon=4.5\mathit{h}^{-1}ckpc}$. In total, 200 snapshots
spaced in logarithmic intervals of growth factor from redshift $z=24$ to $z=0$
were produced using the Planck Collaboration et al. (2016) $\Lambda$CDM
cosmology with cosmological parameters:
$\mathrm{\mathit{H_{0}}=\mathit{h}\times 100Mpc\,kms^{-1}}$,
$\mathrm{\mathit{h}=0.6751}$; $\mathrm{\Omega_{m}=0.3121}$,
$\mathrm{\Omega_{b}=0.0491}$ and $\mathrm{\Omega_{\Lambda}=0.6879}$ being the
Hubble’s constant, matter density, baryon density and $\mathrm{\Lambda}$
density respectively.

Dark matter halos and their substructure are identified using the halo finding
code, VELOCIraptor (Cañas et al., 2019; Elahi et al., 2019a), and merger trees
of dark matter halos are constructed using the halo merger tree builder,
TreeFrog (Elahi et al., 2019b).

The dark matter halo catalogues from SURFS provide the static skeletons that
Shark uses to simulate galaxies. Since links between halo descendants are
found for only up to 4 snapshots into the future, it can be the case that
discontinuities exist in the merger trees across snapshots. To smooth over
these discontinuities and enforce continuity over the equations of galaxy
formation, Shark places subhalos between the snapshots of the current subhalo
and its descendent, if snapshots were skipped by TreeFrog. The dark matter
reservoir in which galaxies are embedded are assumed to have an NFW profile
(Navarro et al., 1997).

Galaxies are evolved by first finding the halos that first appear in the
snapshot and do not have a progenitor. The most massive subhalo of these first
generation halos are assigned a gas supply of mass

$\mathrm{M_{gas}=\Omega_{b}/\Omega_{m}\times M_{halo}},$ (4)

and these artificial galaxies are then evolved forward in time. The
differential equations that determine the mean mass and metallicity of the
baryonic components that are simulated in Shark are given by Equations 49-58
in Lagos et al. (2018). In Shark, the physical modelling includes: baryonic
matter accretion onto halos, radiative gas cooling, star formation in discs
and bulges, stellar feedback, reincorporation of ejected gas, chemical
enrichment, galaxy mergers, disc instabilities, photoionisation feedback,
black hole growth and AGN feedback and environmental effects. As Shark lacks
the ability to model complex morphologies of galaxies, the basic morphological
description that is used is a bulge-disc distinction. This distinction is
important because the assembly of stellar mass, and subsequent chemical
enrichment, in either the disc or the bulge is not necessarily identical. This
distinction is reflected in the star formation law that is a function of the
molecular-to-atomic gas fraction, surface density of gas and pressure in
either the disc or the bulge, depending on where the star formation is
occurring. In addition, it is assumed that star formation proceeding in the
bulge in the form of starbursts converts the molecular gas into stars more
efficiently (by a factor of 10) than star formation occurring in the disc.

For extended details on this modelling, the reader is referred to Lagos et al.
(2018). Shark adopts a universal Chabrier IMF (Chabrier, 2003).

Figure 1: Obscured rest-frame FUV LFs for Flares (orange), Shark-default
(blue) and Shark-steep (light blue) at redshifts $z=6,8,10$, as labelled. The
dotted lines show the intrinsic LFs for Flares and Shark. Grey symbols with
error bars show observational constraints of the obscured rest-frame FUV LFs
from Finkelstein et al. (2015); Oesch et al. (2018); Bouwens et al. (2021). It
is apparent that Flares and Shark-steep offer the best fits to the
observations, while Shark-default predicts number densities that are too low.

####  2.2.2 Dust attenuation calculation in Shark

The work horse behind the lighting of Shark galaxies is ProSpect (Robotham et
al., 2020), a multiwavelength spectral energy distribution package capable of
generating and fitting spectral energy distributions under the consideration
of many different astrophysical mechanisms. ProSpect by default utilises the
SPS libraries of Bruzual & Charlot (2003), though the highly flexible nature
of ProSpect allows for the use of many different SPS libraries (such as
e-MILES, Vazdekis et al. 2016, and BPASS, Eldridge et al. 2017). Dust
attenuation is described using the model of Charlot & Fall (2000). Star
formation and metallicity histories from Shark galaxies are fed to ProSpect
with sensible assumptions about how dust attenuates galaxy emission and how
energy is then reradiated in the infrared. Much like Flares, the attenuation
is multicomponent in nature being related to the optical depth of both birth
clouds and the interstellar medium.

There are four distinct methods of dust attenuation that are included in
Shark, and the reader is referred to Lagos et al. (2019) for further details
on all of them. Here we only concentrate on methods 3 and 4 as referred to in
Lagos et al. (2019), and we briefly outline them below.

The optical depth of the diffuse ISM is modelled as

$\tau_{\mathrm{ISM}}=\tau_{\mathrm{ISM,V}}\times(\lambda/550\mathrm{nm})^{\eta_{\mathrm{ISM}}},$
(5)

Trayford et al. (2020) use the radiative transfer code SKIRT (Camps & Baes,
2015) to derive scaling relations between $\tau_{\mathrm{ISM,V}}$,
$\eta_{\mathrm{ISM}}$ and the surface density of dust $\Sigma_{\mathrm{dust}}$
for galaxies between $z=0\to 2$ in the EAGLE simulation, finding the
relationships to be mostly independent of redshift. In Shark the dust surface
density is computed for each galactic component (disc and bulge), and the
attenuation parameters of Equation 5 are determined by sampling the scaling
relations of Trayford et al. (2020). Note that these parameters are found
independently for the bulge and disc.

The opacity of light from the birth clouds is calculated with

$\tau_{\rm
BC}=\tau_{\mathrm{ISM}}+{\tau}_{\mathrm{BC,V}}\times(\lambda/550\mathrm{nm})^{\eta_{\mathrm{BC}}}.$
(6)

The parameter $\tau_{\mathrm{ISM}}$ appears in Equation 6 because light from
stars inside the birth clouds is attenuated by both the birth clouds
themselves and the ISM. The V-band optical depth of the birth cloud is derived
from the following equation:

$\tau_{\mathrm{BC,V}}=\tau_{\mathrm{BC,0}}\frac{\mathrm{DTM}\,Z_{\mathrm{gas}}\,\Sigma_{\mathrm{gas,cl}}}{\mathrm{DTM_{MW}}\,Z_{\odot}\,\Sigma_{\mathrm{MW,cl}}}$
(7)

$\tau_{\mathrm{BC,0}}=1$, DTM is the dust-to-metal ratio, $Z_{\mathrm{gas}}$
is the gas metallicity, $\Sigma_{\mathrm{gas,cl}}$ is the surface density of
gas, $\mathrm{DTM_{MW}}=0.33$ is the dust-to-metal ratio of the Milky Way,
$Z_{\odot}=0.0189$ is the solar metallicity and
$\Sigma_{\mathrm{MW,cl}}=85\rm\,M_{\odot}\mathrm{pc}^{-2}$ is the typical
surface density of molecular clouds in the Milky Way (Krumholz et al., 2009).
The value of $\eta_{\mathrm{BC}}$ is taken to be the default for the Charlot &
Fall (2000) attenuation model: $\eta_{\mathrm{BC}}=-0.7$.

The dust masses, which are used to compute $\Sigma_{\mathrm{dust}}$, are by
default determined from the metallicity of the gas in the galaxy according to
the best fitting $\mathrm{M_{dust}/M_{Z}-Z_{gas}}$ relation computed by Rémy-
Ruyer et al. (2014). An alternative method for calculating dust masses uses a
steeper fit to the $\mathrm{M_{dust}/M_{Z}-Z_{gas}}$ within the errors of the
best-fitting relation, which is more consistent with the more recent data of
De Vis et al. (2019). Shark uses these as the dust-to-metal ratio. The surface
density of the dust is computed for discs and bulges separately according to
the following equations;

$\mathrm{\Sigma_{dust,disc}=\frac{0.5\,M_{dust,disc}}{\pi
r_{50,disc},\,l_{50}}},$ (8)

where $\mathrm{M_{dust,disc}}$ is the dust mass in the disc,
$\mathrm{r_{50,disc}}$ is the half-gas mass radius of the disc and
$\mathrm{l_{50}=sin(i)\times(r_{50,disc}-r_{50,disc}/7.3)+r_{50,disc}/7.3)}$
is the projected minor axis with inclination $\mathrm{i}$. The factor $7.3$
originates from the average ratio between the scale height and scale length
observed in local galaxy discs (Kregel et al., 2002), and the inclination is
determined from the angular momentum vector of the host sub halo or chosen
randomly for orphan galaxies. For bulges,

$\mathrm{\Sigma_{dust,bulge}=\frac{0.5\,M_{dust,bulge}}{\pi
r_{50,bulge}^{2}}},$ (9)

where $\mathrm{M_{dust,bulge}}$ is the dust mass of the bulge and
$\mathrm{r_{50,bulge}}$ is the half-gas mass radius of the bulge. Bulges are
considered to be spherically symmetric so the inclination is unimportant.

The two models are thus referred to as

  * •

EAGLE-$\mathrm{\tau}$ RR14, which is the default for Shark and uses the EAGLE
parametrisation of the Charlot & Fall (2000) model presented in Trayford et
al. (2020) and the best-fitting relation of the
$\mathrm{M_{dust}/M_{Z}-Z_{gas}}$ from Rémy-Ruyer et al. (2014).

  * •

EAGLE-$\mathrm{\tau}$ RR14-steep that is the same as above but uses a steeper
relation for the $\mathrm{M_{dust}/M_{Z}-Z_{gas}}$. This indicates that this
model will have lower dust masses for fixed metallicities than the previous
model.

We choose to concentrate on only two of the possible four dust models included
in Shark here for the purpose of demonstrating the varying effects of dust
models upon a single simulation.

###  2.3 SPS and Dust Models: comparing Flares and Shark

We investigate how different parameters used in the generation of the
artificial spectral energy distributions affect the emission of our simulated
galaxies. Later we use synthetic photometry to apply magnitude cuts in
galaxies to isolate those that would be considered JWST detected.

Figure 1 shows the obscured rest-frame far ultraviolet (FUV, $\mathrm{\sim
1500}$Å) LFs at $z=6,8,10$. We show the effect of using both the default and
steep dust models in Shark as described in Section 2. We do not have the
flexibility to adopt different methods for determining the parameters in
Equations 1 and 2 for Flares in an efficient way as the use of a range of
possible recipes that treat the dust is computationally expensive to run. We
also do not perform any stellar mass cut in these plots. Observational
constraints of the rest-frame UV LFs at these redshifts from Finkelstein et
al. (2015); Oesch et al. (2018); Bouwens et al. (2021) are shown. These use a
combination of optical and near infrared images from the Hubble Space
Telescope. We also show the unobscured LFs of each simulation for reference.
The unobscured LFs are similar to the obscured ones at the faint end but
diverge at the bright end, where there is significant difference between the
intrinsic and obscured brightness per unit volume, showing that the effect of
dust upon the FUV LF is significant, even at $z=10$ in both Flares and Shark.
The kinks in the intrinsic LFs of Flares for $\mathrm{M_{UV}\lesssim-21}$ are
due to the overdense regions that Flares samples at fixed magnitude, which
Shark does not due to its limited volume. Because of the highly flexible
nature of Shark and ProSpect, we are able to investigate the effect of using
the BPASS (Eldridge et al., 2017; Stanway & Eldridge, 2018, version 2.2.1) and
BC03 (Bruzual & Charlot, 2003) SPS libraries on the intrinsic rest-frame UV
and optical emission, and we found that the differences between using these
two libraries is minimal. Therefore, although Flares uses BPASS, we continue
to use the default BC03 library with Shark for consistency with previous Shark
papers.

The default Shark dust model causes a greater deficit in the FUV brightness
than the steep model over nearly all magnitudes at each redshift, as expected
from the higher dust mass at fixed metallicity in the this model than in the
steep model. While both simulations agree reasonably well with the
observational results, such as the UV-FIR emission of galaxies (e.g., Lagos et
al., 2018; Lagos et al., 2019; Vijayan et al., 2021, 2022), it is important to
highlight the large uncertainties that characterise current observational
measurements. Due to the Shark-steep dust model matching these observations
far better than the default model, we will use the steep dust model wherever
we make use of Shark photometry within the subsequent sections. Lagos et al.
(2019) showed that the normalisation of the $z=0$ FUV LF predicted by steep
Shark model was slightly higher around the break, $\mathrm{L^{*}}$, than the
default model by a fraction of a dex. The fact that the steep model provides a
better fit to the observed UV LFs at $z>6$ means that either we require a
redshift dependent dust-to-metal ratio relation or a chemical evolution model
that enriches galaxies more slowly than the current instantaneous recycling
model in Shark. We come back to this discussion later in the context of the
predicted cosmic stellar mass and SFR history.

Figure 2: The JWST NIRCAM F200W apparent magnitude volume weighted
distribution in Flares (orange), Shark-default (blue) and Shark-steep (light
blue) at redshifts $z=6,8,10$. We have performed a stellar mass cut of
$\mathrm{M_{\star}\geq 10^{8.5}M_{\odot}}$. The faint dotted lines show the
LFs for all galaxies in Flares and Shark-steep.

The distribution of galaxies that will be recoverable with the JWST, using
NIRCam and the F200W filter for example, depends on the fraction of
$\mathrm{\lambda\lesssim 0.5\mu m}$ photons that escape from the galaxy at
$z\gtrsim 5$. Figure 2 shows the apparent magnitude distribution per unit
volume in the F200W band for Flares and Shark galaxies. As with the FUV LFs,
both simulations predict similar shapes to the F200W distributions over the
breadth of the magnitude domain and throughout redshift. The Shark steep model
predicts more bright galaxies than the default model, and the cause of this
difference is the same as it is for the UV LFs in Figure 1. These differences
are, however, lesser than those seen in Figure 1 as the rest wavelength traced
by the JWST is less affected by dust than that traced by the Hubble Space
Telescope, for example, at $z\gtrsim 5$ (we elaborate on this further in
Section 3). Also note that both simulations have a peak at around $\sim 30$
mag. This peak is driven by resolution effects becoming significant at higher
magnitudes, as galaxies with stellar masses $\lesssim 10^{8.5}\,\rm M_{\odot}$
start to dominate the number density fainter than $30$ mag as can be seen from
the dotted lines that show the distribution for these galaxies in the
simulations. The sensitivity of the JWST in the F200W filter is predicted to
be 29 mag, as per the exposure calculator tool 111The 29mag cut comes from
using an integration time of $10,000$ seconds and imposing detections to be
above $5\sigma$. The exposure calculator for JWST can be found here
https://jwst.etc.stsci.edu. Hence, throughout the text we will refer to JWST
detected galaxies as those brighter than 29 mag in F200W.

The FUV and Nircam F200W LFs highlight clear differences between Flares and
Shark. Although variations in the dust model used in simulations can alleviate
this tension, fundamental differences in the way that either simulation
assembles stellar mass affect the predicted emission of galaxies (e.g.,
Bellstedt et al., 2020; Koushan et al., 2021; Wilkins et al., 2022b). These
processes are interconnected as dust tracks metal enrichment that is itself
connected to star formation and stellar evolution. As such, it is useful to
investigate how each simulation models the distribution of stellar mass and
SFR in galaxies out to redshift $z=10$ and whether the JWST will be able
improve our understanding of them.

##  3 Stellar mass functions and star formation rate functions

Figure 3: SMFs (top) and SFRFs (bottom) at redshifts $z=6,8,10$. The light
coloured, dotted lines show these quantities for the total population of
galaxies in the Flares (orange) and Shark (blue). The light coloured dashed
lines show these quantities for the galaxies with stellar masses above
$\mathrm{10^{8.5}M_{\odot}}$ in Flares (orange) and Shark (blue). The solid
lines show these quantities for the JWST detected populations in Flares (red)
and Shark (cyan). Points with error bars show observational constraints of the
SMF from González et al. (2011); Song et al. (2016); Stefanon et al. (2021)
and SFRF from Katsianis et al. (2017). Figure 4: Fraction of galaxies that are
brighter than 29 magnitudes in a population selected to have stellar masses
(top panels) or SFRs (bottom panels) above a limit, as a function of that
limit. This is shown for Flares (orange) and Shark (blue) at $z=6,\,8,\,10$,
as labelled. The solid lines with connecting dots show the fraction of
galaxies detected in the F200W filter of NIRCam on JWST, while the dashed
lines with connecting squares shows this for galaxies detected in the F775W
filter of ACS on HST. At $z\geq 8$, HST is detecting $<10$% of even the most
massive, or highest SFR galaxies in both models.

Figure 3 shows the stellar mass functions (SMF) and SFR functions (SFRF) for
galaxies at redshifts $z=6,8,10$ in both Flares and Shark. Note that here we
present the predictions assuming no errors in stellar mass or SFRs. It is
common, however, for predictions to be presented assuming random errors for
stellar mass and SFR, which due to the Eddington bias, tend to shift the high-
mass end towards higher stellar masses or SFRs (e.g. Lagos et al. 2018). We
use the instantaneous SFRs of the galaxies in the results presented for both
simulations. For Flares galaxies we follow Lovell et al. (2021) and define the
stellar mass as the total mass of star particles within a 30 kpc aperture
centred on the potential minimum of the subhalo, as such the SFR of Flares
galaxies is taken to be within the same aperture. We use these definitions
throughout this work.

An important effect for which we must account that motivates our particular
selection is the stellar mass resolution of simulated galaxies. Schaye et al.
(2015) and Furlong et al. (2015) suggest that galaxies are sufficiently well
sampled provided that they consist of at least 100 star particles, which for
the star mass resolution of EAGLE and therefore Flares, results in a mass
resolution of $\mathrm{M_{\star}^{lim}\approx 10^{8.2577}M_{\odot}}$ in the
simulation. In a similar fashion, reasonably sampled dark matter halos include
100 dark matter particles in the SURFs simulations (Elahi et al., 2018) that
precipitates into a galaxy stellar mass resolution of
$\mathrm{M_{\star}^{lim}\approx 10^{8}M_{\odot}}$ when considering the gas
seeding and star formation recipes in Shark (Lagos et al., 2018; Elahi et al.,
2018). We thus show the distributions of galaxies with $\mathrm{M_{\star}\geq
10^{8.5}M_{\odot}}$ with the dashed, coloured lines.

Below $\mathrm{M_{\star}\sim 10^{8.5}M_{\odot}}$, and therefore
$\mathrm{SFR\sim 1\,M_{\odot}yr^{-1}}$, the simulations diverge at all
redshifts. This divergence is below the mass resolution of the simulations, so
it is likely artificial. At $z=6$, the simulations agree up to
$\mathrm{M_{\star}\sim 10^{10}M_{\odot}}$ and $\mathrm{SFR\sim 10^{1}\to
10^{2}M_{\odot}yr^{-1}}$. For higher stellar masses and SFRs the simulations
are in tension, with Flares predicting an excess in the number densities
compared to Shark. At $z>6$, the maximum stellar masses and SFRs at which the
simulations agree is generally lower than at later redshifts. At $z=8$ the
simulations agree up to $\mathrm{M_{\star}\sim 10^{9.0}M_{\odot}}$ and
$\mathrm{SFR\sim 10^{1}M_{\odot}yr^{-1}}$, an order of magnitude lower than
the agreement threshold at $z=6$. By $z=10$ the stellar mass distributions of
the simulations are in tension over all stellar masses, with Flares predicting
a higher number density than Shark across the whole stellar mass range probed.
The extension to higher SFRs and stellar masses in Flares is due to the larger
effective volume of Flares compared to Shark. Interestingly, at $z=10$ the
distribution of SFRs agree over a similar SFR range as that at $z=8$. This
indicates that Shark assembles less stellar mass than Flares for a SFR in the
range where they agree. In section 6 we discuss the causes behind the
differences at the high mass end of the SMF between the two simulations. The
solid cyan and red lines show the SMF and SFRF for the galaxies that will be
detected by JWST at this wavelength. In Shark and Flares there is remarkable
agreement between the distributions calculated for the total population and
the JWST detected population above the resolution limit in the simulations;
this means that both simulations predict that the JWST will be able to observe
all galaxies of $\mathrm{M_{\star}\geq 10^{8.5}M_{\odot}}$ and
$\mathrm{SFR\geq 10^{0}M_{\odot}yr^{-1}}$ out to at least $z=10$.

The top panels of Figure 3 show observational inferences of the SMF. In
particular, the $z=6$ SMF of (González et al., 2011) who use a combination of
rest-frame optical, ultraviolet and infrared fluxes obtained from the HST and
Spitzer Telescope to derive stellar masses; the SMFs at redshifts $z=6$ and
$z=8$ from Song et al. (2016) who also use a combination of HST and Spitzer
observations, which cover a greater area than González et al. (2011); the
redshift $z=6,8,10$ SMFs of Stefanon et al. (2021) use a combination of HST
and deep IRAC observations. Both simulations generally agree with the
available observations within the uncertainties. At the highest redshifts, the
observations exhibit large uncertainties meaning that with current
observations, the high redshift SMF cannot be well constrained. A good example
of the poor constraining power of these observations, is that even though
Shark and Flares predict different number densities of massive galaxies at
$z=8$, different sets of observations appear to prefer one or the other model.
Critically, these observations at these extreme redshifts are almost entirely
derived from short wavelength observations that are sensitive to the methods
that are used to correct for dust attenuation. Additionally, some estimates
require near infrared measurements to make use of the full galaxy SED to
determine the stellar mass for the highest redshift objects, which is
troublesome with previous observatories, such as Spitzer, as a result of poor
NIR sensitivity. This type of systematic uncertainty is not accounted for in
the error bars of the observations in Figure 3, nor in the simulations.

A similar story is also true for the observational results shown for the SFRF
in the lower panels of Figure 3. We show results from Katsianis et al. (2017)
who use the UV LFs of Bouwens et al. (2015) to derive dust-corrected SFRs.
Katsianis et al. (2017) use the $\mathrm{IRX-\beta}$ of Meurer et al. (1999)
and the linear $\mathrm{\beta-M_{UV}}$ of Bouwens et al. (2012) to determine
the absorption at $1600$Å, and thus the optical depth, to correct the UV
luminosities for dust attenuation following the method of Hao et al. (2011).
It can be seen that there is good agreement between these results and the
predictions of Flares and Shark across the range of SFRs probed. Where
possible, we have tried to standardise the cosmologies utilised by the
simulations and observations as differences in cosmologies can introduce
tensions (Croton, 2013); it should however be noted that discrepancies
precipitated by different cosmologies are likely dwarfed by differences
induced by modelling details, light-to-mass derivations and dust corrections
(Speagle et al., 2014). It can also be seen that the agreement with these
observational results persist for the galaxies below our chosen stellar mass
selection. However, this agreement is sensitive to the dust corrections used
to derive SFRs.

Figure 5: $\mathrm{SFR-M_{\star}}$ between redshift $z=5\to 10$ for entire
population in Flares (orange) and Shark (blue). The dashed lines with markers
show the median trend and the coloured shaded regions show the $5^{\rm
th}-95^{\rm th}$ percentiles. The solid orange and blue lines show the fit to
the SFS as per Equation 10 for the population with masses above
$\mathrm{10^{8.5}M_{\odot}}$ in Flares and Shark respectively. The red and
light blue solid lines show the fit for the JWST detected population in Flares
and Shark respectively. The grey shaded regions highlights the stellar mass
resolution limit of $\mathrm{10^{8.5}M_{\odot}}$.

Figure 4 shows the fraction of galaxies brighter than 29 magnitudes above a
certain mass or SFR. We calculate this fraction for galaxies observed with the
2$\mu$m F200W filter of NIRCam on the JWST and the 0.75$\mu$m F775W filter of
the Advanced Camera for Surveys on the HST. Our choice of HST filter for
comparison is motivated by its use as the detection band of HST deep fields
(e.g., Beckwith et al., 2006a).

We see that the fraction of JWST detected galaxies is greater than HST
detected galaxies for almost all stellar mass and SFR limits and at all
redshifts, showing that the rest-wavelengths traced by JWST are much less
affected by dust. At $z=10$ there are virtually zero HST-detected galaxies
predicted by Flares at all stellar masses and SFRs, while only $\sim 10\%$ of
star-forming galaxies (SFR $\mathrm{\gtrsim 10M_{\odot}yr^{-1}}$) are detected
in Shark.

Flares predicts that JWST can detect $\gtrsim 90\%$ of galaxies with
$\mathrm{M_{\star}\gtrsim 10^{8.3}M_{\odot}}$ and $\mathrm{SFR\gtrsim
10^{0.5}M_{\odot}yr^{-1}}$ at $z=10$, whereas Shark predicts that JWST can
detect $\sim 70\%$ above those limits. So, in both simulations a significant
population of undetected JWST galaxies do indeed exist. Though the fraction of
galaxies missed by JWST is not significant enough to hinder a thorough
exploration of the $\mathrm{SFR-M_{\star}}$ as is the case with the HST. The
fact that these fractions are varied between the simulations boils down to
different prescriptions of stellar mass assembly and the generation of
synthetic photometry. The JWST should be able to detect sufficiently many
galaxies to distinguish between these model dependent prescriptions.

##  4 Star forming main sequence

In this section, we present our predictions on the shape and stellar mass-
dependent scatter of the SFS at $z\geq 5$.

Figure 6: Redshift evolution of the parameters of the SFS as per Equation 10;
the normalisation, $\mathrm{S_{0}}$, is calculated at
$\mathrm{10^{9.5}M_{\odot}}$ instead of the maximum value to which Equation 10
approaches at high stellar masses. The dashed lines with shaded regions show
the results and $\mathrm{1\sigma}$ uncertainties for the total population of
galaxies in Flares (orange) and Shark (blue). The solid lines with error bars
show the results and $\mathrm{1\sigma}$ uncertainties for the JWST detected
galaxies in Flares (red) and Shark (cyan). The grey points with error bars
show the observational results and uncertainties obtained from Thorne et al.
(2021).

###  4.1 Fitting the SFS

To fit our simulation’s SFS, we follow the same procedures outlined in Thorne
et al. (2021) who use a two-component power law. The latter allows to fit a
possible turn over in the SFS, which, when present, appears in massive
galaxies. Lovell et al. (2021) presented evidence of a turn over in the SFS in
Flares. The functional form of the SFS is thus

$\mathrm{\mathnormal{f}(\mathcal{M})=S_{0}-\log_{10}\left[\left(\frac{10^{\mathcal{M}}}{10^{M_{0}}}\right)^{-\alpha}+\left(\frac{10^{\mathcal{M}}}{10^{M_{0}}}\right)^{-\beta}\right]},$
(10)

where $\mathrm{\mathcal{M}=\log_{10}(M_{\star}/M_{\odot})}$, $\mathrm{S_{0}}$
is the limiting SFR of the function at high stellar mass in units of
$\mathrm{\log_{10}(M_{\odot}yr^{-1})}$, $\mathrm{M_{0}}$ is the stellar mass
in units of $\mathrm{\log_{10}(M_{\odot})}$ at which the SFS turns over and
$\mathrm{\alpha}$ and $\mathrm{\beta}$ are the respective low and high stellar
mass slopes. Equation 10 from Thorne et al. (2021) is itself an adaptation of
the two component model used in Lee et al. (2015) that, instead of assuming
that the SFS flattens to a slope of zero at high stellar masses as is the case
in Lee et al. (2015), allows for an additional degree of freedom, $\beta$, to
describe the high stellar mass slope. We use Equation 10 as the foundation of
a Bayesian fitting routine. The likelihood function that we maximise is a
student-t distribution with a fixed scale of $0.3$. We use broad uniform
priors on the parameters of the SFS: $\mathrm{S_{0}\sim U(0.01,4.01)}$;
$\mathrm{M_{0}\sim U(8.0,19.0)}$; $\mathrm{\alpha\sim U(0.5,2.0)}$;
$\mathrm{\beta\sim U(0.1,0.61)}$. We are thus hypothesising that the SFS does
in fact exist out to $z=10$ and that its shape is consistent with our
function. We estimate the parameters and $\mathrm{1\sigma}$ uncertainties by
performing a Markov-Chain-Monte-Carlo (MCMC) analysis using the Python package
EMCEE (Foreman-Mackey et al., 2013). We have elected to not perform any kind
of SFR or specific SFR
($\mathrm{\log_{10}\left[\frac{SFR}{M_{\star}}yr^{-1}\right]}$) selection in
defining the SFS here, as is the case in many works that focus on the SFS at
lower redshifts in the literature (Davies et al., 2016, 2019, 2022; Katsianis
et al., 2019). Our choice is motivated by the fact that we do not expect a
significantly quenched and off-SFS galaxy population contained in our high
redshift sample. We note, however, that there are passive galaxies in Flares
up to $z\sim 8$ though they constitute $\mathrm{<3\%}$ of the total population
(Lovell et al. in prep).

Figure 5 shows the $5^{\rm th}-95^{\rm th}$ percentiles and medians between
redshifts $z=5\to 10$ for the total galaxy populations in both simulations. It
is clear that both simulations predict the existence of a tight SFS up to
$z=10$, indicating that the process of self-regulation in galaxy growth is
present and efficient in the very early universe. Despite differences in the
stellar mass assembly, the two simulations predict similar SFSs. Flares
exhibits a slightly clearer high stellar mass turn over than Shark. Shark
exhibits approximately uniform variance in the SFS with stellar mass, whereas
Flares shows more puffed up variance at low and high stellar masses and a
minimum in the variance around $\mathrm{\sim 10^{9}M_{\odot}}$ particularly at
$z\lesssim 7$. Note, however, that the increase in the variance at low stellar
masses in Flares happens mostly at the regime where resolution is expected to
significantly affect our results. Therefore, the increased scatter is likely
dominated by star formation stochasticity and requires higher resolution
simulations to probe robustly. We return to the subject of variance about the
SFS in Section 4.2. The skew in the SFS at high stellar masses seen in the
$5^{\rm th}-95^{\rm th}$ percentiles is due to the population of massive
galaxies in Flares undergoing quenching, the cause of which is likely energy
injection from their AGN driving a suppression of star formation in massive
galaxies.

We also show the fitted relations as solid lines, and it can be seen that the
fits agree well with the predicted relations. We only fit galaxies with
$\mathrm{M_{\star}\geq 10^{8.5}M_{\odot}}$. We do see that the lower stellar
mass slope, $\alpha$, is steeper in Flares than it is in Shark at $z\lesssim
7$, which may originate from the highly overdense, starburst galaxy
populations that Flares samples but Shark does not; this is likely to be
strongly related to the enhancement in stellar mass and SFR seen in Figure 3.

Most notably at $z\lesssim 7$ the Flares fits are able to better capture the
full extent of the turn over in the SFS observed in the median trend, with the
median showing a slightly higher normalisation around $\mathrm{\sim
10^{9.5}M_{\odot}}$ than the fit. This is likely caused by quenching that is
driving a cessation of star formation in massive galaxies. The fact that this
turn over is not seen in Shark galaxies below $\mathrm{M_{\star}\sim
10^{10}M_{\odot}}$ indicates that massive galaxy quenching occurs earlier in
Flares (Lovell et al. in prep). Above $z\sim 7$ the SFSs of both simulations
are remarkably similar, indicating that self-regulation is occurring in a
similar fashion between these two simulations.

We have calculated the same parameters of the SFS for the population that will
be detected by the JWST in both simulations. We have included Gaussian
uncertainties of $0.2$ dex and $0.4$ dex on the stellar masses and SFRs of the
JWST populations respectively to account for potential uncertainties in
deriving those quantities from actual observations. We investigate the effect
of uncertainties more closely in Appendix A. Remarkably, the JWST population
is very closely congruent with the total population in the simulations over
the entire $\mathrm{SFR-M_{\star}}$ plane at all redshifts shown in Figure 5.
This indicates that planned JWST observations are deep enough to recover the
SFS.

Redshift | $\mathrm{S_{0}}$ | $\mathrm{M_{0}}$ | $\alpha$ | $\beta$ |
---|---|---|---|---|---
Flares |  |  |  |  |
5.0 | $2.40\pm 0.77$ | $10.38\pm 0.75$ | $1.31\pm 0.16$ | $0.46\pm 0.17$ |
6.0 | $2.18\pm 0.77$ | $10.01\pm 0.77$ | $1.36\pm 0.23$ | $0.44\pm 0.17$ |
7.0 | $2.32\pm 0.92$ | $10.17\pm 0.97$ | $1.28\pm 0.29$ | $0.47\pm 0.17$ |
8.0 | $2.18\pm 1.00$ | $10.09\pm 1.14$ | $1.22\pm 0.36$ | $0.47\pm 0.17$ |
9.0 | $2.11\pm 1.04$ | $9.88\pm 1.23$ | $1.30\pm 0.42$ | $0.42\pm 0.17$ |
10.0 | $1.98\pm 0.98$ | $9.58\pm 1.33$ | $1.06\pm 0.42$ | $0.40\pm 0.17$ |
Flares-JWST |  |  |  |  |
5.0 | $2.02\pm 0.71$ | $9.90\pm 0.72$ | $1.40\pm 0.21$ | $0.48\pm 0.18$ |
6.0 | $2.35\pm 0.79$ | $10.11\pm 0.76$ | $1.37\pm 0.23$ | $0.45\pm 0.17$ |
7.0 | $2.28\pm 0.85$ | $9.98\pm 0.89$ | $1.33\pm 0.29$ | $0.41\pm 0.17$ |
8.0 | $1.93\pm 0.96$ | $9.69\pm 1.27$ | $1.12\pm 0.40$ | $0.46\pm 0.17$ |
9.0 | $2.07\pm 0.98$ | $9.73\pm 1.49$ | $0.97\pm 0.41$ | $0.39\pm 0.17$ |
10.0 | $2.22\pm 0.95$ | $9.69\pm 1.27$ | $1.09\pm 0.43$ | $0.43\pm 0.17$ |
Shark-JWST |  |  |  |  |
5.0 | $3.90\pm 0.32$ | $12.62\pm 0.40$ | $0.93\pm 0.02$ | $0.37\pm 0.24$ |
6.0 | $3.36\pm 0.88$ | $11.93\pm 1.11$ | $0.92\pm 0.05$ | $0.50\pm 0.21$ |
7.0 | $2.35\pm 0.94$ | $10.58\pm 1.20$ | $0.92\pm 0.07$ | $0.62\pm 0.17$ |
8.0 | $2.67\pm 1.03$ | $10.83\pm 1.34$ | $0.90\pm 0.10$ | $0.59\pm 0.17$ |
9.0 | $2.42\pm 0.91$ | $10.25\pm 1.12$ | $0.98\pm 0.17$ | $0.50\pm 0.18$ |
10.0 | $2.05\pm 0.97$ | $9.72\pm 1.38$ | $0.89\pm 0.27$ | $0.51\pm 0.17$ |
Shark |  |  |  |  |
5.0 | $3.84\pm 0.19$ | $12.22\pm 0.19$ | $1.01\pm 0.01$ | $0.23\pm 0.14$ |
6.0 | $3.51\pm 0.94$ | $11.75\pm 1.05$ | $1.01\pm 0.09$ | $0.37\pm 0.21$ |
7.0 | $2.77\pm 0.85$ | $10.73\pm 0.95$ | $1.05\pm 0.11$ | $0.56\pm 0.19$ |
8.0 | $1.76\pm 0.73$ | $9.44\pm 0.80$ | $1.18\pm 0.17$ | $0.61\pm 0.16$ |
9.0 | $1.70\pm 0.84$ | $9.30\pm 0.98$ | $1.14\pm 0.24$ | $0.53\pm 0.17$ |
10.0 | $1.86\pm 0.82$ | $9.30\pm 0.91$ | $1.24\pm 0.29$ | $0.49\pm 0.17$ |

Table 1: Redshift evolution of the SFS fit parameters and associated
$1-\sigma$ uncertainties as described in Section 4 and shown in Figure 6.

Figure 6 shows the redshift evolution of the parameters of the SFS for Shark
and Flares, for both the total and JWST galaxy detected populations. Tabular
data of these fits is recorded in Table 1. The shaded regions and error bars
show $\mathrm{1\sigma}$ uncertainties on the parameter estimates from the MCMC
sampling for the total and JWST detected populations respectively. It is
interesting that above $z\sim 7$ the simulations predict very similar SFS
considering the different modelling processes, and that only $z=0$
observations were used for the tuning of free parameters. The normalisation
evolves similarly in the simulations for all considered redshifts. Below
$z\sim 7$, the fitting to the SFS in Shark indicates a turn over mass that is
an order of magnitude higher mass than that in Flares and a shallower low
stellar mass slope. The fitted turn over mass in Shark is actually higher than
the most massive galaxies that the simulations predict, indicating that there
is no turn over in the SFS in Shark at $z\gtrsim 5$. This, in combination with
the low stellar mass slope tending to unity, suggests that the SFS in Shark is
better described by a single power law at $z\gtrsim 5$. The main difference
between the SFS in the simulations is thus the turn over as a result of
quenching, and the redshift at which that quenching becomes prevalent. We
elaborate on this in Section 6.

Shown as well are observationally derived estimates on SFS parameters from
Thorne et al. (2021) who fit the SFS using the same two-component power law in
20 redshift bins with an equal width of $\mathrm{\sim 0.75Gyrs}$. SFRs and
stellar masses are derived by SED fitting galaxies between
$\mathrm{\lambda\sim 0.154\to 504\mu m}$ in the D10-COSMOS field of DEVILS
(Davies et al., 2018, 2021). These observations have a coarse temporal
resolution at high redshift (with the highest redshift bin being $z=5\to 9$)
and a high stellar mass cut of $\mathrm{\gtrsim 10^{10.4}M_{\odot}}$ for $z>5$
meaning the range of stellar masses probed is small. Despite these caveats,
the simulations agree somewhat with the observations. The SFS normalisation is
consistent with the simulations at redshift $z\lesssim 9$. The turn over mass
is more consistent with Shark than Flares, with Flares being lower, however,
the turnover mass in Shark happens above the range of simulated stellar masses
and hence is not well constrained. This suggests that the observations favour
a single power law SFS at $z\gtrsim 5$, which was a similar conclusion reached
in Thorne et al. (2021) who found negligible turn overs beyond $z\sim 2$. To
assert that the observations favour one model or the other is disingenuous,
however, because the lack of observations here hinder a rigorous exploration
of this parameter space. Furthermore, mass completeness limit used in Thorne
et al. (2021) is $\mathrm{\gtrsim 10^{10.4}M_{\odot}}$ for $z>5$ meaning the
range of stellar masses probed is small.

We thus predict that the JWST will be able to sufficiently detect all
intrinsically bright galaxies and offer for the first time reliable fits to
the SFS for the first galaxies to have formed in the Universe. It is
encouraging then that we predict that the JWST will be able to recover the SFS
for redshifts up to $z=10$ as these observations will be necessary to decouple
slight differences in the simulations. This will elucidate the origins of
these kinds of fine distinctions between the two simulations (and potentially
other simulations), which are embodied in galaxy formation theory.

###  4.2 Stellar mass dependence of the SFS scatter

The usefulness of the SFS is that it encodes information about astrophysical
process that may be driving the shape and normalisation. In particular, the
stellar mass dependent scatter is an indicator of how effective different
feedback mechanisms are in regulating star formation in galaxies. At lower
redshifts, $z\lesssim 0.7$, the $\mathrm{\sigma_{SFS}-M_{\star}}$ has been
observed to exhibit a minimum vertex parabolic shape, with its turning point
in the neighbourhood of $\mathrm{\sim 10^{9}M_{\odot}}$ (Davies et al., 2019,
2022). The minimum point of the stellar mass dependent scatter is interpreted
as the galaxy phase at which gas inflow is balanced by both star formation and
feedback. Galaxies in this phase experience little variation about the SFS,
and so minimal scatter, as gas compaction and subsequent depletion events are
highly self regulated for these galaxies (Tacchella et al., 2016).

Figure 7: $\mathrm{\sigma_{SFS}-M_{\star}}$ relation for the total galaxy
populations in Flares (orange) and Shark (blue) at redshifts $z=5\to 10$. Also
shown is the same relation but for the JWST detected galaxies in Flares (red)
and Shark (cyan). The grey shaded regions shows the stellar mass resolution
limit, $\mathrm{M_{\star}=10^{8.5}M_{\odot}}$.

Figure 7 shows the $\mathrm{\sigma_{SFS}-M_{\star}}$ for Shark and Flares,
where $\mathrm{\sigma_{SFS}}$ is the $16^{\rm th}-84^{\rm th}$ percentile
range. It can be seen that both simulations are fairly consistent about the
shape of the relation with stellar mass over most stellar masses and redshifts
probed. Below $\mathrm{M_{\star}=10^{8.5}M_{\odot}}$ the simulations disagree,
with Shark predicting a decrease in the scatter at $z\gtrsim 9$ and a uniform
scatter at $z\lesssim 8$, while Flares predicts a significant uptick of the
scatter at all redshifts. The decreased scatter in low stellar mass Shark
galaxies at $z\gtrsim 9$ may suggest that they are better self regulated as
most are not far off from the SFS. Stochastic star formation then becomes more
important at $z\lesssim 8$ to drive up the scatter. Conversely the uptick seen
in Flares throughout all redshifts is likely an effect of stochastic star
formation due to the poor resolution at these masses. For both simulations,
these effects occur below the stellar mass resolution limit of the simulations
however, and so these results cannot be considered robust.

A striking feature of the $\mathrm{\sigma_{SFS}-M_{\star}}$ is the lack of
significant scatter at the massive end, $\mathrm{M_{\star}\gtrsim
10^{10}M_{\odot}}$, that both simulations predict at redshifts $z=5\to 10$.
The likely reason for this is that the massive galaxies at these redshifts in
the simulations have not yet had enough time to sufficiently grow their
supermassive black holes. AGN feedback is thought to increase the asymmetric
variance from the SFS as feedback continuously inhibits star formation and
drives massive galaxies below the SFS (Davies et al., 2022). This is
consistent with other studies that investigate this effect in simulations
(Katsianis et al., 2019). Although there is a general lack of significant
scatter at these large stellar masses, with decreasing redshift, Flares starts
to show an increase in $\mathrm{\sigma_{SFS}-M_{\star}}$ at stellar masses
$\gtrsim 10^{10.3}\rm M_{\odot}$ at $z\sim 5$, which shows the initiation of
quenching in these galaxies. This is consistent with lower redshift results
from simulations that show that the scatter in the SFS for galaxies
$\mathrm{M_{\star}\gtrsim 10^{10}M_{\odot}}$ increases by 0.05 dex from z=5 to
z=0 (Matthee & Schaye, 2019). This is not seen in Shark. This difference
between the simulations together with the JWST-like samples following the same
SFS as all the simulated galaxies, indicates that the JWST will be able to
place strong constraints on the onset of quenching in massive galaxies from
the characterisation of the SFS at high redshift.

While there is the issue of stellar mass resolution we can conclude that faint
dwarf galaxies below the detection sensitivity of the JWST are likely, if not
completely, responsible for the disparity between the total and JWST detected
populations observed for stellar masses $\mathrm{M_{\star}\lesssim
10^{8.5}M_{\odot}}$. JWST will not detect many galaxies with low star
formation rates, as was confirmed in Figure 3. This means that the
contribution to the scatter about the SFS for low stellar mass objects is
biased to those on or above the SFS, which would lead to an underestimation of
the scatter. This trend is seen in Figure 7 for both simulations, although
again this is the mass range that falls below the resolution of the
simulations.

Most importantly, above stellar masses of $\mathrm{10^{8.5}M_{\odot}}$ the
JWST detected population in each simulation is in agreement with the total
population up to the highest stellar masses probed and for all redshifts shown
up to $z=10$. We thus predict that the JWST should be able to sufficiently
detect all intrinsically bright galaxies to understand the
$\mathrm{\sigma_{SFS}-M_{\star}}$ of the Universe up to redshift $z=10$.

##  5 Cosmic stellar mass and star formation rate history

Figure 8: The CSMH (left) and CSFRH (right) at redshifts $z=5\to 10$. The
light coloured, dashed lines show these quantities for the entire population
of galaxies in Flares (orange) and Shark (blue). The solid lines show the
results for the population that will be detected by the JWST in Flares (red)
and Shark (cyan). The dot dashed grey lines show the fits to the CSFRD and
CSMD from Madau & Dickinson (2014). Data points with error bars in the left
panel show observational results of the CSMH (Yabe et al., 2009; González et
al., 2011; Lee et al., 2012; Labbé et al., 2013; Stefanon et al., 2021). In
the right panel observational estimates on the CSFRH (Bouwens et al., 2012;
Schenker et al., 2013; Bouwens et al., 2020) are shown as data points with
error bars.

Figure 8 shows the cosmic stellar mass and cosmic star formation rate
densities (CSMD/CSFRD) as a function of redshift, or cosmic stellar
mass/cosmic star formation rate history (CSMH/CSFRH). The dashed lines show
the results for all galaxies in the simulations. Both simulations predict
consistent shapes of both the CSMH and CSFRH over all redshifts shown. Flares,
however, exhibits a noticeably higher normalisation of both the CSMH and CSFRH
compared to Shark. Flares is $\mathrm{\sim 4.4}$ ($\mathrm{\sim 2.0}$) times
higher in CSMD (CSFRD) than SHARK averaged over redshift, but up to
$\mathrm{\sim 9.6}$ ($\mathrm{\sim 2.9}$) times larger at $z=10$. The
difference in CSMD becomes systematically smaller with decreasing redshift,
with Flares being $\mathrm{\sim 2}$ times larger than Shark by $z=5$. The
difference in CSFRD is lowest at $z=7$ where Flares is $\mathrm{\sim 1.6}$
times larger than Shark but becomes $\mathrm{\sim 1.9}$ times larger by $z=5$.
This is potentially a result of the stellar mass and SFR enhancement captured
by Flares as a result of the contribution from overdense environments sampled
in Flares that are not present in Shark. It can also be seen that the
difference in the normalisation between the two simulations is larger in the
CSMH than in the CSFRH, which appears to be the result of the different
chemical enrichment models used in the simulations. Shark instantaneously
recycles and enriches the interstellar medium with metals following star
formation (Lagos et al., 2019), while there is a delay between star formation
and enrichment in Flares to account for the evolution of type Ia and II
supernovae and AGB stars (Schaye et al., 2015; Crain et al., 2015). We also
inspected the gas phase metallicities in the simulations and found that Shark
predicted higher metallicities per unit stellar mass and SFR at $z=5$ and
$z=10$ than Flares. This can explain the larger offset in the CSMH, as more
mass is diverted from stars toward metal pollution of the ISM in Shark than
Flares, particularly at earlier times. The fact that the differences in the
CSMH become smaller at later times is the result of the instantaneous
recycling approximation becoming a better approximation as the dynamical time
of halos becomes longer.

To ensure that resolution effects are not biasing these results, we also
examined the CSFRH and CSMH only for galaxies that satisfy
$\mathrm{10^{8.5}M_{\odot}<M_{\star}\leq M^{9.0}M_{\odot}}$. The reason for
this particular selection is to be above the resolution limit but below the
stellar mass range where Flares predicts a much greater number density of
stellar mass and SFR compared to Shark, which would bias our results. We found
that the same trends as the total population persist in this case except at
$z=10$ where the difference in normalisation between the CSMD and CSFRD are
similar. As such, while chemical enrichment is an important factor in
influencing the normalisation of the CSMH and CSFRH, we cannot rule out
different star formation histories also somewhat driving this difference.

We have also included some observational estimates on both the CSMH and CSFRH.
The observations on the right panel of Figure 8 show the CSFRH derived by
Madau & Dickinson (2014) who curated UV LFs of Bouwens et al. (2012) and
Schenker et al. (2013), and then converted them into a CSFRH. The
observational data points shown in the left panel are UV-derived estimates of
the CSMD (Yabe et al., 2009; González et al., 2011; Lee et al., 2012; Labbé et
al., 2013; Stefanon et al., 2021). Where possible, we have standardised the
cosmologies and IMFs used for these observational results to be most
comparable with our simulations. The observations for both the CSMH and CSFRH
are consistent with the shapes of the predicted relations from both
simulations up to $z\sim 10$ and $z\sim 8$ for the CSMH and CSFRH
respectively. The Shark predictions of the CSMH and CSFRH are in better
agreement with the observations, with Flares generally predicting a higher
normalisation than the observations in both the CSMH and CSFRH over the
relevant redshift range. We must again stress, however, that all of these
observational estimates are entirely from UV sources that have been corrected
for dust attenuation using the infrared-excess technique, and are thus subject
to the uncertainties of the method (e.g., Shivaei et al., 2020). We showed in
Figure 4 that HST-dark systems are significant in number in the simulations up
to $z\sim 10$, implying that dust must be an important factor when determining
SFR from rest-frame UV observations.

The solid lines show the cosmic densities for the population that will be
detectable with the NIRCam F200W filter of the JWST. We only show the result
with a photometric cut, truncating at $\mathrm{29\,ABmag}$, and not a cut in
stellar mass. Noticeably, both simulations predict that the JWST detected
populations do not account for the entire distribution of stellar mass and SFR
in the universe over all these redshifts. Perhaps unsurprisingly, this
indicates that there are faint galaxies that fall below the detection
sensitivity of the JWST and are not accounted for in the calculation of the
CSMH and CSFRH. Fortunately, this is well understood as the Malmquist bias for
which there are numerous correction methods (Weigel et al., 2016). As only the
normalisation, and not the shape, of the CSMH and CSFRH for the JWST detected
population is different from that of the total population it shows that the
JWST should only be sensitivity limited, and not miss a significant population
of heavily obscured systems compared to other shorter wavelength instruments
such as the HST.

##  6 Discussion

Figure 9: Top: $\mathrm{M_{subhalo}-M_{\star}}$ at $z=6,8,10$, as labelled,
for Flares (orange) and Shark (blue). The dotted lines show the median
quantity while the shaded regions show the 5-95 quantiles. We only show
central galaxies here. Bottom: $\mathrm{sSFR-M_{subhalo}}$. The colours and
labels are the same as the top row. We only show central galaxies here.

We have shown that both simulations, Flares and Shark, predict that the JWST
will be able to observe galaxies with a large enough range of stellar masses
and SFRs to allow a thorough study of SFS, including the onset of quenching in
massive galaxies, and thus stellar mass assembly during the early Universe.
Additionally, we predict that the JWST will be a vital tool to better
understand cosmic reionisation as a large census of galaxies that span a high
dynamic range of stellar masses and star formation rates will will likely
encompass ionising sources. Prior observational studies (e.g., Thorne et al.,
2021) indicate the prevalence of a SFS around $z\sim 5$, and here we show that
both the simulations analysed here predict the existence of a SFS up to at
least $z=10$. The existence of the SFS is a product of self-regulation of star
formation in Flares and Shark where galaxies relax to the locus of the SFS
through a balance of star formation and feedback. Each of these simulations
implement these prescriptions differently so it is an interesting result that
self-regulation appears in both independently. Flares for example uses a gas
density floor above which stars form, whereas this is not the case for Shark.
Furthermore, Flares implements stellar feedback stochastically to neighbouring
SPH gas particles directly (Dalla Vecchia & Schaye, 2012), whereas stellar
feedback is implemented with the calculation of the mass-ejection rate in
Shark, which depends on the maximum circular velocity of the galaxy (e.g.,
Lagos et al., 2018). We show in Section 6.2 that the implementation of
feedback from AGN between the simulations can explain the different shapes of
the SFS. In this work, we quantified the predicted SFS’s shape and scatter as
a function of stellar mass in each simulation. We expect the JWST to be able
to test these predictions in the near future. Figure 3 shows that both Flares
and Shark predict similar intrinsic distributions of stellar mass and SFR
which indicates that, although these simulations differ in many details, the
sub-grid methods of stellar mass assembly used in each must act similarly. The
relative agreement with the existing observations indicate that our
understanding of stellar mass assembly must not be so far removed from how the
process functions in nature; though, it will be a task for the JWST to confirm
this suspicion with updated, near infrared photometry.

Despite this overall agreement, there are important differences between the
two simulations that we discuss below.

###  6.1 The abundance of massive galaxies

A glaring difference between the simulations is the excess number density of
high stellar mass galaxies with high star formation rates at fixed stellar
mass in Flares compared to Shark. To get to the bottom of this difference, we
show in the top panels of Figure 9 the subhalo-stellar mass relation for both
Flares and Shark at $z=6,8,10$ for central galaxies. We only use centrals in
these calculations to simplify the comparison; including satellite galaxies
only slightly increases the higher mass, $\mathrm{M_{\star}\gtrsim
10^{10}M_{\odot}}$, scatter of these trends. Flares extends the
$\mathrm{M_{subhalo}-M_{\star}}$ 0.2 dex and 1 dex beyond Shark at $z=6$ and
$z=10$ respectively. By the fact that the simulations broadly agree on the
shape of the relation over all redshifts, it is possible to conclude that the
difference in the predicted SMFs is due to Flares covering a wider dynamic
range, extending to rare, large over-densities (i.e. cosmic variance).

Despite the broad agreement between the simulations in the top panel of Figure
9, it is worth highlighting that Flares predicts slightly more stellar mass
for a fixed subhalo mass than Shark, at $\mathrm{M_{subhalo}\gtrsim
10^{12}\,\rm M_{\odot}}$ at $z=6$ and $z=8$. The lower panels of Figure 9 show
the specific SFR as a function of subhalo mass for the simulations at
$z=6,8,10$. Especially at $z=6,8$, Flares predicts a higher sSFR than Shark
for $\mathrm{M_{subhalo}\gtrsim 10^{11}M_{\odot}}$, indicating that star
formation is more efficient in massive halos in Flares than in Shark.
Therefore it is likely a combination of both Flares sampling greater
overdensities than Shark and the star formation prescription that are driving
differences in the resulting stellar mass and SFR distributions.

###  6.2 The onset of quenching in massive galaxies

Another important difference between the simulations is their predictions on
the onset of quenching in massive galaxies, $\mathrm{\geq 10^{10}M_{\odot}}$
as indicated by the prevalence of the turn-over and increased scatter in the
SFS shown in Figures 5, 6 and 7. Davies et al. (2022) interpreted the
increased scatter at the massive end of the SFS for the DEVILS sample as being
caused by AGN feedback inhibiting star formation and driving galaxies below
the SFS. This is corroborated by earlier EAGLE simulations’ results of
Katsianis et al. (2019), who also attributed the increased scatter at the
massive end to AGN feedback.

Figure 10: $\mathrm{\Delta SFS-M_{BH}}$ at redshifts $z=5,6,7,8,9,10$ as
labelled for Flares (orange) and Shark (blue). We use the sSFR per unit black
hole mass shifted up by the median sSFR between $\mathrm{M_{BH}=10^{6.75}\to
10^{7.25}M_{\odot}}$ as a proxy measure for distance to the SFS. The solid
lines indicate the median and the shaded regions indicate 16-84 quantiles. The
dashed and dotted lines show the results at $z=0$ and $z=2$. We show $z=0$ and
$z=2$ results from EAGLE-AGNdT9 (Schaye et al., 2015; Crain et al., 2015;
McAlpine et al., 2016) with brown lines in place of Flares as Flares stops at
$z\approx 5$. Because Flares uses the same code and model as EAGLE-AGNdT9,
this is an appropriate comparison to make. The X symbols indicate the black
hole mass at which $\mathrm{\Delta SFS}$ experiences its greatest decline at
$z=2$, while the diamond symbols indicate the same thing but for $z=0$. Note
that for EAGLE, the Xs and diamonds positions are identical.

Here, we explore much higher redshifts than Davies et al. (2022)
($0.1<z<0.85$) and Katsianis et al. (2019) ($z<4$). To investigate the effect
of AGN feedback on quenching and the increased scatter in the SFS at the high-
mass end we show the distance from the SFS as a function of super massive
black hole (SMBH) mass for the two simulations in Figure 10. We use the sSFR
per unit mass shifted up by the median sSFR between
$\mathrm{M_{BH}=10^{6.75}\to 10^{7.25}M_{\odot}}$ as a proxy measure for
distance to the SFS. We use this measure instead of our fits to the SFS to
promote comparisons with lower redshift results for which we do not have SFS
fits. We see that the $\mathrm{\Delta SFS\sim 0}$ is approximately constant,
in Shark over all black hole masses at all redshifts. We see at $z\lesssim 9$
that the $\mathrm{\Delta SFS}$ between Shark and Flares begins to diverge with
Flares experiencing a downturn in $\mathrm{\Delta SFS}$ for SMBH masses
$\mathrm{M_{BH}\gtrsim 10^{7}M_{\odot}}$. The divergence for
$\mathrm{M_{BH}\gtrsim 10^{7}M_{\odot}}$ is greatest at $z\lesssim 7$ that
lines up well with the similar redshifts where the turn-over and increased
scatter become more prevalent in Flares compared to Shark, indicating a
connection between forming very massive SMBHs and quenching, congruous with
the lower redshift conclusions in the literature. Bower et al. (2017) showed
that in EAGLE, galaxies being quenched by AGN feedback are characterised by a
strongly non-linear SMBH growth phase. The SMBH mass is thus a good predictor
of the regime in which AGN feedback is efficient. Terrazas et al. (2020)
showed that in another cosmological hydrodynamic simulation, Illustris-TNG, a
similar behaviour is seen, and galaxies above a given SMBH mass experience AGN
feedback quenching.

The dashed and dotted lines in Figure 10 show the $z=0$ and ${z=2}$ relations
respectively for EAGLE AGNdT9-50 (Schaye et al., 2015; Crain et al., 2015;
McAlpine et al., 2016) and Shark. We remind the reader that Flares adopted the
same model as EAGLE meaning the latter is the right comparison dataset for
Flares at ${z<5}$. The EAGLE Ref-100 simulation uses a different set of
parameters meaning that it is unsuitable for comparison despite being a larger
volume. At ${z=0}$, both EAGLE and Shark predict a sharp downturn in sSFR
around $\mathrm{M_{BH}\gtrsim 10^{7}M_{\odot}}$ indicating that above this
SMBH mass AGN feedback becomes efficient at inhibiting star formation in
galaxies. This SMBH mass threshold for quenching does not evolve in EAGLE and
Flares, staying around $\mathrm{M_{BH}\gtrsim 10^{7}M_{\odot}}$ at all times;
i.e., as soon as a galaxy hits that SMBH, its star formation activity starts
to quench due to AGN feedback. This is highlighted by the orange X and diamond
symbols that correspond to the black hole mass at which the $\mathrm{\Delta
SFS}$ experiences its sharpest down turn at ${z=2}$ and ${z=0}$ respectively.
Note, that the position of these symbols are identical. However, Shark behaves
very differently, with the SMBH threshold mass above which quenching driven by
AGN feedback happens evolves towards higher SMBH masses as the redshift
increases. By ${z=2}$, the SMBH mass threshold for quenching in Shark is
$\sim$1 dex higher than at ${z=0}$ as indicated by the blue X and diamond
symbols respectively.

In Shark the SMBH accretion rate due to hot halo cooling is what regulates the
efficiency of star formation quenching in central galaxies (Lagos et al.,
2018). We observe a greater variation in turn-over mass over time in Shark
that is a reflection of the AGN entering the appropriate feedback mode at
different SMBH masses at different redshifts. This is not the case in EAGLE
and Flares, and quenching mostly occurs at a fixed mass. This implies that the
AGN in Shark for ${z\geq 5}$ are mostly growing in the QSO mode, and are yet
to have significant accretion rates from the hot-halo cooling mode. This means
they cannot quench galaxies, while all that matters in Flares is that the AGN
have grown to sufficient masses, which happens at ${z\lesssim 7}$. Therefore,
the cause of the difference seen in the SFS at the massive end between the
simulations is a result of different models of AGN feedback. The JWST will be
able to probe the SFS up to $\mathrm{M_{\star}\sim 10^{10}M_{\odot}}$,
providing stringent constraints on theoretical models of AGN feedback.

###  6.3 The volume probed by JWST surveys

So far we have not accounted for the survey volume required to ascertain a
sample representative of the Universe. Uncertainties in pencil beam surveys
such as the Hubble Ultra Deep Field (Beckwith et al., 2006b) are dominated by
cosmic variance where the probability of obtaining a survey sample
representative of the entire Universe scales with area coverage (Driver &
Robotham, 2010).

Figure 11: Top:  Weighted-averaged $\mathrm{\chi^{2}}$ values for each of the
SFS fits where we have resampled galaxy stellar masses and SFRs according to
different survey volumes 20 times for Flares (left) and Shark (right).
$\mathrm{\chi^{2}=\sum^{4}_{i=0}\frac{(x_{i}-x^{true}_{i})^{2}}{x^{true}_{i}}}$,
where $\mathrm{x_{i}}$ represents the derived parameters of the resampled SFS
and $\mathrm{x_{i}^{true}}$ are the derived parameters for the entire
population in the simulations. When calculating the average
$\mathrm{\chi^{2}}$, we weight each of the $\mathrm{\chi^{2}}$ values by the
inverse of the sum of squared uncertainties of the four parameters. The Filled
circles show the $\mathrm{\chi^{2}}$ for the JWST-detected galaxies; upright,
left-pointing and right-pointing triangles show the $\mathrm{\chi^{2}}$ for a
JWST surveys of areas $\mathrm{10arcmin^{2}}$ (i.e. a single pointing),
$\mathrm{220arcmin^{2}}$ (i.e. a medium deep field) and $\mathrm{0.6deg^{2}}$
(i.e. the Web COSMOS survey), respectively, as labelled. In cases where the
number of data points from the sample is less than the number of fit
parameters, e.g., at $z=10$ for single-pointing surveys, we set the averaged
$\mathrm{\chi^{2}}$ to an upper limit of 2. Bottom:  SFS plots at ${z=5}$
(left) and ${z=10}$ (right), for Flares (red and orange colours) and Shark
(blue and cyan colours). We use the same symbols and colour gradient as the
top panels to reflect different fits for each sample. Orange and blue shaded
regions show the $5^{\rm th}-95^{\rm th}$ percentile ranges for all galaxies
in the simulations. We find the average set of parameters for each of our 20
samples, weighting by the inverse of the error of each of the four SFS
parameters. The dashed, black line shows the SFS for our initial guess of
parameters ($\mathrm{[1.0,10.0,1.0,0.5]=S_{0},M_{0},\alpha,\beta}$). The
offsets of some of the fitted relations, e.g., single-pointing surveys at
$z=10$, are due to the fits defaulting to the initial guess because of
insufficient number statistics, which biases the average SFS relation toward
the initial guess.

To better assess the impact of survey area on the ability to constrain the
shape of the SFS, we refit the SFS after resampling the stellar mass and star
formation rate distributions assuming some survey volume. Specifically, we
calculate the 2D histogram of SFR and stellar mass per unit volume and for a
given redshift for JWST-selected galaxies in both simulations, and calculate
the cumulative distribution function per unit volume at
${z=5,\,6,\,7,\,8,\,9\,\mathrm{and}~{}10}$. Multiplying this by a given survey
volume gives the expected cumulative number distribution of galaxies in the
survey at a given redshift. Here, we do not assume any random errors on
stellar masses or star formation rates to aid in comparison between fits from
the same simulation. We repeat this resampling and fitting process twenty
times and average the derived SFS parameters to mitigate randomness caused by
sampling the joint stellar mass and SFR distributions. We calculate
$\mathrm{\chi^{2}=\sum^{4}_{i=0}\frac{(x_{i}-x^{true}_{i})^{2}}{x^{true}_{i}}}$,
where $\mathrm{x_{i}}$ represents the derived parameters of the resampled SFS
and $\mathrm{x_{i}^{true}}$ are the derived parameters for the entire
population in the simulations, to assess how well each survey can reconstruct
the SFS of the true population. We average this over the 20 samples. When
calculating the average, we weight each sample by the inverse of the sum of
squared uncertainties of each of the four parameters. We also impose a
narrower prior on the turn-over mass ($\mathrm{M_{0}\sim U(8,13)}$) when
fitting our samples to prevent highly unphysical solutions to the high-mass
end, especially in cases where high mass data points are missing as is the
case for surveys that cover small areas. We note that the upper bound of this
prior straddles the derived turn-over masses shown in Table 1.

We show the results in Figure 11 for the entire population of galaxies in both
simulations, the JWST-detected galaxies (determined from their F200W apparent
magnitude), galaxies detected in a single pointing of
$\mathrm{10\,arcmin^{2}}$, galaxies detected in $\mathrm{220\,arcmin^{2}}$ and
galaxies detected in $\mathrm{0.6\,deg^{2}}$. The choice of
$\mathrm{220\,arcmin^{2}}$ and $\mathrm{0.6\,deg^{2}}$ survey areas comes from
the JWST Medium Deep Field (WMDF) Survey (Proposal ID: 1176, PI: R. Windhorst)
and the Web COSMOS survey (Proposal ID: 1727, PI: J. Kartaltepe)
222https://www.stsci.edu/jwst/science-execution/program-information.html.

We see that in general it becomes more difficult to constrain the shape of the
SFS with smaller area surveys. This is true in Shark over all redshifts shown,
while only at ${z>7}$ in Flares do single pointings struggle to constrain the
SFS. The fact that all survey areas can constrain the SFS in Flares at ${z\sim
5}$ while Shark can not is a reflection of the two component SFS function
being unsuitable for Shark at these lower redshifts. The wide offsets of the
SFS fits for the average set of parameters of the narrowest surveys at $z=10$
seen in the bottom right panel of Figure 11 are due to the fits essentially
failing and defaulting to the initial guess of the parameters, which biases
the SFS to the initial guess. We do see that both simulations predict that a
Web COSMOS-like survey without any assumed random errors will be able
constrain the true SFS parameters at ${z\geq 5}$. We further explore the
effect of random errors on $\mathrm{M_{\star}}$ and SFR in appendix A.
Nevertheless, we advocate for wider area surveys that will be able to get
accurate SFRs and stellar masses needed to thoroughly study the SFS.

Interestingly, at ${z=6,10}$ in Flares the average $\chi^{2}$ for the Web
COSMOS-like survey is higher than the narrower WMDF-like survey. This is
likely a result of the fitting technique. At ${z=10}$ for example, the turn-
over mass is $\mathrm{10^{10}M_{\odot}}$ for the WMDF-like survey, which goes
beyond the range of resampled galaxies implying that the turn-over is
unphysical for this sample, despite the new imposed prior. For the Web COSMOS-
like survey, we predict more galaxies at higher stellar mass, but not so many
that the turn-over in the SFS can be constrained. In this case, the SFS
favours a single power law fit, and the average $\chi^{2}$ increases. This
just further shows why it is necessary to have accurate SFRs and stellar
masses over the breadth of the $\mathrm{SFR-M_{\star}}$.

Our resampling method does not take into account the underlying structure of
the Universe that will affect our predicted constraints on the SFS. As Flares
is comprised of 40 unique subvolumes that have each been accordingly weighted,
a smooth mass distribution cannot be reconstructed, and thus cannot be used to
test cosmic variance directly. Therefore, to test the effect of cosmic
variance, we limit our attention to Shark, where we can partition the box into
smaller sub volumes of equal weighting and refit the SFS in each. We divide
the whole Shark box in 3 dimensions into 8 independent, contiguous and equal
size subvolumes. We, once again, choose not to use random errors to emphasise
the effect of cosmic variance only.

Figure 12: Left:  $\mathrm{\chi^{2}}$ values (as defined in Figure 11), for
the SFS fit performed in each subvolumes of the total Shark box. Galaxies
selected based on their F200W magnitude, are shown with filled circles. We
show each set of $8$ subvolumes with different symbols, as labelled. Middle
and right:  SFS plots at ${z=5}$ (middle) and ${z=10}$ (right). We use the
same symbols and colour gradient as the left panels to show the SFS fits of
each galaxy sample.

Figure 12 shows the resulting SFS of this exercise. Cosmic variance mildly
impacts the ability to constrain the SFS at all redshifts. At ${z=5}$ cosmic
variance negligibly affects the SFS fit. By ${z=10}$ cosmic variance affects
the whole appearance of the SFS, with the effect on the constraint at high
masses being significant. Therefore, in addition to our advocacy of wider
surveys we further advocate for surveys distributed over a range of fields and
sight lines.

##  7 Conclusions

We have presented predictions on the ability for the JWST to quantify the
details of star formation and stellar mass assembly of galaxies between
${z=5\to 10}$. By using two distinct simulations, Flares and Shark, we have
provided a theoretical framework, unmarred by potential systematic biases
contained in a single simulation, in which the emerging observations of the
JWST can be tested to better understand the astrophysics of the early
Universe.

The key results are:

  * •

Stellar mass and SFR functions:  We predict that throughout redshifts ${z=5\to
10}$, the JWST will be able to detect all intrinsically bright galaxies up to
redshift $\rm z=10$, and thus will be able to parametrise the distributions of
stellar mass and SFRs of the first galaxies in the Universe. Conversely, both
Flares and Shark predict that very dust-obscured galaxies are completely
absent from current HST surveys. This is expected to affect current
measurements of the SMF and SFRF from HST observations (Katsianis et al.,
2017; McLure et al., 2010; Bouwens et al., 2015; Song et al., 2016; González
et al., 2011; Thorne et al., 2021; Stefanon et al., 2021). We predict that the
JWST observations are only limited by the detection sensitivity of the
instrument and survey area, and will not miss highly obscured galaxies.

  * •

Star forming main sequence: We predict that the JWST will be able to
parametrise the shape and the stellar mass dependent scatter of the SFS
between ${z=5\to 10}$ for the first time by obtaining observations over a
sufficient wide range of stellar masses and SFRs. As the SFS encodes
information about the astrophysics driving the shape and scatter, the JWST
will prove to be an invaluable tool in improving our understanding of galaxy
formation up to the onset of cosmic reionisation when the first galaxies
formed, and on the rising of massive, AGN-quenched galaxies towards the end of
reionisation.

  * •

Cosmic stellar and SFR histories: Because we predict the JWST to be able to
observe every intrinsically bright galaxy out to redshift ${z\sim 10}$, these
observations should be able to constraint the CSMH and CSFRH between ${z=5\to
10}$ robustly. The main issue here will be the contribution of low-mass
galaxies, $M_{\star}\lesssim 10^{8}\,\rm M_{\odot}$ that are fainter than
F200W 29mag, and which is expected to be significant at $z\gtrsim 5$. This
could be resolved if the SMF and SFRF are extrapolated.

  * •

The dependence of predictions on the baryon physics models:  By using two
independent simulations we have shown that each prediction on the ability of
the JWST is model dependent. Both Flares and Shark show similar distributions
of stellar mass and SFR indicating that their sub-grid recipes are similar,
but not exact. An important delineator between the two simulations is how each
models dust that is itself sensitive to how each simulation deals with star
formation and chemical enrichment. As we predict that the JWST will be able to
provide a sufficient census of galaxies out to redshift ${z\sim 10}$ we
suggest that the JWST observations, in combination with multi-wavelength
follow up observations, will be able distinguish between the two models.

  * •

The effect of survey area and cosmic variance:  By sampling from our simulated
galaxy populations to obtain particular survey samples we find that wide area
surveys are necessary to constrain the shape of the star-forming main sequence
from $z=5$ to $z=10$. Furthermore, by dividing the Shark box into 8
independent and contiguous subvolumes and refitting the SFS in each we show
that cosmic variance is an important factor in influencing the shape of the
SFS at ${z\sim 10}$, but less so at lower redshifts.

With continued observations, the JWST will open the door to the next milestone
toward our understanding of galaxy formation. We have presented predictions on
the ability of the JWST so that we are well poised to compare them against the
JWST observations as they continue to arrive. A deeper study of the specific
sub-grid recipes used in simulations will allow us to interpret the
astrophysics that the telescope recovers. Although Flares and Shark broadly
agree with each other, it is a key objective of the JWST to distinguish, and
comprehend the origin of, the slight differences between the independent
models.

## Acknowledgements

We thank the anonymous reviewer for their time and constructive comments. JCJD
is supported by an Australian Government Research Training Program (RTP)
Scholarship. We thank the entire Flares team for their support and feedback.
CL has received funding from the ARC Centre of Excellence for All Sky
Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013.
LJMD acknowledges support from the Australian Research Councils Future
Fellowship scheme (FT200100055). CCL acknowledges support from the Royal
Society under grant RGF/EA/181016. Cosmic Dawn Centre is funded by the Danish
National Research Foundation. This work was supported by resources provided by
The Pawsey Supercomputing Centre with funding from the Australian Government
and the Government of Western Australia. This work used the DiRAC@Durham
facility man- aged by the Institute for Computational Cosmology on be- half of
the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by
BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1,
ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant
ST/R000832/1. DiRAC is part of the National e-Infrastructure.

## Data Availability

Figures, scripts and additional data is available upon reasonable request to
the authors. Flares is hosted here https://github.com/flaresimulations and
Subfind/photometry outputs are available here
https://flaresimulations.github.io/. Shark is hosted here
https://github.com/ICRAR/shark. Shark Subfind/photometry outputs are available
upon reasonable request.

## References

  * Abramson et al. (2014) Abramson L. E., Kelson D. D., Dressler A., Poggianti B., Gladders M. D., Oemler Augustus J., Vulcani B., 2014, ApJ, 785, L36
  * Bahé et al. (2017a) Bahé Y. M., et al., 2017a, MNRAS, 470, 4186–4208
  * Bahé et al. (2017b) Bahé Y. M., et al., 2017b, MNRAS, 470, 4186–4208
  * Barnes et al. (2017) Barnes D. J., et al., 2017, MNRAS, 471, 1088
  * Baugh (2006) Baugh C. M., 2006, Reports on Progress in Physics, 69, 3101
  * Beckwith et al. (2006a) Beckwith S. V. W., et al., 2006a, AJ, 132, 1729
  * Beckwith et al. (2006b) Beckwith S. V. W., et al., 2006b, AJ, 132, 1729
  * Bellstedt et al. (2020) Bellstedt S., et al., 2020, Monthly Notices of the Royal Astronomical Society, 498, 5581
  * Benson (2010) Benson A. J., 2010, Phys. Rep., 495, 33
  * Bouwens et al. (2012) Bouwens R. J., et al., 2012, ApJ, 752, L5
  * Bouwens et al. (2015) Bouwens R. J., et al., 2015, ApJ, 803, 34
  * Bouwens et al. (2020) Bouwens R., et al., 2020, The Astrophysical Journal, 902, 112
  * Bouwens et al. (2021) Bouwens R. J., et al., 2021, The Astronomical Journal, 162, 47
  * Bower et al. (2017) Bower R. G., Schaye J., Frenk C. S., Theuns T., Schaller M., Crain R. A., McAlpine S., 2017, MNRAS, 465, 32
  * Brinchmann et al. (2004) Brinchmann J., Charlot S., White S. D. M., Tremonti C., Kauffmann G., Heckman T., Brinkmann J., 2004, MNRAS, 351, 1151
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Cañas et al. (2019) Cañas R., Elahi P. J., Welker C., del P Lagos C., Power C., Dubois Y., Pichon C., 2019, MNRAS, 482, 2039
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  * Camps & Baes (2015) Camps P., Baes M., 2015, Astronomy and Computing, 9, 20
  * Casey et al. (2014) Casey C. M., Narayanan D., Cooray A., 2014, Physics Reports, 541, 45–161
  * Casey et al. (2018) Casey C. M., et al., 2018, ApJ, 862, 77
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Charlot & Fall (2000) Charlot S., Fall S. M., 2000, ApJ, 539, 718
  * Cole et al. (2000) Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168
  * Cook et al. (2020) Cook R. H. W., Cortese L., Catinella B., Robotham A., 2020, MNRAS, 493, 5596
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937–1961
  * Croton (2013) Croton D. J., 2013, Publ. Astron. Soc. Australia, 30
  * Croton et al. (2016) Croton D. J., et al., 2016, ApJS, 222, 22
  * Curtis-Lake et al. (2021) Curtis-Lake E., Chevallard J., Charlot S., Sandles L., 2021, MNRAS, 503, 4855
  * Dalla Vecchia & Schaye (2012) Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
  * Davies et al. (2016) Davies L. J. M., et al., 2016, MNRAS, 461, 458
  * Davies et al. (2018) Davies L. J. M., et al., 2018, MNRAS, 480, 768
  * Davies et al. (2019) Davies L. J. M., et al., 2019, MNRAS, 483, 1881
  * Davies et al. (2021) Davies L. J. M., et al., 2021, MNRAS, 506, 256
  * Davies et al. (2022) Davies L. J. M., et al., 2022, MNRAS, 509, 4392
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  * De Vis et al. (2019) De Vis P., et al., 2019, Astronomy & Astrophysics, 623, A5
  * Draine (2003) Draine B. T., 2003, ARA&A, 41, 241
  * Driver & Robotham (2010) Driver S. P., Robotham A. S. G., 2010, MNRAS, 407, 2131
  * Elahi et al. (2018) Elahi P. J., Welker C., Power C., Lagos C. d. P., Robotham A. S. G., Cañas R., Poulton R., 2018, MNRAS, 475, 5338
  * Elahi et al. (2019a) Elahi P. J., Cañas R., Poulton R. J. J., Tobar R. J., Willis J. S., Lagos C. d. P., Power C., Robotham A. S. G., 2019a, Publ. Astron. Soc. Australia, 36
  * Elahi et al. (2019b) Elahi P. J., Poulton R. J. J., Tobar R. J., Cañas R., Lagos C. d. P., Power C., Robotham A. S. G., 2019b, Publ. Astron. Soc. Australia, 36, e028
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Australia, 34, e058
  * Finkelstein et al. (2015) Finkelstein S. L., et al., 2015, ApJ, 810, 71
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fudamoto et al. (2021) Fudamoto Y., et al., 2021, Nature, 597, 489–492
  * Furlong et al. (2015) Furlong M., et al., 2015, MNRAS, 450, 4486
  * Gardner et al. (2006) Gardner J. P., et al., 2006, Space Sci. Rev., 123, 485
  * González et al. (2011) González V., Labbé I., Bouwens R. J., Illingworth G., Franx M., Kriek M., 2011, ApJ, 735, L34
  * Grogin et al. (2011) Grogin N. A., et al., 2011, ApJS, 197, 35
  * Hao et al. (2011) Hao C.-N., Kennicutt R. C., Johnson B. D., Calzetti D., Dale D. A., Moustakas J., 2011, ApJ, 741, 124
  * Katsianis et al. (2017) Katsianis A., et al., 2017, MNRAS, 472, 919
  * Katsianis et al. (2019) Katsianis A., et al., 2019, ApJ, 879, 11
  * Koekemoer et al. (2007) Koekemoer A. M., et al., 2007, ApJS, 172, 196
  * Koushan et al. (2021) Koushan S., et al., 2021, Monthly Notices of the Royal Astronomical Society, 503, 2033
  * Kregel et al. (2002) Kregel M., van der Kruit P. C., de Grijs R., 2002, MNRAS, 334, 646
  * Krumholz et al. (2009) Krumholz M. R., McKee C. F., Tumlinson J., 2009, ApJ, 699, 850–856
  * Labbé et al. (2013) Labbé I., et al., 2013, ApJ, 777, L19
  * Lagos et al. (2011) Lagos C. D. P., Lacey C. G., Baugh C. M., Bower R. G., Benson A. J., 2011, MNRAS, 416, 1566
  * Lagos et al. (2014) Lagos C. D. P., Baugh C. M., Zwaan M. A., Lacey C. G., Gonzalez-Perez V., Power C., Swinbank A. M., van Kampen E., 2014, MNRAS, 440, 920
  * Lagos et al. (2018) Lagos C. d. P., Tobar R. J., Robotham A. S. G., Obreschkow D., Mitchell P. D., Power C., Elahi P. J., 2018, MNRAS, 481, 3573
  * Lagos et al. (2019) Lagos C. d. P., et al., 2019, MNRAS, 489, 4196
  * Lee et al. (2012) Lee K.-S., et al., 2012, ApJ, 752, 66
  * Lee et al. (2015) Lee N., et al., 2015, ApJ, 801, 80
  * Leslie et al. (2020) Leslie S. K., et al., 2020, ApJ, 899, 58
  * Lilly et al. (2013) Lilly S. J., Carollo C. M., Pipino A., Renzini A., Peng Y., 2013, ApJ, 772, 119
  * Lovell et al. (2021) Lovell C. C., Vijayan A. P., Thomas P. A., Wilkins S. M., Barnes D. J., Irodotou D., Roper W., 2021, MNRAS, 500, 2127
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Matthee & Schaye (2019) Matthee J., Schaye J., 2019, MNRAS, 484, 915
  * McAlpine et al. (2016) McAlpine S., et al., 2016, Astronomy and Computing, 15, 72
  * McLure et al. (2010) McLure R. J., Dunlop J. S., Cirasuolo M., Koekemoer A. M., Sabbi E., Stark D. P., Targett T. A., Ellis R. S., 2010, MNRAS, 403, 960
  * Meurer et al. (1999) Meurer G. R., Heckman T. M., Calzetti D., 1999, ApJ, 521, 64
  * Michałowski et al. (2014) Michałowski M. J., Hayward C. C., Dunlop J. S., Bruce V. A., Cirasuolo M., Cullen F., Hernquist L., 2014, A&A, 571, A75
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493–508
  * Noeske et al. (2007a) Noeske K. G., et al., 2007a, ApJ, 660, L43
  * Noeske et al. (2007b) Noeske K. G., et al., 2007b, ApJ, 660, L47
  * Oesch et al. (2018) Oesch P. A., Bouwens R. J., Illingworth G. D., Labbé I., Stefanon M., 2018, ApJ, 855, 105
  * Pearson et al. (2018) Pearson W. J., et al., 2018, A&A, 615, A146
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A1
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Rémy-Ruyer et al. (2014) Rémy-Ruyer A., et al., 2014, A&A, 563, A31
  * Rigby et al. (2022) Rigby J., et al., 2022, Characterization of JWST science performance from commissioning, http://arxiv.org/abs/2207.05632
  * Robotham et al. (2020) Robotham A. S. G., Bellstedt S., Lagos C. d. P., Thorne J. E., Davies L. J., Driver S. P., Bravo M., 2020, MNRAS, 495, 905
  * Salim et al. (2007) Salim S., et al., 2007, ApJS, 173, 267
  * Schaller et al. (2015) Schaller M., Dalla Vecchia C., Schaye J., Bower R. G., Theuns T., Crain R. A., Furlong M., McCarthy I. G., 2015, MNRAS, 454, 2277
  * Schaye et al. (2010) Schaye J., et al., 2010, MNRAS, 402, 1536
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schenker et al. (2013) Schenker M. A., et al., 2013, ApJ, 768, 196
  * Shen et al. (2021) Shen X., Vogelsberger M., Nelson D., Tacchella S., Hernquist L., Springel V., Marinacci F., Torrey P., 2021, arXiv e-prints, p. arXiv:2104.12788
  * Shivaei et al. (2020) Shivaei I., et al., 2020, ApJ, 899, 117
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Song et al. (2016) Song M., et al., 2016, ApJ, 825, 5
  * Speagle et al. (2014) Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105–1134
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Stanway & Eldridge (2018) Stanway E. R., Eldridge J. J., 2018, Monthly Notices of the Royal Astronomical Society, 479, 75
  * Stefanon et al. (2021) Stefanon M., Bouwens R. J., Labbé I., Illingworth G. D., Gonzalez V., Oesch P. A., 2021, ApJ, 922, 29
  * Tacchella et al. (2016) Tacchella S., Dekel A., Carollo C. M., Ceverino D., DeGraf C., Lapiner S., Mandelker N., Primack Joel R., 2016, MNRAS, 457, 2790
  * Terrazas et al. (2020) Terrazas B. A., et al., 2020, MNRAS, 493, 1888
  * Thorne et al. (2021) Thorne J. E., et al., 2021, MNRAS, 505, 540
  * Trayford et al. (2020) Trayford J. W., Lagos C. d. P., Robotham A. S. G., Obreschkow D., 2020, MNRAS, 491, 3937
  * Vazdekis et al. (2016) Vazdekis A., Koleva M., Ricciardelli E., Röck B., Falcón-Barroso J., 2016, MNRAS, 463, 3409
  * Vijayan et al. (2019) Vijayan A. P., Clay S. J., Thomas P. A., Yates R. M., Wilkins S. M., Henriques B. M., 2019, MNRAS, 489, 4072
  * Vijayan et al. (2021) Vijayan A. P., Lovell C. C., Wilkins S. M., Thomas P. A., Barnes D. J., Irodotou D., Kuusisto J., Roper W. J., 2021, MNRAS, 501, 3289
  * Vijayan et al. (2022) Vijayan A. P., et al., 2022, MNRAS, 511, 4999
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Vogelsberger et al. (2020) Vogelsberger M., et al., 2020, MNRAS, 492, 5167
  * Weigel et al. (2016) Weigel A. K., Schawinski K., Bruderer C., 2016, MNRAS, 459, 2150
  * Whitaker et al. (2012a) Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012a, ApJ, 754, L29
  * Whitaker et al. (2012b) Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012b, ApJ, 754, L29
  * Wilkins et al. (2022a) Wilkins S. M., et al., 2022a, First Light And Reionisation Epoch Simulations (FLARES) V: The redshift frontier, http://arxiv.org/abs/2204.09431
  * Wilkins et al. (2022c) Wilkins S. M., et al., 2022c, First Light And Reionisation Epoch Simulations (FLARES) VI: The colour evolution of galaxies $z=5-15$, http://arxiv.org/abs/2207.10920
  * Wilkins et al. (2022b) Wilkins S. M., et al., 2022b, First Light And Reionisation Epoch Simulations (FLARES) VII: The Star Formation and Metal Enrichment Histories of Galaxies in the early Universe, http://arxiv.org/abs/2208.00976
  * Wuyts et al. (2011) Wuyts S., et al., 2011, ApJ, 742, 96
  * Yabe et al. (2009) Yabe K., Ohta K., Iwata I., Sawicki M., Tamura N., Akiyama M., Aoki K., 2009, ApJ, 693, 507

##  Appendix A Effect of uncertainties on the SFS fitting

Figure 13: Relative error on each parameter, as labelled, used for fitting the
SFS for the JWST population as a function of equal Gaussian uncertainty upon
the stellar masses and star formation rates. Orange colours show the results
for Flares and blue colours show the results for Shark. Each redshift,
${z=6,8,10}$, is shown by a shading gradient with ${z=10}$ being the darkest
shade. The black line highlights no relative error, and the shaded region
shows the 10% level.

To investigate the effect of measurement uncertainties on the ability to
recover the SFS parameters of the true galaxy population we refit the SFS of
the JWST population with varied Gaussian uncertainties added to the stellar
masses and star formation rates. Figure 13 shows the relative error of the
parameters of the SFS as a function of the assumed uncertainty on the stellar
masses and star formation rates of the JWST detected populations in the
simulations. The relative error is

$\mathrm{Relative\,error=\left|\frac{X_{JWST}-X_{All}}{X_{All}}\right|,}$ (11)

where $\mathrm{X_{JWST}}$ are the parameters calculated for the JWST
population with the absorbed uncertainties and $\mathrm{X_{All}}$ are the
calculated parameters for the entire galaxy population in the simulations.
When using the same prior distributions as those we used for the fitting in
Figures 5 and 6 we found that for varied uncertainty combinations the relative
error was largest in different parameters, making it difficult to see just how
additive normal noise affected the fitting. For this reason, we fit the
parameters independently, keeping the remaining three parameters fixed, with
each uncertainty combination. We achieve this by using a uniform prior on the
desired parameter and narrow normal priors centred on the true population
parameters for the remaining three parameters, cycling through this process
until all parameters have been considered. So for 7 uncertainty combinations,
4 parameters and 3 selected redshifts we perform $\mathrm{7\times 4\times
3=84}$ MCMC fitting routines on each simulation.

Both simulations predict a similar trend of increased relative error on all
parameters of the SFS as the uncertainty on the stellar masses star formation
rates increases, as expected. Even with zero uncertainty there is a non-zero
relative error indicating systematic offsets between the JWST detected and
total galaxy populations in the simulations. This is perhaps unsurprising as
the JWST will miss faint galaxies, and so the relationship between the total
and JWST populations is not one-to-one. There is, in general, a larger
increase in the relative error when the uncertainty on the stellar mass is
increased from 0.2 dex to 0.4 dex compared to when the uncertainty on the SFR
is increased from 0.2 dex to 0.4 dex to 0.6 dex. We thus suggest that accurate
constraints on the stellar mass are more necessary to describe the SFS than
constraints on SFR. The grey shaded regions shows the 10% margin. Beyond
$\mathrm{\sim(0.2,0.4)}$ dex uncertainties on $\mathrm{(M_{\star},SFR)}$ both
simulations fail to recover the normalisation of the SFS to within 10% at all
redshifts shown. For almost any uncertainty combination, the simulations are
able to recover the turn over mass to within $10$% at all redshifts except
${z\sim 6}$ in Shark. The intermediate and high stellar mass slopes require
far more accurate constraints on the stellar mass and SFR. Beyond
$\mathrm{\sim 0.2\,dex}$ and $\mathrm{\lesssim 0.6\,dex}$ uncertainty the
intermediate stellar mass slope cannot be recovered to within 10% at all
redshifts except ${z=6}$ in Flares. For most uncertainty combinations used in
Flares the high stellar mass slope can be recovered except beyond
$\mathrm{\sim(0.2,0.6)}$ dex uncertainties on $\mathrm{(M_{\star},SFR)}$. On
the other hand, for almost every uncertainty combination Shark cannot recover
the high stellar mass slope at all redshifts suggesting that the high stellar
mass slope is subtle in Shark by default. Interestingly, the apparent redshift
dependencies on the recovery of these parameters appear to be reversed between
the simulations, with Flares showing decreasing relative error for the same
uncertainty combination between ${z=10\to 6}$, whereas this is true between
${z=6\to 10}$ in Shark. We do not comment any further on this in this work.

Some larger combinations of uncertainties result in a lower relative error
than smaller uncertainties. For example, the $\mathrm{\sim(0.4,0.6)}$ dex
combination is closer to the true population value of the turn over at all
redshifts in Shark than the $\mathrm{\sim(0.4,0.2)}$ combination is. This is
more than likely a statistical effect or a symptom of the fitting routine as,
intuitively, a larger uncertainty should result in worse agreement; though we
do not investigate this further here. Ultimately, these trends show that
poorer constraints on $\mathrm{M_{\star}}$ and $\mathrm{SFR}$ result in a
worse recovery of the true parameters of the SFS. These trends also provide an
indication on the level of accuracy that would be needed to adequately
determine the true SFS to within 10%.