First Light And Reionisation Epoch Simulations (FLARES) V: The redshift frontier

The James Webb Space Telescope (JWST) is set to transform many areas of astronomy, one of the most exciting is the expansion of the redshift frontier to $z>10$. In its first year alone JWST should discover hundreds of galaxies, dwarfing the handful currently known. To prepare for these powerful observational constraints, we use the First Light And Reionisation Epoch (FLARES) simulations to predict the physical and observational properties of the $z>10$ population of galaxies accessible to JWST. This is the first time such predictions have been made using a hydrodynamical model validated at low redshift. Our predictions at $z=10$ are broadly in agreement with current observational constraints on the far-UV luminosity function and UV continuum slope $\beta$, though the observational uncertainties are large. We note tension with recent constraints $z\sim 13$ from Harikane et al. 2022 - compared to these constraints, FLARES predicts objects with the same space density should have an order of magnitude lower luminosity, though this is mitigated slightly if dust attenuation is negligible in these systems. Our predictions suggest that in JWST's first cycle alone, around $600$ galaxies should be identified at $z>10$, with the first small samples available at $z>13$.


INTRODUCTION
Identifying and characterising the first generation of galaxies is one of the core aims of modern extragalactic astronomy. Doing so will provide the essential constraints to galaxy formation models, helping us elucidate the key physics of early galaxy formation and evolution (Dayal & Ferrara 2018;Robertson 2021).
Over the last 15 years, remarkable progress has been made in studying the distant Universe, with > 1000 candidate galaxies now identified at z = 6 − 10 (e.g. Bouwens et al. 2021). These have come predominantly from the analysis of deep near-IR observations from Wide Field Camera 3 (WFC3) on the Hubble Space Telescope (e.g. Bunker et al. 2010;Bouwens et al. 2010;McLure et al. 2010;Wilkins et al. 2011;Finkelstein et al. 2015;Bouwens et al. 2021), with a small but growing sample of bright sources, more amenable to multi-wavelength follow-up, identified from ground-based imaging (e.g. Bowler et al. 2015). A number of these candidates have now been confirmed by spectroscopy, targeting Ly α (e.g Stark et al. 2010;Curtis-Lake et al. 2012;Caruana et al. 2014; E-mail: s.wilkins@sussex.ac.uk 2014; Schenker et al. 2014;Stark et al. 2017;Mason et al. 2019) or far-IR lines (e.g Knudsen et al. 2016;Pentericci et al. 2016;Hashimoto et al. 2019).
While the vast majority of sources are at lower-redshift, a handful of objects have now been detected at z > 10, often combining Hubble and Spitzer observations. These include the surprisingly bright galaxy GN-z11 (Oesch et al. 2016) and more recently a pair of candidates at z ∼ 13 (Harikane et al. 2021). With the successful launch of the James Webb Space Telescope (JWST), these sources look set to be only the first of many identified at z > 10 (Robertson 2021). JWST offers the sensitivity, survey efficiency, and wavelength coverage to push well beyond the current redshift frontier. In addition, JWST will provide rest-frame UV spectroscopy, allowing the confirmation of sources and the accurate measurement of many key physical properties. With the first results from JWST imminent, it is essential that we have theoretical predictions in place to allow us to interpret these revolutionary observations.
The main drawback of hydrodynamical simulations -especially those including RT -is that they are computationally expensive, limiting their volume, resolution, and/or redshift end point. This is a particular problem at very high redshift, where the source density is so low that simulations comparable to, and ideally much larger, observational surveys are essential to yield useful statistical predictions. Indeed, flagship cosmological simulations that have been validated at low-redshift, like EAGLE (Schaye et al. 2015;Crain et al. 2015), SIMBA , ILLUSTRIS (Vogelsberger et al. 2014a,b;Genel et al. 2014;Sijacki et al. 2015), and ILLUSTRIS-TNG (Naiman et al. 2018;Nelson et al. 2018;Marinacci et al. 2018;Springel et al. 2018;Pillepich et al. 2018) fail to pass this threshold, with only a small number of observationally accessible sources at z > 10.
One solution is to carry out much larger simulations, but limited to high-redshift. This is a strategy implemented by, e.g., MASSIVEBLACK (Khandai et al. 2012), BLUETIDES (Feng et al. 2016;Wilkins et al. 2017), and ASTRID (Ni et al. 2022). An alternative strategy is to re-simulate a range of environments drawn from a very large low-resolution parent simulation (e.g. Crain et al. 2009). This has the advantage of more efficiently allowing us to extend the dynamic range (see discussion in Lovell et al. 2021). These individual re-simulations are also much smaller than single large boxes, allowing wider and more efficient use of HPC systems. The chief disadvantages are the loss of most clustering information and the requirement to carefully understand the weightings of the individual simulations to obtain the correct statistical properties of the galaxy population. Machine learning approaches may allow us to overcome some of these limitations (e.g. Lovell et al. 2022;Bernardini et al. 2022).
Re-simulations of multiple regions is utilised in the FLARES: First Light And Reionisation Epoch Simulations project. FLARES combines the validated at z = 0 EAGLE physics model with a resimulation strategy yielding a much larger effective volume and dynamic range. Compared to the (100 Mpc) 3 EAGLE reference simulation, FLARES contains 10-100× as many high-redshift galaxies. In this article we use FLARES to make predictions for the galaxy population at z > 10, building on earlier work focused at z = 5 − 10 Vijayan et al. 2021Vijayan et al. , 2022Roper et al. 2022).
This article is organised as follows: in Section 2 we briefly describe the FLARES project. In Section 3 we explore the physical properties of galaxies at z > 10 and in Section 4 explore their observational properties, including the UV luminosity function and forecasts for upcoming surveys ( §4.1), the UV continuum slope β ( §4.3), UV emission lines ( §4.4), and the UV -optical colours ( §4.5). Finally, in Section 5 we present our conclusions.

THE FIRST LIGHT AND REIONISATION EPOCH SIMULATIONS
In this study, we make use of the First Light And Reionisation Epoch Simulations (FLARES; Lovell et al. 2021;Vijayan et al. 2021). FLARES is a suite of 40 spherical re-simulations, 14 h −1 Mpc in radius, of regions selected from a large (3.2 Gpc) 3 dark matter only simulation. The regions selected to re-simulate span a range of en-vironments: (at z ≈ 4.7) log 10 (1 + δ 14 ) = [−0.3, 0.3] 1 with overrepresentation of the extremes of the density contrast distribution. FLARES adopts the AGNdT9 variant of the EAGLE simulation project (Schaye et al. 2015;Crain et al. 2015). The AGNdT9 variant implements a higher heating temperature from active galactic nuclei (AGN) compared to the reference EAGLE run, thus producing more energetic, less frequent feedback events. We adopt identical resolution to the fiducial EAGLE simulation, i.e. dark matter and initial gas particle masses of m dm = 9.7 × 10 6 M and m g = 1.8 × 10 6 M , respectively, and a softening length of 2.66 ckpc. As with the original EAGLE simulations, we assume Ω m = 0.307, Ω Λ = 0.693, h = 0.6777 based on results from Planck Collaboration et al. (2014).
Galaxies in FLARES are first identified as groups via the Friends-Of-Friends (FOF, Davis et al. 1985) algorithm, and subsequently subdivided into bound groups with the SUBFIND (Springel et al. 2001;Dolag et al. 2009) algorithm. When measuring properties we use 30 kpc radius apertures, centred on the most bound particle of each subgroup.

Spectral Energy Distribution Modelling
Key to making observational predictions is the spectral energy distribution (SED) modelling applied to the simulation outputs. The SED modelling approach is presented in Vijayan et al. (2021), broadly following the approach developed by Wilkins et al. (2013bWilkins et al. ( , 2016cWilkins et al. ( , 2018Wilkins et al. ( , 2020, with modifications to the dust treatment. In short, we begin by associating each star particle with pure stellar spectral energy distribution (SED) using v2.2.1 of the Binary Population and Spectral Synthesis (BPASS; Stanway & Eldridge 2018) stellar population synthesis model assuming a Chabrier (2003) initial mass function (IMF). We then associate each star particle with HII region giving rise to nebular continuum and line emission. Specifically, we follow the approach detailed in Wilkins et al. (2020), in which the pure stellar spectrum is processed with the cloudy photo-ionisation code (Ferland et al. 2017). We account for the effects of dust attenuation in both the birth clouds of young stellar populations and the interstellar medium (ISM). The latter is accounted for using a line-of-sight (LOS) model similar to that described in Wilkins et al. (2018). In this model we treat stellar particles as emitters along a line of sight and account for the attenuation due to gas particles which intersect this LOS. To do this we determine the metal column density and convert this to a dust optical depth using the fitting function for the dust-to-metal ratio presented in Vijayan et al. (2019). For the attenuation due to the birth cloud component, we scale it with the star particle metallicity, thus assuming a constant dust-to-metal ratio. The proportionality factors for the two components are fixed to match the z = 5 UV LF from Bouwens et al. (2015), z = 5 UV-continuum slope relations and the [OIII]λ 4959, 5007 + Hβ line luminosity and EW relations at z = 8 from De Barros et al. (2019a).

Comparison to EAGLE
The core objective of FLARES is to expand both number and dynamic range of simulated galaxies at z > 5 compared to the EAGLE reference simulation. The number of galaxies in both FLARES and EAGLE with stellar mass greater than ∈ {10 8 , 10 9 , 10 10 , 10 11 } M at z > 5 are shown in Fig. 1. At z = 10 FLARES contains 100× (10×) as many galaxies as EAGLE with M > 10 9 M (M > 10 8 M ). At z = 10 FLARES contains ∼ 1000 galaxies at M > 10 8 M , dropping to 10 by z = 15.

Environmental dependence
A key feature of FLARES is the explicit simulation of a wide range of environments with log 10 [1+δ 14 (z = 4.7)] ≈ [−0.3, 0.3]. In Lovell et al. (2021) we showed that the galaxy stellar mass function, and thus the total number of galaxies above a mass threshold, was extremely sensitive to the environment. A consequence of this, and the low number density of galaxies at z > 10, is that the vast majority of our simulated galaxies are in very over-dense regions as shown in Fig. 2. At z = 15 only one simulation with δ < 0.5 contains any galaxies. Since each simulation is appropriately weighted (see Lovell et al. 2021) the lack of any galaxies in many density contrast bins should not affect distribution functions like the galaxy stellar mass function or UV luminosity function. However, if galaxy scaling relationships are sensitive to the environment, even the appropriately weighted relations could be biased. Fortunately, Lovell et al. (2021) found no significant evidence of environmental dependence in the key scaling relations.

PHYSICAL PROPERTIES
We begin by exploring a selection of key physical properties of galaxies at z > 10. In Fig. 3 we show the relationship between stellar mass and the specific star formation rate, age, gas-phase metallicity, and far-UV dust attenuation. At present, there are no observational constraints for these properties but this should soon change.
The relationship between stellar mass and specific star formation rate is predicted to be fairly flat, though with significant redshift evolution of the normalisation. Similarly, the average age -defined here as the time since its stellar mass was half its current value -is flat with stellar mass but evolves strongly with redshift. At z = 10, the average age is predicted to be ≈ 20 Myr rising to almost ≈ 100 Myr at z = 10. On the other hand, the gas phase metallicity Z g shows a significant trend with stellar mass but little redshift evolution. Similarly we see a strong trend between the far-UV attenuation and stellar mass but little redshift evolution. Given the gas phase metallicity trends this is unsurprising since the attenuation in FLARES is related to surface density of metals.

Mass-to-light ratio
One of the most fundamental properties is the mapping between stellar mass and the (dust-attenuated) far-UV luminosity. This encodes the star formation and metal enrichment history of each galaxy in addition to reprocessing by dust and gas. This relationship is shown in Fig. 4. At M < 10 8.5 M this relationship is close to linear; at highmasses, however, the luminosity falls below the linear expectation. This is predominantly due to the effects of dust, but is also affected by the higher stellar metallicities in the most massive galaxies. The scatter in this relationship is ≈ 0.2 dex independent of redshift and stellar mass.

OBSERVATIONAL PROPERTIES
We now turn our attention to some of the properties of galaxies that can be observed at z > 10. JWST will, for the first time, provide deep >2 µm imaging and spectroscopy allowing us to measure several properties including the rest-frame UV luminosity (and thus luminosity function), the UV continuum slope, UV emission lines, and potentially even UV -optical colours via MIRI imaging. Model spectral energy distributions of star forming galaxies at z = 10 and z = 15 are shown, alongside the JWST filter transmission functions in Fig. 5.

Far-UV luminosity function
The rest-frame far-UV luminosity function (LF) is one of the key statistical descriptions of the galaxy population at high-redshift. This is predominantly due to its accessibility but also the fact that the observed UV light traces both unobscured star formation (Kennicutt & Evans 2012;Wilkins et al. 2019) and the production of ionising photons (Wilkins et al. 2016b). The far-UV LF has now been measured extensively to z ∼ 8 with tentative constraints at ∼ 10 (e.g McLeod et al. 2016;Oesch et al. 2018;Finkelstein et al. 2021) and more recently at z ∼ 13 (Harikane et al. 2021).
In Fig. 6 we show the far-UV luminosity predicted by FLARES at z = 15 → 10 alongside both these observational constraints and other model predictions. We show both the observed (dust-attenuated) and intrinsic distribution functions. The far-UV LF follows the familiar steep decline with luminosity seen at low-redshift, dropping by roughly a factor of 10 3 from L = 10 28 → 10 29 erg s −1 Hz −1 . The number density of sources also evolves strongly, increasing by ∼ 10× from z = 14 → 10 with stronger evolution at the bright end. At low-luminosity (L < 10 28.5 erg s −1 Hz −1 ) the impact of dust is small leaving the intrinsic and observed LF similar. However, as noted previously brighter, more massive galaxies, increasingly have stronger dust attenuation leading to a divergence in the predictions. Shown on the z = 13 panel are the recent observational constraints from Harikane et al. (2021). This study used observations from Hyper Suprime-Cam, VISTA, and Spitzer of the COSMOS and SXDS fields to identify a pair of bright sources with spectral energy distribution consistent with z ∼ 13 star forming galaxies. In addition, one source has a tentative line detection consistent with [OIII]88µm at z = 13.27, in-line with its photometric redshift. If real these sources suggest little evolution in the bright end of the UV luminosity function from z ∼ 10 → 13. FLARES contains objects with a similar space density as these sources but with dust -attenuated luminosities around a factor of 10× smaller. As noted, at z > 10 the FLARES dust model is likely to become increasingly unreliable since it is based on modelling calibrated at lower-redshift. However, even using in-  trinsic luminosities, sources with a similar space density in FLARES are still around 5× fainter, suggesting significant remaining tension. Fig. 6 also shows a comparison with other model predictions at z ≥ 10 including the semi-empirical model of Mason et al. (2015), the Santa Cruz semi-analytical model (Yung et al. 2019), the large volume cosmological hydrodynamical simulation BLUETIDES (Wilkins et al. 2017), and the high-resolution FIRE-2 simulations (Ma et al. 2019). While there is good agreement at z = 10, the agreement with Mason et al. (2015) and BLUETIDES breaks down at highredshift; FLARES predicts a similar density of bright galaxies but consistently predicts more faint galaxies.

JWST Cycle 1 forecasts
As noted in the introduction the z > 10 galaxy population will soon be accessible via several deep imaging surveys conducted by JWST. Using our luminosity function predictions we can now forecast the number of sources expected for each survey and the eventual constraints on the z = 10 − 15 LF.
We begin, in Fig. 7, by presenting forecasts for the cumulative number of sources accessible to various JWST Cycle 1 GO, ERS, and GTO programmes. These include: CEERS, NG-DEEP 2 , PRIMER, COSMOS-Web 3 , JADES, the Northern Ecliptic Pole element of PEARLS 4 , and PANORAMIC. These predictions assume 100% completeness down to the 10σ point-source F277W depth, which in reality is likely to be difficult to achieve. The approximate depths/areas for each survey were provided by each programme PI except for JADES which are taken from Williams et al. (2018). Across all 7 programmes we predict ∼ 500, 85, and 6 galaxies at z > 10, z > 12, and z > 14 respectively. JADES is predicted to dominate these numbers contributing around half of expected sources at these redshifts.
In Fig. 8 we then make forecasts for the z = 15 → 10 luminosity function for the combination of the Cycle 1 programmes. The result is strong constraints at z = 10, comparable to the current z = 7 con- straints from Hubble's entire campaign (Bouwens et al. 2021). While these constraints progressively weaken toward higher redshift subsequent observations throughout JWST's tenure should ultimately enable useful constraints to z ∼ 15 and potentially beyond. Crucially this will allow us to differentiate between the different model predictions in this era.

Sizes
In Fig. 9 we present measurements of galaxy half light radius in the far-UV for z = 10 − 12. The intrinsic size measurements are derived from the particle distribution, while the observed size measurements use a non-parametric pixel approach in which size is derived from the non-contiguous pixel area containing half the galaxy's total luminosity. The latter approach well encompasses the clumpy natures of high redshift galaxies (e.g Jiang et al. 2013;Bowler et al. 2017). The redshift range is limited by the number of galaxies in FLARES with a sufficient number of stellar particles to make robust size measurements (N ≥ 100). As in Roper et al. (2022) we impose a 95 percent completeness limit in each redshift bin. Intrinsically the high redshift galaxy population is extremely compact with a negative size-luminosity relation. However, the bright central regions of these compact galaxies are also efficiently seeding their surroundings with metals, even at this early epoch. This seeding leads to strong dust obscuration in the brightest regions, and thus a large increase in the observed size and a positive observed size-luminosity relation. The normalisation of the observed size-luminosity relation quickly evolves at these high redshifts, with an increase of ∼ 0.1 dex from z = 12 to z = 10 driven by increasing obscuration of the brightest regions due to the formation of dust in the highly star forming cores of these bright galaxies. JWST will not only be able to probe this obscured size-luminosity relation in the rest frame far-UV, the reddest NIRCam filters will also be able to probe deeper into the increasingly unobscured size-luminosity relation at longer wavelength. This will provide a valuable view into the intrinsic size-luminosity relation and its negative slope.

The UV continuum slope
The most accessible spectral diagnostic available at high-redshift is the UV continuum slope β : as can be seen in Fig. 5 the rest-frame UV continuum to λ = 350 nm is accessible to NIRCam to z ≈ 15 and can be measured with 3-4 of NIRCam's wide filters. While primarily an indicator of dust attenuation the UV continuum slope is also sensitive to the star formation and metal enrichment history, the Lyman continuum escape fraction, the initial mass function, and our understanding of stellar evolution and atmospheres 5 (e.g Wilkins et al. 2013b,a). Fig. 10 shows predictions for β from z = 15 → 10 colour coded by the rest-frame far-UV attenuation. At L FUV < 10 28.5 erg s −1 Hz −1 (M FUV > −19.5) the slope is ≈ −2.4 with little evolution with redshift. At L FUV > 10 28.5 erg s −1 Hz −1 the slope progressively reddens due to the increasing dust attenuation. At L FUV > 10 29 erg s −1 Hz −1 the average slope has reddened to ≈ −2. Fig.  10 also shows observational constraints from Wilkins et al. (2016a); while the uncertainties are large these observations are consistent with our predictions.
The origin of the UV continuum slope is explored in more detail in Fig. 11. Here we show the median β for unprocessed starlight (dotted line), starlight with reprocessing by gas (dashed line), and reprocessing with gas and dust (solid line, the same as that shown in Fig. 10). With no reprocessing the predicted slopes are ≈ −2.8 at low-luminosity rising to ≈ −2.5 at L FUV ≈ 10 29.5 erg s −1 Hz −1 with some weak (∆β = 0.1) evolution with redshift z = 15 → 10. These trends reflect variation in the star formation and metal enrichment history with the brightest galaxies typically having higher metallicities. The addition of nebular (continuum) emission reddens β . As the impact of nebular emission is strongest at low metallicity this has the effect of flattening the previous trend with luminosity and redshift leaving galaxies with β ≈ 2.5. The addition of dust has an impact at all luminosities though this is most pronounced at L FUV > 10 29 erg s −1 Hz −1 where dust is predicted to redden the slope by ≈ 0.5.
As noted previously the impact of dust attenuation in FLARES is particularly uncertain at z > 10. Recent luminosity function constraints (i.e. Finkelstein et al. 2021;Harikane et al. 2021) tentatively suggest the presence of too much dust attenuation in FLARES. Precise constraints on the UV continuum slope from JWST -and ideally ALMA dust continuum observations -in bright (L FUV > 10 29 erg s −1 Hz −1 ) should allow us to determine whether this is the case.

UV emission lines
In addition to deep near-IR imaging, JWST will also be able to obtain deep near-IR spectroscopy utilising NIRSpec, NIRCam, and NIRISS 6 . NIRSpec provides both a multi-object and integral field  Table 1. unit spectroscopy at 0.6-5.3µm, while NIRCam and NIRISS together provide wide field slit-less spectroscopy across the near-IR. At z < 10 this enables the observation of various strong optical lines. However, at z > 10 most of the strong lines fall outside the accessible range leaving a handful of weaker UV lines. Most prominent amongst these is the [CIII],CIII]λ λ 1907, 1909Å doublet for which a handful of detections are already available at z > 6 ( Stark et al. 2015Stark et al. , 2017Topping et al. 2021).
Predictions from FLARES for the rest-frame equivalent width distribution of [CIII],CIII] are presented in Fig. 12 as a function of the observed (dust attenuated) far-UV luminosity. Equivalent widths show a weak decline to higher luminosity and lower-redshift. The median value is ≈ 10Å with the tail of the distribution reaching beyond 20Å. Unsurprisingly, given the younger age, these predictions are offset to higher equivalent widths than observations at lower redshift (e.g Maseda et al. 2017;Llerena et al. 2022). Constraints at z > 6 include a handful of detections and upper-limits and are likely biased due to their selection method. However, two of the detections: EGS-zs8-1 (Stark et al. 2017) and AEGIS-33376 (Topping et al. 2021) have equivalent widths at the upper extreme of the predicted distribution suggesting some possible tension. This may reflect some of the simplifying assumptions used in our modelling. For example, Wilkins et al. (2020) showed that the equivalent width of [CIII],CIII] is strongly sensitive to the ionisation parameter and hydrogen density, for which we adopted single fiducial values.

UV -Optical colours
The near-UV -optical colour is another key spectral diagnostic of galaxies, its measurement providing insights into the star formation and metal enrichment history, dust attenuation, and nebular emission of galaxies. In the context of z > 6 galaxies the near-UV -optical colour is chiefly impacted by nebular line emission and can be used to infer the strength of combination of the Hβ and [OIII] lines (e.g De Barros et al. 2019b;Endsley et al. 2021). Where a spectroscopic redshift is available it is sometimes possible to avoid strong line emission providing "clean" constraints on the strength of the Balmer break (e.g Hashimoto et al. 2018) and thus a more accurate constraint on the star formation and metal enrichment history. For a wider introduction to the break in the context of the high-redshift Universe see Wilkins et al. in-prep. In principle JWST can observe the rest-frame optical to z = 15 and beyond. However, at z > 11 the optical falls beyond the range accessible to JWST's near-IR instruments and thus requires MIRI observations. Several blind-field cycle 1 surveys will simultaneously collect both multi-band NIRCam and MIRI imaging. For MIRI F770W is the most popular choice and with this in mind we present the predicted F444W-F770W colour in Fig. 13. We do this for both the pure stellar photometry and the photometry including nebular emission and dust. In both cases there is little trend with the UV luminosity. The addition of nebular emission (and dust) has the impact of significantly reddening the colour by ≈ 0.4 mag at z = 10 increasing to ≈ 0.7 mag at z = 15. This reflects both the increasing optical line equivalent widths but also the lines that fall within the F770W band. While in principle, MIRI observations, when combined with NIR-Cam, should thus allow us to constrain optical line emission at z > 10, unfortunately MIRI has both a lower sensitivity 7 and smaller field-of-view than NIRCam resulting in a much reduced survey efficiency. Blind field cycle 1 programmes (i.e. PRIMER, COSMOS-Web, CEERS) obtaining both NIRCam and MIRI F770W observations typically reach depths F770W 3 magnitudes shallower than the deepest NIRCam observations and only over less than half the total area. Combined with the expected surface densities and flux distributions it is then unlikely, even with stacking, that MIRI will yield useful constraints at z > 10. However, there is the possibility of obtaining deep MIRI imaging of individual bright targets in later cycles.

CONCLUSIONS
In this article we have presented theoretical predictions for the properties of galaxies at z = 10 − 15 from the FLARES: First Light And Reionisation Epoch Simulations. These are amongst the first predictions from a cosmological hydrodynamical simulation calibrated at z = 0 at these redshifts enabled by the unique simulation strategy adopted by FLARES. Our major findings are: • Specific star formation rates and ages show little trend with stellar mass though evolve strongly with redshift. However, gas-phase metallicities and dust attenuation rapidly increase with stellar mass but show little redshift evolution.
• The far-UV luminosity function continues its evolution from lower-redshift with the luminosity density predicted to drop by ∼ 10× from z = 10 → 14. At z = 10 the predictions are consistent with observational constraints from McLeod et al. (2016), Oesch et al. (2018), and Finkelstein et al. (2021) though favour less dust attenuation. FLARES contains galaxies with a similar space density to those recently identified by Harikane et al. (2021) but ∼ 5 − 10× fainter depending on whether dust attenuation is included. Similarly agreement with other models is good at z = 10 but diverges to higherredshift with FLARES predicting more faint galaxies than other models. Based on these predictions in cycle 1 alone JWST should identify ∼600, 100, and 6 galaxies at z > 10, 12, and 14 respectively providing robust constraints on the LF to z ∼ 13.
• FLARES predicts little redshift evolution in the relationship between the UV continuum slope β and UV luminosity. The brightest galaxies are predicted to be moderately reddened (∆β ≈ 0.5) by dust.
• UV-optical colours probed by NIRCam and MIRI are likely to be dominated by nebular emission though will be hard to measure due to MIRI's much lower sensitivity.

ACKNOWLEDGEMENTS
We dedicate this article to healthcare and other essential workers, the teams involved in developing the vaccines, and to all the parents who found themselves having to home-school children while holding down full-time jobs. We thank the EAGLE team for their efforts in developing the EAGLE simulation code.  (Hunter 2007). This research made use of ASTROPY http://www.astropy.org a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013. Parts of the results in this work make use of the colormaps in the CMASHER package (van der Velden 2020).

DATA AVAILABILITY STATEMENT
Binned data for making easy comparisons is available in ascii formats at https://github.com/stephenmwilkins/flares_ frontier_data. Data from the wider FLARES project is available at https://flaresimulations.github.io/data.html. If you use data from this paper please also cite Lovell et al. (2021) and Vijayan et al. (2021). . The redshift evolution of the size -luminosity relation predicted by FLARES in the redshift range z = 10 − 12. The lines represent the weighted 50th percentile of the galaxy distribution weighted using the FLARES weighting scheme. Solid lines show the observed size including the effects of dust attenuation, while dashed lines show the intrinsic stellar emission. Observed sizes are measured using a non-parametric pixel approach while intrinsic sizes use a particle based approach.