Jacopo Chevallard
email: chevallard@iap.fr

Emma Curtis-Lake

Stephane Charlot

Pierre Ferruit

Giovanna Giardino

Marijn Franx

Michael V Maseda

Ricardo Amorin

Santiago Arribas

Andy Bunker

Stéphane Charlot

Stefano Carniani

Bernd Husemann

Peter Jakobsen

Roberto Maiolino

Janine Pforr

Timothy D Rawle

Hans-Walter Rix

Renske Smit

Chris J Willott

Simulating and interpreting deep observations in the Hubble Ultra Deep Field

Keywords: galaxies: evolution -galaxies: formation -galaxies: ISM -HII regions -dark ages, reionization, first stars -telescopes

INTRODUCTION

Our picture of the formation and evolution of galaxies across cosmic time has improved significantly in the past 15 years.

Galaxy surveys have provided researchers with a wealth of spectrophotometric data, at both low (e.g. SDSS, York et al. 2000) and high (e.g., GOODS, Stanway et al. 2003;GLARE, Stanway et al. 2004, 2007;COSMOS, Scoville et al. 2007; HUDF, Bunker et al. 2004;Beckwith et al. 2006;[START_REF] Bouwens | [END_REF]Bunker et al. 2010;[START_REF] Mclure | [END_REF]Ellis et al. 2013;Illingworth et al. 2013;CANDELS, Grogin et al. 2011;Koekemoer et al. 2011;VUDS, Le Fèvre et al. 2015;VANDELS, McLure et al. 2017) redshifts. On the theoretical front, computer simulations of galaxy formation (e.g. EAGLE, Schaye et al. 2015;IllustrisTNG, Pillepich et al. 2017) can now accurately predict several key galaxy properties, such as the redshift evolution of the galaxy stellar mass function and star formation rate density (e.g. Genel et al. 2014;Furlong et al. 2015). Yet, many details of the physical processes driving the evolution of galaxies remain unknown. AGN feedback is a key ingredient of galaxy formation models (see Somerville & Davé 2015 and references therein), but observational evidence for AGN-driven 'quenching' of star formation is ambiguous (e.g. Carniani et al. 2016;Suh et al. 2017). Similarly, although several correlations exist between the physical properties of AGNs and their host galaxies (see Kormendy & Ho 2013 and references therein), a causal connection among these properties indicating a 'co-evolution' of galaxies and AGNs has not been clearly demonstrated. The apparent existence of a tight relation between galaxy masses and star formation rates (i.e., 'main sequence', Brinchmann et al. 2004;Noeske et al. 2007), and perhaps of a more 'fundamental' relation involving gas-phase metallicity as well (Mannucci et al. 2010), can be an indication of self-regulation of star formation within galaxies (e.g. Lilly et al. 2013), but the significance and redshift evolution of these relations remains unclear (e.g. Yabe et al. 2015;Telford et al. 2016). The increasing relative abundance of UV-faint galaxies at redshift z 6 (e.g. Bouwens et al. 2015a;Finkelstein et al. 2015) would support a picture in which low-mass star-forming galaxies provide the bulk of H-ionizing (LyC) photons required for cosmic reionization (e.g. Wilkins et al. 2011;Finkelstein et al. 2012;Robertson et al. 2013;Bouwens et al. 2015b), but this conclusion relies on several assumptions about the ionizing emissivity of galaxies and escape fraction of LyC photons from galaxies.

Advancing our understanding of the above processes, and of many others, that drive the evolution of galaxies from the reionization epoch to the present day, requires measuring physical properties -such as stellar masses, star formation rates, stellar and gas metallicities, stellar ages, gas ionization state, the dynamics of gas and stars -for large samples of galaxies, across a broad range of redshifts and galaxy stellar masses. While multi-band photometric campaigns have collected high quality SEDs for thousands of galaxies spanning a wide range of masses out to z 8 (e.g. HUDF, CANDELS), the availability of spectroscopic observations at z 1 is much more limited. At those redshifts, the strongest optical emission lines (e.g. [O ii]λλ3726,3729, Hβ, [O iii]λλ4959,5007, Hα) are shifted to near infrared wavelengths, where observations with ground-based telescopes are challenging because of regions of low atmospheric transmittance, bright sky background and contamination from bright, variable sky lines. Optical emission lines have been measured out to z ∼ 3 for large samples of galaxies using low-resolution slitless spectroscopy with HST (e.g. WISP, Atek et al. 2010;3D-HST, Brammer et al. 2012;GLASS, Treu et al. 2015;FIGS, Pirzkal et al. 2017). Higher-resolution spectroscopic surveys from the ground (e.g. VUDS, VANDELS) have mainly relied on multi-object spectrographs operating at optical wavelengths (e.g. VIMOS at VLT, Le [START_REF] Fèvre | Proc. SPIEVol. 4841, Instrument Design and Performance for Optical/Infrared Groundbased Telescopes[END_REF], hence are limited, at high redshift, to the measurement of rest-frame UV emission lines. These lines are intrinsically weaker than the optical ones (e.g., see table 5 of Steidel et al. 2016), and their interpretation is complicated by radiative transfer effects, dust attenuation and potential contamination from outflows.1 Multi-object, ground-based spectrographs operating at near infrared wavelengths (e.g. MOSFIRE at Keck, McLean et al. 2008;KMOS at VLT, Sharples et al. 2013) permit the extension of rest-frame optical-emission-line measurements out to z ∼ 4, albeit only for galaxies with relatively bright emission lines (e.g. MOSDEF, Kriek et al. 2015;KBSS, Steidel et al. 2014). These measurements enabled, for example, new constraints on the dust attenuation properties of galaxies at z ∼ 2 (e.g. Reddy et al. 2015) and on the conditions of ionized gas (metal abundances, ionization state) at z ∼ 2 -4 (e.g. Strom et al. 2017;Shapley et al. 2017).

In the near future, the limitations of existing observatories to study the high redshift Universe will largely be overcome by the launch of the James Webb Space Telescope (JWST, Gardner et al. 2006). JWST is expected to revolutionise our view of galaxies at z 4, and of low-mass galaxies at 1 z 4, by providing a wide range of imaging and spectroscopic capabilities in the wavelength range λ ∼ 0.6 -28 µm. In particular, the Near Infrared Spectrograph (NIRSpec, Bagnasco et al. 2007;[START_REF] Birkmann | Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave[END_REF]) onboard JWST will increase by more than an order of magnitude the emission line sensitivity in the wavelength region λ ∼ 1 -2.5 µm covered by existing ground-based instruments, while reaching an even greater sensitivity at wavelengths λ 2.5 µm inaccessible with existing observatories. The multi-object spectroscopic capabilities of NIRSpec will enable the simultaneous measurement of up to ∼ 200 galaxy spectra, allowing the observation of standard optical emission lines for large samples of galaxies out to z 10 (see Section 2.4). Combined with stellar mass measurements based on JWST/NIRCam imaging (Near Infrared Camera, [START_REF] Horner | Proc. SPIEVol. 5487, Optical, Infrared, and Millimeter Space Telescopes[END_REF]), this will provide researchers with unique data to constrain the physical processes driving the evolution of galaxies at z 4, the formation of the first galaxies at z 10, the way cosmic reionization proceeded in space and time and the contribution of different sources to the cosmic reionization budget.

The uniqueness of JWST/NIRSpec, however, also poses new challenges for the planning and interpretation of observations of high redshift galaxies. Unlike ground-based multi-object spectrographs, which typically require the user to define a 'mask' of apertures on the sky, NIRSpec is equipped with a micro-shutter array (MSA, [START_REF] Kutyrev | Space Telescopes and Instrumentation 2008: Optical, Infrared, and Millimeter[END_REF]) providing a fixed grid of apertures on the sky. Optimising the NIRSpec MSA usage for a given scientific goal therefore requires a careful prioritization of the targets, as well as calculations to study spectral overlaps and truncations. The different filters and dispersers available on NIRSpec, covering the wavelength range λ ∼ 0.6 -5.3 µm and spectral resolutions R ∼ 100, ∼ 1000 and ∼ 2700, provide complementary information, but the optimal choice of filters, dispersers and exposure times depends on the scientific case, target properties and redshift range. Given the paucity of high quality galaxy spectra at z 4, NIRSpec data will open a largely unexplored space of observables, for both emission lines and stellar continuum studies. Maximising the information extracted from such data therefore requires models and approaches adapted to describing galaxies with a wide range of stellar populations and interstellar medium properties, likely ex-tending well beyond those measured with existing observatories. Having adequate models and simulations is therefore critical, especially for the planning of Cycle 1 and Cycle 2 observations, as these will be largely based on spectroscopic follow-up of existing, pre-JWST imaging data.

In this paper, we build a physically-motivated, semi-empirical framework to simulate integrated galaxy spectra as could be observed at low spectral resolution with JWST/NIRSpec in the Hubble Ultra Deep Field (HUDF). We then use these simulated spectra to study our ability to constrain different galaxy physical parameters, such as star formation rates, ages, mass-to-light ratios, metallicities and properties of dust attenuation and ionised gas, for objects at redshift z ∼ 4 -8. We adopt a self-consistent physical model which accounts for stellar and ionised gas emission, integrated in the stateof-the-art beagle tool (Chevallard & Charlot 2016), to generate the input mock catalog of galaxy spectra and to fit the simulated observed spectra to retrieve the galaxy properties.

In Section 2, we describe the catalogue of z 3 dropouts of Bouwens et al. (2015a) and our approach to fit the UV-to-near infrared photometry of these sources with the beagle tool. We also present our method for associating model spectra with each dropout galaxy, and demonstrate how such a method produces spectra consistent with several, independent observables. In Section 3 we present our simulations of JWST/NIRSpec observations and the full-spectrum fitting of these simulations with beagle. In Section 4 we examine the constraints on different galaxy physical parameters obtained by fitting the simulated NIRSpec data. In Section 5 we discuss our results in the light of NIRSpec observational programs, in particular of the NIRSpec Guaranteed Time Observations (GTO) program, and future extensions of this work to cover a broader range of NIRSpec observing modes. Finally, we summarize our conclusions in Section 6.

Throughout the paper, we express magnitudes in the AB system, adopt a zero-age solar metallicity Z = 0.017 (corresponding to a present-day metallicity of 0.01524, see Table 3 of Bressan et al. 2012) and the latest constraints on cosmological parameters from Planck, i.e. Ω Λ = 0.6911, Ω m = 0.3089 and H 0 = 67.74 (see last column 'TT,TE,EE+lowP+lensing+ext' of Table 4 of Planck Collaboration et al. 2015). All the emission line equivalent widths in the text refer to rest-frame values.

A SEMI-EMPIRICAL CATALOGUE OF GALAXY SPECTRA IN THE HUDF

In this section we detail our approach to create a semi-empirical catalogue of high-redshift galaxy spectra. This catalogue is constructed by matching the predictions of a spectral evolution model to a large sample of z 3 galaxies with multi-band HST photometry, and represents the first of several steps involved in the process of creating and interpreting simulated NIRSpec observations. The adoption of a semi-empirical approach is motivated by the need to minimize the model dependence of our analysis, decoupling it from a particular theoretical approach (e.g. hydro-dynamic simulation vs semi-analytic model) and its specific implementation. Also, anchoring our simulations to existing HST photometry allows us to create and study NIRSpec observations which more closely resemble those that will be obtained early on in the JWST mission, i.e. based on HST-selected targets. Adopting a purely empirical approach based on existing spectroscopic observations, on the other hand, is not viable because of the widely different wavelength coverage and sensitivity of current observatories compared

N B = 362 N V = 151 N I = 102 N Z = 66 N Y = 34
Figure 1. Distributions of the WFC3/F160W magnitude for the B (z ∼ 4, purple line), V (z ∼ 5, cyan), I (z ∼ 6, green), Z (z ∼ 7 objects, orange) and Y (z ∼ 8 objects, red) dropout galaxies selected by Bouwens et al. (2015a) in the HUDF.

to JWST/NIRSpec. Existing spectroscopic surveys targeting high redshift galaxies are often limited to observed wavelengths λ 1 µm (e.g. VVDS, Le Fèvre et al. 2013, VUDS, VANDELS), while surveys extending to λ ∼ 2.5 µm (e.g. MOSDEF) target relatively bright (m F160W 25) galaxies. We note that the semi-empirical approach adopted here, which relies on observations of UV-selected galaxies, can miss populations of galaxies with little UV emission, such as very dusty star-forming galaxies and passively evolving ones. This UV selection reflects what will actually be used to prepare early JWST/NIRSpec observations, before NIRCam observations of depth comparable to that of existing HST observations become available. In future work, we will use the recently published mock galaxy catalogue of Williams et al. (2018) to analyse NIR-Spec targets selected from deep NIRCam observations.

Multi-band HST photometry of dropout galaxies in the HUDF

In this work, we focus on galaxies at redshift z 3, for which JWST/NIRSpec will enable measurements of standard optical emission lines (e.g. Hβ, [O iii]λλ4959,5007, Hα, [N ii]λλ6548,6584, [S ii]λλ6716,6731), which are largely inaccessible with existing ground-based telescopes. As shown by the initial works of Steidel et al. (1996) and Madau et al. (1996), high redshift star-forming galaxies can be effectively selected from broad-band photometric data by using the Lyman break 'dropout' technique. This technique exploits the (nearly) complete absorption by neutral hydrogen of any light emitted by a galaxy blue-ward the Lyman limit (912 Å). As the Lyman limit is redshifted to redder wavelengths with increasing redshift, this makes an object become undetectable ('dropout') from a given optical/near-infrared band. Several groups over the years have used this technique to identify high-z star-forming galaxies, initially using optical data from the HST Advanced Camera for Surveys (ACS) to select objects out to z ∼ 6 (e.g. Steidel et al. 1999;Dickinson et al. 2004;Bunker et al. 2004), and later exploiting near-infrared observations with the Wide Field Camera 3 (WFC3) to obtain large samples of z 6 galaxies (e.g. [START_REF] Wilkins | [END_REF]Bouwens et al. 2011;McLure et al. 2013;Schenker et al. 2013).

Here we use the high-z galaxy candidates selected with the dropout technique by Bouwens et al. (2015a) in the Hubble Ultra Deep Field (HUDF). Since we are interested in galaxies out to the highest redshifts, we only consider sources selected from the 4.7 arcmin2 region in the HUDF with deep HST/WFC3 nearinfrared observations. We adopt the multi-band HST catalogue of Illingworth et al. (2013), built from a combination of all available HST/ACS and WFC3 observations in the HUDF (for a complete list see table 2 of Illingworth et al. 2013), and refer to this as the 'eXtreme Deep Field' (XDF) catalogue. The XDF catalogue includes observations in 9 bands, 5 in the optical, based on the ACS/WFC filters F435W, F606W, F775W, F814W and F850LP, and 4 in the near-infrared, taken with the WFC3 filters F105W, F125W, F140W and F160W. The 5σ depth of a point source is m F160W ∼ 29.8, computed within a circular aperture of 0.35 arcsec diameter, while the ACS + WFC3 combined image used for the source extraction reaches a 5σ depth of 31.2 within the same circular aperture (Illingworth et al. 2013). We do not use Spitzer/IRAC data, since the vast majority of our galaxies are too faint to be detected in existing Spitzer images.

We consider dropout galaxies in the filters B 435 (z ∼ 4), V 606 (z ∼ 5), I 775 (z ∼ 6), Z 850 (z ∼ 7) and Y 105 (z ∼ 8), which we will indicate as B, V, I, Z and Y dropouts in the remainder of this paper. Details on the Lyman Break selection in each band can be found in sec 3.2 of Bouwens et al. (2015a) (see also their table 2). Fig. 1 shows the F160W magnitude distribution of galaxies selected in the different dropout bands. The vast majority (∼ 80 per cent) of the galaxies at z 4.5 (V dropouts and above) have m F160W > 28, while the 50 per cent completeness magnitude varies between a minimum of ∼ 29.3 (F105W filter, I dropouts) to a maximum of ∼ 29.7 (F160W filter, Y dropouts). The adopted Lyman Break selection produces, for each dropout band, the redshift distributions shown in Fig. 2. Galaxy redshifts are centred around the average values reported above, with Full Width Half Maximum (FWHM) of the redshift distribution set by the filter widths, i.e. δz ∼ 1 for the B, V, I and Z dropouts, increasing to δz ∼ 2 for the Y dropouts. The absolute UV magnitude of the sources, computed from the median m F160W magnitude in each dropout class, assuming a flat f ν spectrum, and using the central redshift of each dropout band, varies from ∼ -19.1 (B dropouts), to ∼ -19.5 (V), ∼ -19.65 (I), ∼ -19.75 (Z), ∼ -20.3 (Y).

Broad-band SED fitting of high-redshift galaxies in the HUDF

Similarly to Chevallard & Charlot (2016, see their sec 4), we use the beagle tool to fit the XDF photometry of 715 dropout galaxies with a self-consistent physical model that includes stellar emission, continuum+line emission from H ii regions and diffuse ionized gas and dust attenuation. We do not model the emission from an AGN potentially contaminating the HST photometry, as the expected number of z > 3 type-1 AGNs in the 4.7 arcmin 2 field here considered is consistent with zero (e.g. see section 4.2.5 of Grazian et al. 2015). We let the redshift free to vary, and adopt a two-component star formation history constituted by a 'smooth' function and a burst. The smooth component is described by a delayed exponential function ψ(t ) ∝ t exp (-t /τ sfr ), where τ sfr is the star formation timescale and t the age of the galaxy, taken to lie between 10 7 yr and the maximum time allowed since the onset of star formation at the galaxy redshift (see below). The burst (with constant star formation rate) covers the last 10 7 yr of star formation, the timescale over which ∼ 99.9 per cent of the H-ionizing Fall 1993;Binette et al. 1994), and it is parametrised in terms of the 'current star formation rate' ψ, i.e. the SFR averaged over the past 10 Myr. Decoupling the 'past' star formation history and the 'current' star formation rate allows us to obtain galaxy spectra with any contribution of emission lines relative to stellar continuum. We fix the maximum redshift for the formation of the first stars in a galaxy at z max form = 15, so that at any redshift z the maximum allowed time since the onset of star formation is t max form (z) = t u,zt(z max form ), where t u,z and t(z max form ) refer to the age of the Universe at redshift z and z max form , respectively. We approximate the distribution of stellar metallicities in a galaxy, including the metallicity of young stars (with ages t 10 7 yr, Z young ) and the interstellar metallicity (Z ism ), with a single metallicity Z, i.e. Z young = Z ism = Z. 2 Following Gutkin et al. (2016, see also Charlot & Longhetti 2001) we describe the properties of gas ionized by young stars by means of galaxy-wide ('effective') parameters. The ionization parameter log U S determines the ratio of H-ionizing photons to gas density at the Strömgren radius of an effective star cluster, while the dust-to-metal mass ratio ξ d ('depletion factor') sets the amount of depletion of heavy elements onto dust grains (see sec 2.3 of Gutkin et al. 2016 for a discussion of depletion factors). 3 We model the effect of dust attenuation on stellar and gas emission by appealing to the two-component model (diffuse ISM + birth clouds) of Charlot & Fall (2000), parametrized in terms of the total attenuation optical depth τV and the fraction of attenuation arising in the diffuse ISM µ. We account for the effect of absorption from the intergalactic medium (IGM) by means of the average prescription of Inoue et al. (2014).

Following the Bayesian approach adopted in beagle, we define the posterior probability distribution of the model free param-

P(Θ | D, H) ∝ π(Θ) L(Θ) , (2.1)
where D indicates the data, H the adopted model, π(Θ) the prior distribution and L(Θ) the likelihood function. We adopt independent priors for all parameters, uniform for the parameters z, log(M/M ), log(τ sfr /yr), log(Z/Z ), log U S and ξ d , Gaussian for log(t/yr) and log(ψ/M yr -1 ) and exponential for τV (see Table 1). 4 We consider a multi-variate Gaussian likelihood function with independent errors σ i on each measurement y i

-2 ln L(Θ k ) = i y i -ŷi (Θ k ) σ i 2 , (2.2)
where the summation index i runs over all observed bands (i.e., bands with positive errors, and negative or positive fluxes), y i is the observed flux, σ i = σ obs, i 2 + (σ 0 y i ) 2 , where σ obs, i is the observational error and σ 0 = 0.02 is an additional error term that we add to avoid obtaining results dominated by systematic uncertainties (see sec. 4.2 of Chevallard & Charlot 2016) and ŷi ( Θ k ) indicates the fluxes predicted by our model for a set of parameters Θ k . We adopt the Nested Sampling algorithm (Skilling et al. 2006) as implemented in multinest (Feroz et al. 2009) to sample the posterior probability distribution of the 9 free parameters

Θ xdf = [z, log(M/M ), τ sfr , t, Z, ψ, log U S , ξ d , τV ].
It is worth briefly pausing to discuss the potential risk of over-fitting our data by using such a flexible physical model. The 9 photometric bands used in the fitting cover the wavelength range 0.4 λ/µm 1.6, hence they mainly probe the rest-frame UV emission of z 3 galaxies, and only bands redward the Lyman break (4 bands for the Y dropouts) provide constraints on the galaxy physical properties. In star-forming galaxies, the UV continuum emission is mostly sensitive to the recent ( 10 8 yr) star formation history and to dust attenuation, while other parameters, such as the mass of older stars and the physical conditions of gas, have little influence on the emission at these wavelengths. Since in this work we aim at simulating NIR-Spec observations covering a large variety of intrinsic galaxy spectra, extending beyond those observed in relatively bright galaxies at z 4, we must also account for the variation of model parameters largely unconstrained by existing photometric observations. For this reason we adopt the flexible, 9-parameters model described above, and combine the weak constraints provided by HST photometry on some model parameters with well established relations among galaxy physical quantities to obtain physically-motivated combinations of parameters (see Section 2.3). We then validate this approach by comparing our model predictions with external (photometric and spectroscopic) data-sets at redshift z ∼ 2 -8 (Section 2.4).

The results of the beagle fitting of the XDF data are summarised in Fig. 2 in terms of the stacked posterior probability distribution of redshift. This is computed by combining the redshift probability distribution of each galaxy using a kernel-density estimator. Similarly to an histogram, Fig. 2 represents a density plot, where the solid curves depend on both the density of galaxies at each redshift and on the redshift probability distribution of each individual source. Fig. 2 provides an interesting visual comparison among the redshift distributions obtained by Bouwens et al. (2015a) from the analysis of Monte Carlo simulations of artificial sources (see their section 4.1) and those computed in this work. Given the very different approaches adopted to compute the two sets of lines of Fig. 2, the figure highlights a good agreement between the two methods. The differences among the two approaches are larger for the Y dropouts, for which the peaks of the two distributions are separated by (δz ∼ 0.5), possibly because of the low constraining power of data for the Y dropouts (low S/N, few bands red-ward the continuum break) combined with a low number of sources (N Y = 34).

Linking model spectra and observed galaxies in the HUDF

The beagle tool provides us with the posterior probability distribution of the 9 free parameters of the model (reported in column 5 of Table 1), along with sets of selected observables (e.g., full spectrum, absolute and apparent magnitudes) and derived quantities (e.g. rates of H-and He-ionizing photons). These are obtained by sampling the posterior probability distribution of model parameters using the multinest algorithm (see section 3.3 of Chevallard & Charlot 2016 for more information on the beagle output). From a purely statistical perspective, we could associate with a given XDF source any spectrum from the set drawn using multinest in this way. However, since the data are consistent with a broad range of model spectra, and a correspondingly broad range of model parameters, this approach would not ensure that the model spectra match other, independent observables (e.g. galaxy colours at longer wavelengths, emission line measurements), nor that the selected combinations of model parameters are consistent with measured relations among galaxy physical quantities (e.g. the massmetallicity relation). Obtaining a physically-motivated distribution of emission line strengths is particularly relevant for XDF galaxies, which exhibit faint continuum emission and will likely be detected only via emission lines with JWST/NIRSpec. By obtaining such a physically-motivated distribution of emission lines, we can produce simulated NIRSpec observations with realistic S/N ratios. This, in turn, should ensure that the constraints on the galaxy physical parameters derived in Section 3 can be applied to future 'deep' NIRSpec observations of HST-selected high-redshift galaxies.

After some experimentation, we find that requiring for the fitted galaxies to follow three well-established relations between galaxy physical parameters produces observables in agreement with different types of observations. These relations include a relation between (i) stellar mass and star formation rate ('main sequence' of star-forming galaxies); (ii) stellar mass, star formation rate and gas-phase metallicity ('fundamental' metallicity relation); (iii) interstellar metallicity and ionization parameter. We adopt the (redshift-dependent) relation between stellar mass and star formation rate of Speagle et al. (2014), which they derive from a broad compilation of stellar mass and star formation rate measurements at redshift z ∼ 0 -6. Their relation (see their equation 28) can be expressed as

log(ψ/M yr -1 ) =[0.84 -0.026 (t /Gyr)] log(M * /M )- [6.51 -0.11 (t /Gyr)] ,
(2.3)

where t indicates the age of the Universe at redshift z.

Following Williams et al. (2018), we consider the relation between stellar mass, star formation rate and gas-phase metallicity of [START_REF] Hunt | [END_REF], which they obtain from the analysis of the UVto-infrared photometry and optical emission lines of a sample of ∼ 1000 star-forming galaxies at redshift z ∼ 0 -3.7. The relation of [START_REF] Hunt | [END_REF] can be expressed as 12 + log(O/H) = -0.14 log (ψ/M yr -1 ) + 0.37 log (M * /M ) + 4.82 .

(2.4)

Finally we adopt the relation between interstellar metallicity and ionization parameter of Carton et al. (2017), which they derive by adopting the 'theoretical' dependence of ionization parameter on interstellar metallicity of Dopita et al. (2006) (log U S ∝ 0.8 log Z ism ), and computing the zero-point of this relation by using metallicities and ionization parameters of SDSS DR7 galaxies estimated by Brinchmann et al. (2004). The relation between log U S 2017) at z ∼ 5 -6 from ∼ 50 galaxies selected from 4 gravitationally-lensed HST Frontier Fields, while the dashed line is an extrapolation of this relation at lower stellar masses. and Z ism can then be expressed as log U S = -0.8 log (Z ism /Z ) -3.58 .

(2.5)

In practice, we enforce these relations by multiplying the posterior probability distribution of equation (2.1) by three 'weight' functions which depends on the above parameters. In order to more densely populate the tails of these relations, we adopt a Student's-t weight function with 3 degrees of freedom, which is significantly broader than a Gaussian function. The Student's-t function is expressed as

f (x) = 6 √ 3 π 3 + x 2 2 1 σ x , (2.6)
where x is the standardised variable

x = x -x σ x , (2.7)
and x indicates the value of the model parameter in the posterior distribution computed with beagle and x the value of the parame- ter predicted by equations (2.3)-(2.5). We fix the scatter to σ x = 0.3 for the mass-star formation rate relation (see, e.g., Speagle et al. 2014;Shivaei et al. 2015), and to σ x = 0.2 for the 'fundamental' metallicity relation and for the metallicity-ionization parameter relation. Equations (2.3)-(2.6) provide the conditional distribution of a model parameter (star formation rate, gas-phase metallicity, ionization parameter) given other model parameters, for example equations (2.3) and (2.6) define the conditional probability P(ψ | M * , z). This conditional distribution does not depend on the data D, hence it can be readily interpreted as a conditional prior distribution linking different model parameters.

It is well known that the photometric degeneracy between the Lyman and Balmer breaks can make a galaxy photometric SED be equally well described by two different models with widely different redshifts (e.g. Ilbert et al. 2006). As the beagle tool can identify these multiple redshift solutions, in this work we restrict the analysis to redshifts consistent with the expected dropout redshift of Bouwens et al. (2015a). To achieve this, among the combinations of model parameters Θ obtained with beagle, we only consider those satisfying the condition z > z min , where z min = 2 for the B dropouts, z min = 3 for the V, z min = 4 for the I and z min = 5 for the Z and Y ones (Fig. 2). We then randomly draw, for each galaxy, a set of parameters among those sampled by multinest, where the probability of drawing any set is proportional to the re-weighted posterior probability distribution. We can repeat the random draw (with replacement) N times to obtain N model spectra consistent with the observed HST photometry, and corresponding to combinations of [M * , ψ] satisfying the adopted mass-star formation rate relation. This is particularly relevant to evaluate selection effects and to study the consequences of deriving population-wide relations among physical parameters (e.g. mass-metallicity, mass-star formation rate) from a finite number of observations. Both applications will be part of successive works, while in the remainder of this paper we consider a single random draw to associate a model spectra to each XDF source, and refer to the ensemble of 715 model spectra as a 'mock catalogue'.

We show in Fig. 3 by means of (grey) density contours on a logarithmic scale the two-dimensional stacked posterior probability distribution of stellar mass M * and star formation rate ψ, along with the pairs [M * , ψ] drawn from the re-weighted posterior probability distribution (circles of different colours). By analogy with the onedimensional stacked posterior probability distribution of redshift (see Section 2.2), we compute the two-dimensional stacked posterior probability distribution of M * and ψ summing the individual probability distributions computed with a kernel-density estimator. Fig. 3 shows the effect of varying in our model both the current star formation rate and total stellar mass when only rest-frame UV observations are available. The observed HST photometry of a large fraction of XDF galaxies can be reproduced equally well by our model for widely different combinations of M * and ψ, which makes the stacked distribution occupy a broad region of Fig. 3. We note that the sharp upper envelope created by the coloured circles in Fig. 3 is caused by the fixed duration (10 Myr) of the current burst of constant SF, which limits the maximum ψ attainable at fixed M * .

Fig. 4 shows the distribution of the main model free parameters of the mock catalogue. The low constraining power of HST photometry with respect to most model parameters makes these distributions reflect the adopted priors, i.e. Gaussian for ψ, exponential for τV and uniform for the other parameters (see Table 1). Also, no strong correlations are visible except for the relation between log M * and log ψ that we enforced. It is worth briefly discussing the absence in Fig. 4 2012) stacks and our mock spectra. For each band, the shaded regions of different colours are slightly shifted horizontally for clarity. (b) Same as (a), but for V dropout galaxies (redshift ∼ 5). We did not plot the median SED for the F850W ∼ 26 bin because of the low (4) number of galaxies in our mock catalogue falling in the bin.

et al. 2004), while observations targeting galaxy populations similar to those used in our analysis, i.e. faint high-redshift galaxies, depict a less clear picture. Recent results from the VUDS survey suggest that highly star-forming galaxies with stellar masses comparable to those of our mock galaxies (log M/M * ∼ 6.5 -9.5) span a broad range of gas-phase metallicities 12 + log(O/H) ∼ 7.5 -8.5 at fixed stellar mass (Calabrò et al. 2017). Similar results have been obtained by Izotov et al. (2015) and Jimmy et al. (2015) from the analysis of low-mass local 'analogues' of high redshift galaxies. These findings can have a physical origin, since the metal content of star-forming gas results from the complex interaction between enrichment from previous stellar generations, AGN and stellar feedback and gas accretion, or can be driven by systematic uncertainties in the metallicity estimators (e.g. Kewley & Ellison 2008).

Similarly, an anti-correlation between gas-phase metallicity and ionization parameter has been observed at low redshift (Carton et al. 2017), and can be expected on theoretical grounds (e.g. Dopita et al. 2006), but the redshift evolution and dispersion of the relation for galaxy populations as those considered in this work are largely unknown. While future observations may show that certain combinations of model parameters (e.g. high metallicity and large ionization parameter) are not well suited to describe the conditions of photoionized gas in high-redshift star-forming galaxies, we have decided not to impose any relation between these parameters in our mock galaxy catalogue. This will enable users to adopt any relation between mass-metallicity and metallicity-ionization parameter to select sub-samples of sources from the different realisations of the mock catalogue.

Comparison of model spectra with independent data

Our approach to build a catalogue of mock galaxy spectra guarantees, by construction, that these match the HST photometry of XDF dropouts. In this section, we compare the mock spectra with other observables, namely near infrared photometry from the Infrared Ar-ray Camera (IRAC) onboard the Spitzer Space Telescope, and measurements of the [C iii]λ1907+C iii]λ1909 emission line performed with different instruments. The lower sensitivity of Spitzer with respect to HST only enables the detection in the IRAC band 1 (3.6 µm filter, IR 1 ) and band 2 (4.5 µm filter, IR 2 ) of a minority of XDF dropouts with m F160W 27 (see Fig. 1). For this reason, we have to rely on Spitzer observations of stacked galaxies and of individual (brighter) galaxies over a much wider area than the XDF region.

We compare in Fig. 5 the predicted optical-to-near infrared SEDs (HST +Spitzer bands) of our mock catalogue with the stacked SEDs computed by González et al. (2012) at z ∼ 4 and z ∼ 5 galaxies. González et al. (2012) consider HST/ACS, HST/WFC3 and Spitzer/IRAC observations in the GOODS-South region (Giavalisco et al. 2004;Bouwens et al. 2012), and compute the median SEDs of B, V, I and Z dropouts in bins of observed UV luminosity. In practice, they compute the median HST fluxes of the sources in each UV luminosity bin, and perform aperture photometry on the median-combined image in each IRAC band. They evaluate the uncertainties on the median fluxes through bootstrap resampling. We note that when comparing model predictions with Spitzer/IRAC fluxes of z 4 galaxies, one must consider the effect of optical emission lines contaminating the IR 1 and IR 2 bands. At z ∼ 3 -5 (B dropouts), Hα + [N ii] + [S ii] contaminates the IR 1 band, while at z 5 both IRAC bands can be contaminated by either Hβ

+ [O iii] or Hα + [N ii] + [S ii].
This makes the IRAC colours computed from the stacked SEDs of z 5 galaxies highly dependent on the emission line properties and redshift distribution of the sources entering the stacks, and this dependence is exacerbated in stacks computed from a low number of sources (e.g. see fig. 7 of González et al. 2012). For this reason, we only consider the z ∼ 4 and z ∼ 5 stacked SEDs of González et al. (2012) IRAC fluxes, unlike the HST ones, are not matched to any observation during the mock catalogue creation. Fig. 5 hence demonstrates that, on average, the mass and ages of evolved stars in our mock galaxies, traced by the IRAC fluxes probing the rest-frame opticalto-near infrared emission of galaxies, are consistent with the observed values.

The Spitzer/IRAC colour IR 1 -IR 2 can be used to further test the agreement between observations and mock spectra. As we noted above, the redshift evolution of this colour is sensitive to the presence of optical emission lines, especially the two groups Hβ + [O iii] and Hα + [N ii] + [S ii], contaminating the IR 1 and IR 2 bands (e.g. observationally, Stark et al. 2013;de Barros et al. 2014;theoretically, Wilkins et al. 2013). We therefore compare in Fig. 6 the evolution of the IR 1 -IR 2 colour in the redshift range 2 z 8 predicted by our mock catalogue with data from different sources. These include, from low to high redshift, observations from 3D-HST (Skelton et al. 2014), from the spectroscopic sample of Smit et al. (2016) and from the photometric samples of [START_REF] Rasappu | [END_REF], Smit et al. (2015) and Smit et al. (2014). Since our mock catalogue is based on observations from the small, 4.7 arcmin 2 XDF area, while data are drawn from much larger areas ( 100 arcmin 2 ), we encode in the circle sizes the rest-frame UV luminosity M uv of each object: the larger the circle, the more luminous the galaxy. Fig. 6 highlights the ability of the mock spectra to reproduce the sharp IRAC colour change at z ∼ 3.8 (z ∼ 5) caused by the entry of Hα in the IR 1 (IR 2 ) band, including the most extreme blue and red colours, which are caused by EW(Hα + [N ii] + [S ii]) ∼ 1000. We note that the different extremes reached by the IRAC colours at z ∼ 3.8 -5 (IR 1 -IR 2 ∼ -1) and at z ∼ 5 -6 (IR 1 -IR 2 ∼ 0.8) are likely caused by the different widths of the IRAC filters, ∼ 0.75 µm for IR 1 filter and ∼ 1 µm for IR 2 one. This difference translates into a different contamination of the integrated broad-band flux of each band, at fixed emission line equivalent width (EW). At redshifts z ∼ 5.5 -6.6 the IR 1 band is contaminated by the group of lines Hβ + [O iii] and the IR 2 band by Hα + [N ii] + [S ii], therefore the IR 1 -IR 2 colour depends on the relative intensities of H-Balmer lines vs [O iii]λλ4959,5007, i.e. on the physical conditions of ionized gas, especially metallicity and ionization parameter. In the narrow redshift window z ∼ 6.6 -6.9, the IR 2 band is free of strong emission lines, while IR 1 is contaminated by Hβ + [O iii]. Smit et al. (2015) exploit this property to select extreme emission lines galaxies with accurate photometric redshifts. They search for such extreme objects in the 5 CANDELS fields (∼ 900 arcmin 2 ), and find ∼ 20 sources with IR 1 -IR 2 -1. These are relatively bright galaxies falling in a narrow redshift range, hence, not surprisingly, they do not appear in our mock catalogue, which reaches at most IR 1 -IR 2 ∼ -0.6 at z ∼ 6.6 -6.9. Overall, Fig. 6 indicates that the strengths of the Hα + [N ii] + [S ii] and Hβ + [O iii] lines in our mock catalogue are consistent with those inferred from Spitzer/IRAC colours of z ∼ 2 -8 galaxies. This validation is particularly important since the strength of the optical emission lines in our mock catalogue will translate into a distribution of S/N in the NIRSpec simulations (see Section 3.2), therefore directly affecting the constraints on galaxy physical parameters form NIRSpec spectra discussed in Section 4.

The above tests only provide indirect constraints on spectroscopic features through the contamination of broad band filters by emission lines. As we have already noted, existing observatories do not allow the detection of optical emission lines at z 4, hence to study the evolution of galaxy spectral features across the widest redshift range one must rely on UV emission lines. The Lyα line is the most luminous emission line at UV wavelengths, but radiative transfer effects caused by its resonant nature make the comparison of model predictions and observations non trivial (e.g., Verhamme et al. 2006). We therefore consider another bright UV line, the doubly-ionized carbon line [C iii]λ1907+C iii]λ1909, which has been observed both at low (e.g. Berg et al. 2016;Senchyna et al. 2017) and high redshift (e.g. Erb et al. 2010;Stark et al. 2015). We compare in Fig. 7 the redshift evolution of the [C iii]λ1907+C iii]λ1909 equivalent width from literature with that predicted by our mock catalogue. In the same figure, we also show with a histogram the distribution of [C iii]λ1907+C iii]λ1909 equivalent widths in our mock catalogue, and the distribution derived by Le Fèvre et al. (2017) from observations at redshifts z ∼ 2 -3.8. Fig. 7 indicates that while most galaxies in our catalogue exhibit W C iii] 5 Å, a significant number of mock spectra attain substantially larger W C iii] values, reaching the most extreme W C iii] 20 Å found by Stark et al. (2015Stark et al. ( , 2017) ) and Amorín et al. (2017). Interestingly, Fig. 7 also reveals a larger fraction of galaxies with [C iii]λ1907+C iii]λ1909 equivalent widths 10 Å in our catalogue relative to that in the sample of Le Fèvre et al. (2017). This is likely a consequence of the higher redshifts and lower luminosities of our mock galaxies with respect to the sample of Le Fèvre et al. (2017) (z 4, i AB 25), which imply, on average, younger and more metal poor stellar populations, able to power stronger [C iii]λ1907+C iii]λ1909 emission.

The comparisons of our mock spectra with external observables presented in this section validate our semi-empirical approach to build a mock galaxy catalogue based on HST photometry in the HUDF. In particular, the above tests demonstrate the ability of our catalogue to match the rest-frame optical continuum emission of z ∼ 4 and z ∼ 5 galaxies (Fig. 5), the contamination of the strongest optical emission lines to Spitzer/IRAC bands (Fig. 6), and the observed range of [C iii]λ1907+C iii]λ1909 equivalent widths at z ∼ 1 -8 (Fig. 7).

SIMULATING AND FITTING JWST/NIRSPEC OBSERVATIONS

In this study, we focus on simulations of deep (100 ks) observations performed with JWST/NIRSpec. We consider the multi-object spectroscopy (MOS) mode and adopt the low spectral resolution configuration NIRSpec/prism (PRISM/CLEAR spectral configuration), which enables the simultaneous observation of up to ∼ 200 objects over the full spectral range λ ∼ 0.6 -5.3 µm.

Simulating JWST/NIRSpec spectra

A detailed description of the approach adopted to simulate NIR-Spec spectra can be found in Appendix A. Here, we only provide an overview of our approach and highlight a few key features of the simulations. For each galaxy, we start from the (noiseless) mock spectrum obtained using the procedure described in Section 2. This high-resolution spectrum is redshifted to the photometric redshift obtained with the beagle analysis and then rebinned to the spectral pixel size of the NIRSpec/prism configuration. We then use an idealized model accounting for the response of the telescope and of the instrument to derive the number of electrons per-second generated at detector level in each pixel. Combined with a noise model, this allows us to generate a simulated (noisy) spectrum. This approach is very similar to the one used by many exposure-time calculators (ETCs), but simpler than the elaborate two-dimensional scheme used by JWST's official ETC 'Pandeia' [START_REF] Pontoppidan | Observatory Operations: Strategies, Processes, and Systems VI[END_REF].

Hβ + [O iii] λ4959 + [O iii] λ5007 [O iii] λ4959 + [O iii] λ5007 [O iii] λ4959 [O iii] λ5007 Hα + [N ii] λλ6548, 6584 Hα + [N ii] λλ6548, 6584 + [S ii] λλ6716, 6731 [S ii] λλ6716, 6731
In our simulations, we also account for the fact that a significant fraction of the object light may fall outside of the standard 3shutter slitlet, especially for extended objects (the projected size of the aperture of an individual micro-shutter is 0.20 × 0.42 arcsec). In practice, we compute the (wavelength-dependent) throughput of an extended source with effective radius r e , described as an exponential profile with Sérsic index 1.5 (Shibuya et al. 2015), and average over all source positions within the aperture of the central shutter of the slitlet. Averaging over random ellipticities and position angles produces small corrections ( few per cent), which we therefore ignore. Following the approach of Williams et al. ( 2018 where r 0 e = 6.9 × (1 + z) -1.2 kpc and L 0 uv is the UV luminosity corresponding to a UV magnitude M 0 uv = -21. This leads to median sizes of 0.9 (0.13), 0.75 (0.12), 0.63 (0.11), 0.54 (0.10) and 0.48 (0.10) kpc (arcsec) for the B, V, I, Z, Y dropouts, respectively. The total throughput at λ = 2.5 µm corresponding to the median sizes above is between 34 per cent (B dropouts) and 38 per cent (Y dropouts) for the aperture and light profile adopted.

The resolution R = λ/∆λ of the NIRSpec prism has a strong The red line and red shaded region indicate the posterior median and 95 per cent credible interval, respectively, from the beagle fit of the simulated spectrum.

1 2 3 4 5 λ/µm (observed-frame) 0.0 0.8 1.6 2.4 F λ /(erg s -1 cm -2 Å -1 ) ×10 -20 [ O i i ] H δ H γ H β [ O i i i ] H α (
F λ /(erg s -1 cm -2 Å -1 ) ×10 -20 [ O i i ] H δ H γ H β [ O i i i ] (b) 0.4 0.8 1.2 1.6 λeff/µm
(b) same as (a), but for a Z dropout (XDFZ-3312565447, m F160W = 28.74), with z = 7.176, log (M * /M ) ∼ 8.9, ψ/M yr -1 ∼ 5.4 and log (ψ s /yr -1 ) ∼ -8.3. The small inset in each panel shows the observed XDF photometry (blue diamonds, from Bouwens et al. 2015a) along with the model photometry predicted by beagle for the set of physical parameters associated to these sources with the procedure outlined in Sec. 2.3 (orange stars).

wavelength dependence, from a minimum of R ∼ 30 at λ = 1.2 µm to a maximum of R ∼ 300 at λ = 5 µm. This makes groups of neighbouring emission lines appear blended or unblended depending on a galaxy redshift. Since UV and optical emission lines will provide the firmest constraints on the physical parameters of faint galaxies at high redshift, we have summarised in Fig. 8 the visibility and blending/unblending of different lines as a function of redshift. Emission lines in Fig. 8 are considered unblended when their redshifted wavelengths are separated by a ∆λ corresponding to 2.2 detector pixels, i.e. the typical instrumental spectral resolution element size of an extended source which uniformly fills the micro-shutter aperture of width ∼ 0.2 arcsec (along the spectral direction). We note that the ETC simulations performed here do not account for the impact of the object size on the spectral response, implying that emission lines always fall into a single detector pixel (see Appendix A3). This approximation has a minor impact on the results of our analysis, as the only strong lines potentially affected are the two components of [O iii]λλ4959,5007, whose ratio is fixed by atomic physics, and the Hα and [N ii]λλ6548,6584 lines. The low metallicity of our galaxies and absence of an AGN component in our models imply that [N ii]λλ6548,6584 is in any way too weak to provide any constraint on the physical properties of mock galaxies. 10 provides an overview of NIRSpec/prism observations of 10 galaxies at z ∼ 4 -8 covering a broad range of continuum luminosities m F160W ∼ 27.9 -29.7. The figure highlights how the large star formation rates per unit stellar mass and young ages of high redshift galaxies should enable statistically significant detections of emission lines even in m F160W 29 galaxies out to z ∼ 9.

Equivalent width and S/N distribution of emission lines

We show in Fig. 11 the distribution of equivalent widths of Hβ and of the two groups of lines Hβ + [O iii] and Hα + [N ii] + [S ii] computed from the mock spectra. The mean (median) redshift of the galaxies in the mock catalogue is z ∼ 4.8 (∼ 4.5), and the mean (median) equivalent width of Hα + [N ii] + [S ii] ∼ 450 Å (∼ 600 Å), which compares favourably with the value ∼ 400 Å found by Smit et al. (2016) at redshift z ∼ 3.8 -5.0, and of ∼ 550 Å and ∼ 600 Å found by [START_REF] Rasappu | [END_REF] at redshift 5.1 z 5.4 in their photometric and spectroscopic samples, respectively. Fig. 12 shows the relation between the signal-to-noise ratio of the Hβ line S/N Hβ and the observed F160W magnitude for all galaxies in the mock catalogue, while Table 2 indicates the fraction of galaxies in bins of F160W magnitude with S/N Hβ > 3, 5 and 10. The S/N reported in Fig. 12 and Table 2 accounts for the effect of the instrumental spectral response, while our simulations do not account for this effect. The consequence is, on average, a ∼ 35 per cent larger emission line S/N in our simulated spectra with respect to a computation accounting for the spectral response (see Appendix A).

Fig. 12 and Table 2 indicate that NIRSpec/prism 'deep' observations will detect most z 4 galaxies in the HUDF through their emission lines, with typical S/N Hβ ∼ 3 -30. We expect more than 40 per cent of the XDF dropouts with m F160W < 28.5 to exhibit S/N Hβ > 5, and to detect with S/N Hβ > 3 about one third of the faintest (m F160W > 30) galaxies. The reported S/N Hβ values imply an even larger S/N of Hα (visible out to z ∼ 7), and likely of [O iii]λ5007 (out to z ∼ 9.5), given the large [O iii]/Hβ will discuss in Section 4, the simultaneous detection at high S/N of several hydrogen and metal emission lines enables tight constraints on several key galaxy properties, such as star formation rates and gas-phase metallicities.

Spectroscopic redshifts determination

In the fitting of the simulated NIRSpec spectra, which we describe in the next section, we assume a perfect knowledge of the spectroscopic redshift of the galaxies. Here we therefore discuss how well spectroscopic redshifts can be measured from deep NIRSpec/prism observations. To achieve this, we compute for each simulated spectrum the signal-to-noise ratio of the main UV and optical emission lines that NIRSpec can observe (see Fig. 8 for a list of lines). We then split the galaxies in the redshift and magnitude bins reported in Table 3, and report in the same Table the fraction of objects in each bin for which we detect at least 1 or 2 emission lines with S/N 3 and 5. We consider that objects with at least 2 emission line detections with S/N 3 or with at least 1 emission line detection with S/N 5 will have secure spectroscopic redshifts, since in most cases HST photometry will enable us to distinguish single line detections as being either [O iii]λλ4959,5007 or Hα. Table 3 hence shows that out to redshift z ∼ 7 the vast majority ( 90 per cent) of the galaxies with m F160W 30 will have secure spectroscopic redshifts, while this fraction drops to ∼ 70 per cent for fainter objects with m F160W > 30. At z 7 the Hα line is shifted outside the NIRSpec range, hence the fraction of galaxies with secure spectroscopic redshifts drops to ∼ 70 per cent for objects with m F160W 30, while the number of galaxies in our catalogue with m F160W > 30 and z > 7 is too low to provide a meaningful estimate of this fraction. Overall, Table 3 demonstrates that only a few tens of objects in our catalogue may not have secure spectroscopic redshifts, i.e. ∼ 30 per cent of m F160W > 30 galaxies at z 7 and ∼ 30 per cent of m F160W > 28 galaxies at z > 7.

Fitting simulated JWST/NIRSpec spectra

We use the beagle tool to perform a full, pixel-by-pixel fitting of each simulated NIRSpec spectrum. We adopt a model similar to the one used to fit the HST/XDF photometry (see Section 2.2), but fixing the redshift to the estimated value to mimic the way data will be analysed, i.e. with redshift being determined directly from the position of the spectral lines, or by combining single line detections with photometric redshifts (see Section 3.3 above). Unlike in the HST broad-band fitting, we let the interstellar metallicity Z ism and fraction of V-band attenuation arising in the diffuse ISM, µ, free to vary. We adopt the same independent priors as in Section 2.2, i.e. uniform for log(M/M ), log(τ sfr /yr), log(Z/Z ), log U S , ξ d , log(Z ism /Z ) and µ, exponential for τV and Gaussian for log(t/yr) and log(ψ/M yr -1 ) (see Table 1). As in Section 2.2, we consider a multi-variate Gaussian likelihood function, assuming independent errors in each pixel (see however Appendix A6). The likelihood function can therefore be described by equation (2.2), considering the summation over all spectral pixels, and assuming σ i = σ obs, i , i.e. removing any additional error term since, by construction, we have a perfect knowledge of the noise properties of the simulated spectra. We defer to a future work, based on the NIRSpec Instrument Performance Simulator, a more detailed analysis of other sources of random and systematic uncertainties in JWST/NIRSpec data, related, for instance, to detector defects, flatfielding and spectro-photometric calibration (see Appendix A6).

We show in Fig. 9 Table 4. Error thresholds on the retrieved model parameters corresponding to, 0.2, 0.5 and 0.9 of the inverse cumulative distribution function (quantile function) of the errors σ θ , for the different dropout classes. The error σ θ is defined as half of the 68 per cent central credible interval (1σ for a Gaussian distribution).

(b) (b) (b) (b) (b) -1 0 1 log(ψ c /M yr -1 ) (input) -1 0 1 log(ψ c /M yr -1 ) (c) (c) (c) (c) (c) -4 -3 -2 -1 log U S (input) -4 -3 -2 -1 log U S (d) (d) (d) (d) (d) 7 8 9 log[O/H] (input) 7 8 9 log[O/H] (e) (e) (e) (e) (e) 0 1 2 3 τv (input) 0 1 2 3 τv (f) (f) (f) (f) (f)
Each cell indicate the error threshold s θ on the parameter θ such that the cumulative distribution C(σ θ s θ ) = 0.2, 0.5 and 0.9. cent central credible interval computed from the same distribution of fluxes. Fig. 9(a) and (b) show how the high-S/N detections of several emission lines constrain the model predictions to high precision. As we will discuss in the next section, this in turn induces tight constrains on several physical quantities describing the conditions of photoionized gas in these galaxies.

Parameter s 0.2 θ s 0.5 θ s 0.9 θ B V I Z Y B V I Z Y B V I Z Y log(M/M ) 0.

JWST/NIRSPEC CONSTRAINTS ON GALAXY PROPERTIES

We summarise in Fig. 13 and Table 4 the main results of our analysis of simulated 'deep' NIRSpec observations of high redshift galaxies. We report the constraints on the 6 quantities stellar mass log(M/M ), galaxy age log(t/yr), current star formation rate log(ψ/M yr -1 ), ionization parameter log U S , gas-phase metallic-ity 12 + log(O/H) and V-band attenuation optical depth τV . Fig. 13 shows in different panels the relation between retrieved and input value for the 6 quantities above. The large errorbars visible in Fig. 13(a) in galaxies with log(M/M ) 8 , along with the 'flattening' of the relation between input and retrieved value, highlight the large uncertainty affecting the determination of the stellar masses, for which the posterior distribution is dominated by the prior. The stellar continuum emission in these galaxies, related to the stellar mass, is in fact too faint to be detected in our simulated spectra. Similarly, the large errorbars and spread of the points in Fig. 13(b) show the difficulty of determining the age of the oldest stars from our simulated NIRSpec spectra. The reason is that the UV-to-optical spectra of galaxies with large specific star formation rates are dominated by the emission from young stars, which outshine the older ones. Fig. 13(c), (d), (e) and (f) show that the star formation rate, ionization parameter, gas-phase metallicity and at-tenuation optical depth can be constrained across the entire range spanned by the input parameters. The constraints on these quantities come from emission lines ratios, and are therefore largely unaffected by our ability to detect the stellar continuum emission in faint, high-redshift galaxies.

We quantify in Table 4 the precision of the constraints on the 6 physical parameters shown in Fig. 13. For this, we compute for each parameter the cumulative distribution of the errors σ θ , and report in Table 4 the errors s θ corresponding to 0.2, 0.5 and 0.9 of the inverse cumulative distribution C -1 , i.e. s 0.2 θ , s 0.5 θ and s 0.9 θ . The values s θ therefore indicate the error on the parameter such that 20, 50, or 90 per cent of the objects exhibit errors σ θ s θ . The central column of Table 4 indicates that for ∼ 50 per cent of the simulated galaxies the stellar masses (or mass-to-light ratios) can be constrained within a factor of ∼ 2.5 (∼ 0.4 dex), the age of onset of star formation within a factor of ∼ 3.5 (∼ 0.55 dex), the star formation rate, ionization parameter and gas-phase metallicity within a factor of ∼ 1.5 (∼ 0.2 dex), albeit the precision worsens for galaxies at z 7 (Z and Y dropouts). The V-band dust attenuation optical depth can be constrained with a precision σ τV 0.4.

The model parameters not reported in Fig. 13 andTable 4, namely the timescale of star formation log(τ sfr /yr), dust-to-metal mass ratio ξ d and fraction of dust attenuation arising in the diffuse ISM µ, are largely unconstrained by this analysis. As for the age of the oldest stars, constraining τ sfr and µ would require detecting at high S/N the rest-frame optical-to-near infrared (stellar) continuum emission of the galaxies, which traces older stellar populations. Since the NIRSpec observations here considered mainly provide high-S/N detections of optical emission lines, a way to better constrain these parameters would be to combine NIRSpec spectra with JWST/NIRCam and MIRI photometry. This type of analysis will be included in a future study. Constraining ξ d , on the other hand, requires measuring, with high S/N, emission lines of refractory and non-refractory elements, such as oxygen ([O ii]λλ3726,3729, [O iii]λλ4959,5007), and nitrogen and/or sulphur ([N ii]λλ6548,6584, [S ii]λλ6716,6731, e.g. see figure 3 of Gutkin et al. 2016). The emission lines of non-refractory elements, however, are too weak to be detected in the galaxies studied in this work. We note that although τ sfr , µ and ξ d are unconstrained by the simulated data, their inclusion in the fitting allows us to appropriately propagate into the constraints on the other physical parameters the uncertainty deriving from their poor knowledge, i.e. τ sfr , µ and ξ d in this analysis play the role of 'nuisance' parameters.

We now discuss the effect on the precision of the physical parameters retrieval of the Hα line moving beyond the spectral coverage of NIRSpec at z 7 (see Fig. 8). Table 4 shows that the error threshold s 0.5 θ for the parameters log ψ and 12 + log(O/H) remains constant or slowly increases from the B, to V, to I dropouts, reflecting the presence of fainter objects, on average, when moving from lower to higher redshift sources (see Fig. 1). The Z and Y dropouts, on the other hand, exhibit significantly larger s 0.5 θ values, likely because of the disappearance of the Hα line from the NIR-Spec data at those redshifts. These parameters are in fact mainly constrained by Balmer lines (ψ) or line ratios involving such lines [12 + log(O/H)], and given that Hα is intrinsically more luminous than Hβ, and that it suffers less attenuation, the constraining power of NIRSpec data suddenly drops as Hα is shifted outside the NIR-Spec range. The constraints on the ionization parameter log U S , on the other hand, suffer less from the disappearance of Hα, since the [O ii]λλ3726,3729 line also provides constraints on this parameter. The NIRSpec constraints on the stellar mass (or mass-to-light ratio) do not significantly evolve with redshift, since most galax-ies in our catalogue exhibit low S/N in the continuum at all redshifts. A similar behaviour is observed for τV , likely because of the small absolute values of this parameter, and comparatively large retrieval errors, in our simulations. The age of onset of star formation log(t/yr) is the parameter with the poorest constraints: as for τ sfr , constraining the galaxy age requires high S/N detections of the stellar continuum of evolved stellar populations, which is only possible for the brightest galaxies in the sample.

We note that the constrains on galaxy physical parameters discussed in this section do not account for model systematic uncertainties, since we have adopted the same physical model to create the mock spectra used as inputs for the NIRSpec simulations and to fit the simulated spectra. Physical models describing galaxy emission (e.g. stellar population and photoionization models) are built from the combination of several (uncertain) ingredients and assumptions, whose effect on the model predictions is difficult to quantify (e.g. see Conroy 2013, and references therein). The constraints reported in Figure 13 and Table 4 hence must be considered as lower limits on the precision with which fundamental galaxy properties can be measured with deep JWST/NIRSpec observations.

DISCUSSION

The results presented in the previous section demonstrate how the unique combination of sensitivity and wavelength coverage of the NIRSpec/prism configuration will allow the precise characterisation of the physical properties of star-forming galaxies across a wide range of redshifts and luminosities. Specifically, the ability to observe multiple optical emission lines at high S/N should enable precise measurements of the gas-phase metallicity and starformation rate, which, combined with stellar mass constraints from NIRCam, should significantly boost our ability to understand the mass and metal assembly in galaxies at redshift z ∼ 3 -10. We note that a way to improve the precision of the physical parameter estimates of objects with low S/N emission line detections would be to combine them with galaxies with higher S/N measurements into a Bayesian hierarchical framework. This would enable populationwide constraints on the parameters and to increase the precision of individual constraints. An ongoing effort to study such an approach is currently underway (Curtis-Lake et al., in prep).

The analysis presented so far illustrates the statistical constraints on several galaxy properties that can be obtained from 100 ks NIRSpec/prism observations of high-redshift sources, but it does not provide any insight on the actual number of sources that will be observed with a NIRSpec pointing. Here, we therefore discuss the number of HST-selected galaxies expected to be observed in a single NIRSpec/prism pointing centred on the HUDF. Fig. 14 shows the HUDF/GOODS-South area along with the four NIR-Spec MSA quadrants, where we centred one quadrant on the 4.6 arcmin 2 area of the HUDF with deep WFC3 observations. Fig. 14 shows how HUDF sources will only be observable from one MSA quadrant, while for the remaining three quadrants galaxies will have to be selected from shallower data in GOODS-South. 6 We report in Table 5 the expected number n p of galaxies at z ∼ 4 -5 a For XDF, we include the B and V dropouts. For CANDELS, we consider galaxies at z ∼ 4 and z ∼ 5 from table 4 of Bouwens et al. (2015a). b For XDF, we include the I, Z and Y dropouts. For CANDELS, we consider galaxies at z 6 from the same table as in a . c On-sky target density of the objects in units of arcmin -2 . d Average number of non-overlapping spectra observed with the NIR-Spec/prism in a single pointing. e The number of objects N obj refers to the entire XDF area of ∼ 4.7 arcmin 2 , while to compute n p we consider 1 MSA quadrant. f We consider the total number of objects N obj in CANDELS/GOODS-South 'deep', corresponding to an area of ∼ 64.5 arcmin 2 and z 6 that can be observed with a single NIRSpec pointing in the HUDF/GOODS-South area, along with the on-sky density of galaxies in the HUDF/XDF and GOODS-South regions. Because the faintest targets in the HUDF/XDF may not be detected at a significant level in our deep NIRSpec/prism observations, we also report the expected number of targets observed in the HUDF/XDF with m F160W < 29.5 and < 29. From the on-sky density of the targets, and assuming Poisson-distributed sources, a 3 shutters slitlet configuration, and allowing targets to occupy any position within the open area of a shutter (Ferruit at al., in prep), we can compute the number of expected sources observed in a single NIR-Spec pointing. At z 6, we consider ∼ 200 potential targets in the HUDF/XDF and ∼ 300 in CANDELS/GOODS-South, which implies ∼ 30 non-overlapping galaxy spectra observed, on average, in a single pointing. At z ∼ 4 -6, the larger number of targets, ∼ 500 in the XDF and ∼ 2000 in CANDELS/GOODS-South, translates into ∼ 60 spectra observed in a single pointing. This will leave ∼ 100 unused slitlets which will be occupied by galaxies at lower redshift. Table 5 also indicates that restricting the z 6 targets in the HUDF/XDF to ∼ 80 galaxies with m F160W 29 (which should exhibit S/N Hβ > 5, see Table 2), reduces by 1/3 (from ∼ 30 to ∼ 20) the number of observed galaxies in a pointing. Given the paucity of z 6 targets, it is therefore advisable to prioritise the brighter targets without excluding the fainter ones, as we have previously shown that a significant fraction of m F160W > 29 galaxies should show S/N > 3 detections of Hβ (see Table 2). We note that the main factor limiting the NIRSpec multiplexing of z 6 galaxies is their on-sky density outside the HUDF/XDF region, as this density is only ∼ 5 arcmin -2 in the CANDELS/GOODS-South 'deep' area, where data reach a 5 σ depth of m F160W ∼ 27.5.

Our previous analysis has shown that 'deep' NIRSpec/prism observations will be able to constrains fundamental properties of z 6 galaxies, such as their gas-phase metallicities and star formation rates (Sec 4). The above discussion indicates, however, that obtaining a statistically significant sample of objects with precise measurements of these properties will require several 'deep', 100 ks, pointings centred on the deepest HST/WFC3 fields, later exploiting targets selected from JWST/NIRCam imaging. Figure 14. The four quadrants of the NIRSpec micro-shutter array (MSA) overplotted on galaxies at redshift 4 z 6 (blue symbols) and z 6 (orange symbols) selected by Bouwens et al. (2015a) in the HUDF/XDF (circles) and CANDELS/GOODS-South 'deep' (stars) regions.

SUMMARY AND CONCLUSIONS

In this work, we have presented statistical constraints on different galaxy physical properties obtained from the analysis of simulated NIRSpec/prism observations of high-redshift sources. We have introduced an original approach to build a semi-empirical 'mock' catalogue of z 3 galaxy spectra, based on observed photometric SEDs of a sample of ∼ 700 galaxies identified as B-, V-, I-, Z-and Y-band dropouts in the Hubble Ultra Deep Field. We used the stateof-the-art beagle tool to match the observed UV-to-near infrared broad-band SED of each galaxy with a diversity of model spectra statistically consistent with the HST photometry. We showed how additionally requiring that the models follow a (redshift-dependent) relation between stellar mass and star formation rate reduces degeneracies among model parameters, allowing us to relate, in a statistical way, model spectra with HST-selected galaxies. We validated the above approach by comparing the predicted Spitzer/IRAC band-1 (3.6 µm, IR 1 ) and band-2 (4.5 µm, IR 2 ) fluxes of the resulting mock catalogue of galaxy spectra with observed stacked photometry of galaxies at z ∼ 4 and z ∼ 5, and with the evolution of the IRAC colour IR 1 -IR 2 at redshift z ∼ 2 -7. We also showed that our model predictions cover the entire observed range of [C iii]λ1907+C iii]λ1909 equivalent widths at z ∼ 1 -8.

We have exploited the mock galaxy catalogue described above to simulate JWST/NIRSpec spectra using an 'Exposure Time Calculator'-like (ETC-like) approach. Specifically, we considered NIRSpec observations performed with the low-resolution prism, which covers the wavelength range λ ∼ 0.6 -5.3 µm, assuming a total integration time of ∼ 100 ks. Such an observational strategy is tailored to the detection of faint, distant galaxies out to z ∼ 10. Our simulations show that this strategy enables the detection of several optical emission lines (e.g. Hβ, Hα, [O iii]λλ4959,5007) for the majority of HST-selected targets out to z 8. In particular, we expect a S/N > 3 detection of the Hβ line in most z 3 galaxies brighter than m F160W ∼ 30, and a significant fraction of detections of Hβ in m F160W > 30 objects. Our calculations also imply S/N 20 de-tections of the Hα and [O iii]λλ4959,5007 lines for most sources in the HUDF.

We also used the beagle tool to analyse the simulated NIR-Spec spectra. We could show in this way that the simultaneous detection of several optical emission lines enables tight constraints on the star formation rate, dust attenuation and ionised-gas properties in z 3 galaxies. Our analysis indicates that such observations will allow measurements of star formation rates, ionization parameters and gas-phase metallicities within a factor of ∼ 1.5, mass-to-light ratios within a factor of ∼ 2, galaxy ages within a factor of ∼ 3 and V-band attenuation optical depths with a precision of ∼ 0.3. We also showed how the shift of the Hα line outside the NIRSpec wavelength range at z 7 lowers the precision of the constraints on most physical parameters, especially star formation rate, gas-phase metallicity and ionization parameter. Finally, we have discussed our findings in the context of a NIRSpec/prism pointing centred on the HUDF. We showed that the small area covered by deep HST/WFC3 observations (4.7 arcmin 2 ) and the low on-sky density of high-redshift targets in the CANDELS/GOODS-South region limit the expected number of z 6 galaxies that can be observed in a single NIRSpec pointing. We have shown that, on average, a single NIRSpec pointing, covering an effective (non-contiguous) area of ∼ 9.8 arcmin 2 , will enable the observation of ∼ 20 (∼ 15) galaxies at z 6 (z ∼ 4 -6) selected from the HUDF/XDF region, and ∼ 10 (∼ 45) galaxies at z 6 (z ∼ 4 -6) from the shallower CANDELS/GOODS-South one. This therefore indicates that multiple 'deep' NIRSpec/prism pointings targeting known high-redshift galaxies will be necessary to assemble a statistically significant sample of z 4 galaxy spectra with high-S/N emission-line detections.

In the future, we plan to extend our simulations to account for different observational strategies, in particular considering the combination of low-, medium-and high-spectral resolution configurations used in the NIRCam and NIR-Spec Guaranteed Time Observer teams for their galaxy assembly survey (Rieke et al. 2017). The simulations presented in this work will be made public on the Github repository https://github.com/jacopo-chevallard/Simulating-deep-NIRSpec-observations-in-the-Hubble-Ultra-Deep-Field. cosmic rays, as this is expected to be a minor effect for the short individual exposures that we have adopted.

A2 Extraction of the spectrum

As we simulate integrated, 1-dimensional observed spectra, we must also consider a spectrum extraction strategy. We assume a simple extraction strategy by summing over n spa = 4 detector pixels along the spatial direction, ignoring any light that may be lost in the process. More accurate simulations, based on the NIRSpec Instrument Performance Simulator, show that this assumption is correct within a few percent for a point source. No summation was performed along the spectral direction (n spe = 1).

A3 Electron rates

In order to estimate the noise variance associated to each element of the simulated spectrum, we need to compute the number of electrons accumulated during an exposure, for both the object and the background (zodiacal light + parasitic light from the telescope) spectra. For the object, we start from a high-resolution, spatially integrated, redshifted spectrum F o (λ) in units of W m -2 m -1 (spectral irradiance) generated with the beagle tool as detailed in Section 2.

By linearly rebinning the spectrum, we construct an input spectrum R o (k) in units of Wm -2 on the pixel grid (n k = 365 typically) representative of the NIRSpec PRISM/CLEAR configuration. In this process, we take into account the highly variable spectral resolution of the NIRSpec/prism (R ∼ 30 to ∼ 300), and consider the wavelength range 0.6 -5.0 µm. Note that there is no one-to-one correspondence between the pixels of the rebinned input spectrum and the actual detector pixels, as the former contain the spatially integrated irradiance of the object, i.e. after summation over 4 detector pixels along the spatial direction.

Note that in this approach we do not convolve the input spectrum by the (wavelength-dependent) line-spread function of the instrument. As the emission lines in the input spectrum are spectrally unresolved at the NIRSpec/prism spectral resolution, this causes them to be unrealistically concentrated into single pixels in our rebinned spectra. This yields, on average, a ∼ 35 per cent larger emission line S/N in our simulated spectra with respect to a computation accounting for the spectral response. As highlighted in section A6, this overestimation of the peak signal-to-noise ratio of the emission lines is largely compensated by the conservative values adopted in other steps of our computations (see Table A1).

The number of electrons N o (k) accumulated during an individual exposure and associated to the object is computed using the relation

N o (k) = R o (k) × t final × A telescope × PCE(λ k ) × η o (λ k ) E(λ k ) , (
(k) = F b (λ k ) × δλ(k) × t final × A telescope × PCE(λ k ) × Ω sky E(λ k ) , ( A2
) Where F b (λ k ) is the spectral radiance (surface brightness) of the background in units of W m -2 m -1 arcsec -2 , δλ(k) the size of the current spectrum pixel in units of wavelength (in m) and Ω sky the projected on-sky area of a background element (in arcsec 2 ). The area Ω sky is computed considering, along the spectral direction, the width of a micro-shutter (0.20 arcsec), and 4 detector pixels (0.1 arcsec each) along the spatial direction. For the surface brightness of the zodiacal light, we use a simple model which combines a power law and a 256 K black-body spectrum. We then scale this zodiacal light model to match the typical spectral radiance at the north ecliptic pole, then multiply it by a factor 1.2, as customary done to evaluate the NIRSpec performances. For the contribution of the telescope straylight and background, we use constant values of 0.091 MJy sr -1 below 2 µm and 0.07 Mjy sr -1 above 3 µm, with a linear interpolation in between. These values correspond to the requirements placed on the observatory at 2 and 3 µm, and are conservative levels when compared with the expected background in Hubble Deep Field South region (e.g. [START_REF] Wei | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF].

A4 Noise model

We use a simplified noise noise model in which the total spectrum variance is computed considering two noise sources:

• shot noise of the object, of the background and of the detector dark current (assuming a conservative dark current rate of d = 0.01 electron s -1 pixel -1 );

• detector noise (σ r ), assuming a conservative value of 7 electrons pixel -1 .

The total variance of an individual exposure V(k) is then computed using the relation

V(k) = N o (k)+ 1 + 1 n b × N b (k) + n spa n spe × d t eff + σ 2 r , (A3)
where N o (k) is the number of accumulated 'object' electrons perexposure, n b the number of background spectra available perexposure, N b (k) the number of accumulated 'background' electrons, n spa and n spe the number of detector pixels along the spatial and spectral directions used in the spectrum extraction, d the dark current rate, t eff the effective integration time and σ r the detector noise rate.

A5 Generating the simulated spectrum

After computing the variance, we can obtain the simulated, noisy spectrum S o (k) (in W m -2 m -1 ) assuming a perfect background subtraction (i.e. no residuals other than the noise):

S o (k) = R o (k) δλ(k) + α 2 (k) × V(k) n exp × N k (0, 1) , ( A4
)
where R o (k) is the rebinned (noiseless) spectrum, δλ(k) the size of the spectrum pixel, α(k) the conversion factor from electrons to W m -2 m -1 , V(k) the total variance, n exp the number of exposures and N k (0, 1) a random Gaussian variate drawn from a distribution with zero mean and unit standard deviation. Note that we have also considered additional noise terms not shown in equation (A3) to mimic the effects of a non-ideal detector dark-current subtraction (10 per cent accuracy, 1σ per pixel, no correlation from one pixel to the next) and of a pixel-to-pixel flatfield correction (20 per cent accuracy, 1σ per pixel, no correlation from one pixel to the next). These two contributions have, however, very little impact on the simulated spectra for the dark current and signal-to-noise ratio levels considered in this paper.

A6 Caveats and limitations

The present work is based on simulations of NIRSpec/prism observations computed with a simple ETC-like approach. The flexibility and versatility of this approach allows us to easily explore different observational strategies (e.g. combination of spectral configurations, depth of observations) but it is based on some simplifying assumptions about the NIRSpec instrument. In Table A1, we list various effects that can affect the signal-to-noise ratio and calibration of NIRSpec spectra, and indicate whether these effects are accounted for in our simulations and in the official JWST ETC [START_REF] Pontoppidan | Observatory Operations: Strategies, Processes, and Systems VI[END_REF]. Table A1 shows that, depending on the effect considered, our simulations can be optimistic or pessimistic.

Accounting for these effects, however, would require significantly more complex and computationally-intensive simulations, for in-stance based on the NIRSpec Instrument Performance Simulator, and more realistic data reduction strategies.

APPENDIX B: SLIT-LOSSES FROM THE NIRSPEC MICRO SHUTTER ARRAY (MSA)

A standard NIRSpec MSA mask will consist of a series of 1×3 shutter 'microslits', with the object of interest placed at the center of each shutter as part of a standard nodding scheme. However, the extended size of an object compared to the small (0.20 × 0.46 ) shutters, the bars between shutters, difficulties in perfectly centring every object within the shutters, and the large PSF at long wavelengths ensures that not all of the incident flux from a galaxy will make it onto the detector. In order to quantify the effect of these sources of flux losses, we create simulations based on models of the instrument and models of galaxy light distributions. In order to average over all source position angles and ellipticities, throughout we assume a round source with a Sérsic light profile. For a grid of Sérsic indices and half-light radii, an object is positioned in multiple locations inside the 1×3 microslit and the losses from both geometry and diffraction are calculated as a function of wavelength within the full aperture. These calculations include all of the changes to the PSF as the light propagates to the MSA plane, including the Optical Telescope Element (OTE) and the NIRSpec foreoptics.

Since the specifics of the mask layout and the exact positioning of objects within slits is not known, we calculate the median throughput as a function of wavelength for an object with a given size and Sérsic index, adopting the APT centring tolerances, e.g. accepting that the object before nodding is within the full pitch of the central shutter (including the bars between shutters), accepting that the object is in the open area, etc. We include the effect of nodding within the 1×3 microslit aperture in these calculations, where more light is lost at the upper and lower nodding positions. Since the quantity of interest is the relative throughput with respect to a centred point source, we correct our estimates according to the wavelength-dependent transmission of a centred point source as calculated using the same framework.

This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 3 .

 3 Figure 3. Stacked two-dimensional posterior probability distribution of mass and star formation rate for the B, V, I, Z and Y dropouts (grey density contours, plotted on a logarithmic scale). Circles, colour coded as in Fig. 1, indicate the pairs [M * , ψ] corresponding to a random realisation of model parameters drawn from the re-weighted posterior probability distribution. The solid yellow line indicates the mass-star formation rate relation measured by Santini et al. (2017) at z ∼ 5 -6 from ∼ 50 galaxies selected from 4 gravitationally-lensed HST Frontier Fields, while the dashed line is an extrapolation of this relation at lower stellar masses.

Figure 4 .

 4 Figure 4. Distribution of model free parameters of the mock catalogue of galaxy spectra computed as detailed in Sections 2.2-2.3. The colour coding indicates mock galaxies at different redshifts, and it is the same as in Fig. 1. The off-diagonal panels show the relation between each pair of parameter [log(M/M ), ψ, log U S , 12 + log(O/H), ξ d , τV ], while the diagonal panels display, by means of histograms, the distribution of each parameter.

Figure 6 .

 6 Figure 6. Comparison of the Spitzer/IRAC colour 3.6µm -4.5µm at 2 z 8 from the literature (coloured circles) with that of the galaxies in our mock catalogue (grey circles). In order of increasing redshift, we show galaxies with spectroscopic redshifts at 2 z spec 4 from 3D-HST/GOODS-North and GOODS-South (cyan circles, Skelton et al. 2014), objects from the spectroscopic sample of Smit et al. (2016) (light green) and from the photometric samples of Rasappu et al. (2016) (dark green), Smit et al. (2015) (orange) and Smit et al. (2014) (red). We only plot objects for which the quoted uncertainty on the IRAC colour is 0.30. For clarity, we only show error bars when these are 0.05. As indicated in the inset legend, the area of the circles is proportional to the absolute UV magnitude of each galaxy M uv .

Figure 7 .

 7 Figure 7. Left: comparison of the [C iii]λ1907+C iii]λ1909 (rest-frame) equivalent widths W C iii] from literature (coloured circles) with those measured from the spectra of our mock catalogue (grey circles). Literature values are taken, in order of increasing redshift, from Stark et al. (2014, cyan circles), Maseda et al. (2017, light green), Amorín et al. (2017, blue), Ding et al. (2017, dark green), Stark et al. (2015, orange) and Stark et al.(2017, red). The upper limits indicate 3 σ limits. As in Fig.6, the area of the circles is proportional to the absolute UV magnitude M uv of each galaxy. Right: distribution of [C iii]λ1907+C iii]λ1909 equivalent widths in our mock catalogue and in the sample of z ∼ 2 -3.8 star-forming galaxies of LeFèvre et al. (2017). As LeFèvre et al. (2017) consider only galaxies with W Ciii] /Å 3, we renormalize the distribution of LeFèvre et al. (2017) to the total number of objects in our sample with W Ciii] /Å 3.

  ii λ1640 + O iii] λλ1660, 1666 Si iii] λλ1883, 1892 + C iii] λλ1907, 1909 O ii] λλ3727, 3729 Ne ii] λ3868 Hδ Hγ + [O iii] λ4363 Hγ [O iii] λ4363 Hβ Hβ + [O iii] λ4959

Figure 8 .

 8 Figure 8. Visibility and blending/unblending of the strongest UV and optical emission lines as a function of redshift, for NIRSpec/prism observations. The thin grey lines indicate the lines visibility as a function of redshift, computed assuming the nominal NIRSpec wavelength coverage λ = 0.6 -5.3 µm. The thick black lines show the unblended single lines or blended group of lines which will be observable at different redshifts.

  scheme used by JWST's official ETC 'Pandeia'[START_REF] Pontoppidan | Observatory Operations: Strategies, Processes, and Systems VI[END_REF].In our simulations, we also account for the fact that a significant fraction of the object light may fall outside of the standard 3shutter slitlet, especially for extended objects (the projected size of the aperture of an individual micro-shutter is 0.20 × 0.42 arcsec). In practice, we compute the (wavelength-dependent) throughput of an extended source with effective radius r e , described as an exponential profile with Sérsic index 1.5(Shibuya et al. 2015), and average over all source positions within the aperture of the central shutter of the slitlet. Averaging over random ellipticities and position angles produces small corrections ( few per cent), which we therefore ignore. Following the approach ofWilliams et al. (2018) (see their section 5.2), we associate an effective radius to each galaxy in the catalogue by adopting the (redshift-dependent) relation between galaxy UV luminosity and galaxy size ofShibuya et al. (2015) r e = r 0 e L uv L 0 uv

Figure 9 .

 9 Figure 9. (a) Simulated NIRSpec/prism spectrum (thick grey line, corresponding to a ∼ 100 ks exposure) of a V dropout (XDFV-3948262597, m F160W = 27.70) placed at redshift z = 5.072, with stellar mass log(M * /M ) ∼ 8.18, star formation rate ψ/M yr -1 ∼ 0.7 and specific star formation rate log (ψ s /yr -1 ) ∼ -8.4.

  Fig. 8 indicates that groups of emission lines in the UV such as He iiλ1640 + O iii]λλ1660,1666 and Si iii]λλ1883,1892 + [C iii]λ1907+C iii]λ1909 are blended at all redshifts, while at z 4 Hβ is unblended with respect to [O iii]λ4959 and [O iii]λ5007, and these [O iii] lines are unblended at z 5. The Hα line is blended with [N ii]λλ6548,6584 at 2 z 6.5, while Hγ is blended with [O iii]λ4636 up to z ∼ 8. Note also that all the lines in Fig. 8 but [C iii]λ1907 + C iii]λ1909 and [O ii]λ3726 + [O ii]λ3729 are unblended when using the R ∼ 1000 disperser, while the R ∼ 2700 disperser allows one to resolve also those doublets at most redshifts. Fig. 8 also highlights the ability of NIRSpec/prism observations to provide the simultaneous measurement of the main UV and optical emission lines for galaxies at z ∼ 3 -10: Hα is accessible up to z ∼ 7, Hβ to z ∼ 10 and [O iii]λλ4959,5007 to z ∼ 9.5. We show in Fig. 9 and 10 a few examples of simulated NIRSpec/prism spectra. Fig. 9(a) shows a z = 5.072 galaxy with m F160W = 27.70, and Fig. 9(b) one at z = 7.176, with m F160W = 28.74. Fig. 9(a) shows a high signal-to-noise detection ( S/N Hβ ∼ 8) of the Balmer lines Hα and Hβ, and of the Oxygen lines [O iii]λ4959 and [O iii]λ5007. Fig. 9(b) illustrates that the sensitivity of JWST/NIRSpec will allow us to obtain highly significant detections ( S/N Hβ ∼ 10) of emission lines (Hγ, Hβ, [O iii]λ4959 and [O iii]λ5007) even for a fainter galaxy at redshift z ∼ 7. The inset of Fig. 9(b) shows a flux excess in the HST F140W band, but no strong emission line that can cause such an excess falls in the F140W band at z ∼ 7, making the origin of the observed flux excess unclear. Fig.

Figure 10 .Figure 11 .

 1011 Figure 10. Illustration of simulated NIRSpec/prism spectra (corresponding to a ∼ 100 ks exposure) of 10 galaxies at z ∼ 4 -8. The grey line indicates the noiseless input spectrum, while the blue line shows the simulated one. Each spectrum is normalized to its maximum value, then shifted vertically for clarity. The magnitudes reported in the figure are observed ones and refer to the HST/WFC3 F160W filter, while the redshift is the beagle-based photo-z of the mock galaxy. Shaded regions indicate the location of the main optical emission lines [O ii] (blue), Hβ + [O iii] (green) and Hα + [N ii] + [S ii] (red).

Figure 12 .

 12 Figure 12. Relation between S/N of the Hβ emission line and observed F160W magnitude. The colour coding is the same as in Fig 1 and reflects galaxies in different dropout classes.

[

  O iii]λ5007. The factors λ Hα /λ Hβ and λ [O iii] /λ Hβ account for the different number of photons, and hence of photo-electrons, generated at different wavelengths at fixed emitted energy.

  two examples of the beagle fitting of the simulated spectra. The red lines indicate the spectra obtained by computing, in each pixel, the median of the posterior distribution of predicted fluxes, while the red shaded regions show the 95 per

Figure 13 .

 13 Figure13. Comparison between the input physical parameter (on the x-axis) and the value retrieved by fitting the full NIRSpec/prism simulated spectra with beagle (on the y-axis) for the 6 model parameters stellar mass log(M/M ), log(t/yr), current star formation rate log(ψ/M yr -1 ), ionization parameter log U S , gas-phase metallicity 12 + log(O/H) and attenuation optical depth τV . The different colours correspond to the different dropout classes, and follow the colour coding of Fig.1. The retrieved parameter value corresponds to the posterior median of the marginal posterior distribution, while the error bars to the 68 per cent central credible region. To avoid overcrowding the plot, we only show errorbars for 20 randomly selected objects.

  A1)where λ k is the central wavelength of the k th spectrum pixel (in m), A telescope collecting area of the telescope (25 m 2 ), PCE(λ k ) the photon conversion efficiency of the telescope + instrument (in units of electrons-per-collected-photon), η o (λ k ) the slit losses factor (unitless, see Appendix B) and E(λ k ) = hc/λ k the energy-per-photon (in J photon -1 ). The photon conversion efficiency used in these simulations is consistent with the one used in the official JWST ETC version 1.1.1, available in October 2017.The number of accumulated electrons N b (k) for the background (zodiacal light + parasitic telescope light) is computed us-ing the relation N b

Table 1 .

 1 Priors relative to the different free parameters used in the beagle tool for the broad-band SED fitting of HST/XDF photometry and for the fullspectrum fitting of JWST/NIRSpec simulated spectra. The symbol N(mean; sigma) indicates a Gaussian (Normal) distribution.

	Parameter	Prior	Description	XDF photometry NIRSpec spectra
	z	Uniform ∈ [0, 15]	Redshift
	log(M/M )	Uniform ∈ [5, 12]	Stellar mass (does not account for fraction of mass returned to the ISM by stellar mass loss)
	log(τ sfr /yr)	Uniform ∈ [7, 10.5]	Time scale of star formation in a SFH described by a delayed exponential function.
	log(t/yr)	Gaussian N(8.5; 0.5) truncated ∈ [7, t max form (z)]	Age of onset of star formation in the galaxy
	log(Z/Z )	Uniform ∈ [-2.2, 0.25]	Stellar metallicity
	log(Z ism /Z )	Uniform ∈ [-2.2, 0.25]	Interstellar metallicity
	log(ψ/M yr -1 )	Gaussian N(0; 2) truncated ∈ [-4, 4]	Star formation rate averaged over the last 10 7 yr
	log U S	Uniform ∈ [-4, -1]	Effective gas ionization parameter
	ξ d	Uniform ∈ [0.1, 0.5]	Dust-to-metal mass ratio
	τV	Exponential exp(-τV ) truncated ∈ [0, 5]	V-band attenuation optical depth
	µ	Uniform ∈ [0, 1]	Fraction of V-band attenuation arising in the diffuse ISM
	eters Θ as

  of a correlation between mass and metallicity, and metallicity and ionization parameter. The mass-metallicity relation is well-established at low redshift (e.g. from SDSS data, Tremonti

		24	F775W ∼ 27, N obj = 147 (55) F775W ∼ 28, N obj = 62 (79) F775W ∼ 26, N obj = 108 (25)		Model SEDs Gonzalez+2012	(a)		24	F850LP ∼ 28, N obj = 15 (31) F850LP ∼ 27, N obj = 50 (19)		Model SEDs Gonzalez+2012	(b)
	AB mag	26						AB mag	26
		28							28
	-0.5 0.0 0.5 -f ( f G model )/ f G	3.8	4.0	4.2	4.4	4.6	( f G -f -0.5 0.0 0.5 model )/ f	3.8	4.0	4.2	4.4	4.6
			3.8	4.0	4.2	4.4	4.6			3.8	4.0	4.2	4.4	4.6
				log(λ eff /Å)						log(λ eff /Å)

G Figure 5. (a) In the top panel, diamonds with (1 σ) error bars indicate the stacked photometric SEDs of redshift ∼ 4 galaxies (B dropout) of

González et al. (2012)

in three bins of observed F775W magnitudes, used as a proxy for the rest-frame UV luminosity. Shaded 'violins' indicate the magnitude distribution (95 per cent interval) covered by our mock spectra based on B dropouts in the XDF, split in the same UV luminosity bins as in

González et al. (2012)

. We indicate in the inset legend the number of objects in each

González et al. (2012)

stack, and in parenthesis the number of objects in our mock catalogue entering each bin. In the bottom panel we show the distribution of the flux differences between the

González et al. (

Table 2 .

 2 Fraction of objects in different bins of observed F160W magnitude with S/N Hβ > 3, 5 and 10 (in parentheses the number of objects in each magnitude bin).

	m F160W	S/N 3	S/N 5	S/N 10
	m 28 (199)	0.51	0.42	0.29
	28.0 < m 28.5 (78)	0.65	0.49	0.17
	28.5 < m 29.0 (104)	0.49	0.33	0.07
	29.0 < m 29.5 (160)	0.39	0.14	0.02
	29.5 < m 30 (124)	0.40	0.13	0.02
	m > 30 (44)	0.27	0.16	0.00

Table 3 .

 3 Fraction of objects in different bins of redshift and observed F160W magnitude for which the deep NIRSpec/PRISM observations considered in this work would allow the detection of at least 1 or 2 emission

	lines with S/N	3 and S/N	5. The number of objects in each bin is
	given in parentheses.
	Redshift	m F160W	S/N 3	S/N 5
			N 1	N 2	N 1	N 2
		m 28 (147)	0.98	0.95	0.95	0.90
	z 5	28 < m 29 (108) 29 < m 30 (158)	0.99 0.96	0.94 0.84	0.93 0.91	0.80 0.68
		m > 30 (19)	0.84	0.68	0.68	0.26
		m 28 (45)	0.96	0.89	0.87	0.87
	5 < z 7	28 < m 29 (62) 29 < m 30 (98)	1.00 0.95	0.97 0.87	0.93 0.83	0.87 0.66
		m > 30 (20)	0.95	0.75	0.70	0.40
		m 28 (7)	1.00	1.00	1.00	0.70
	z > 7	28 < m 29 (12) 29 < m 30 (28)	0.83 0.89	0.75 0.68	0.75 0.71	0.42 0.32
		m > 30 (5)	0.60	0.40	0.20	0.00

Table 5 .

 5 Number of galaxies at z ∼ 4 -6 and z 6 expected to be observed with a single NIRSpec/prism pointing centred on the HUDF.

	Catalog	m F160W		z ∼ 4 -6 a				z 6 b
			N obj	σ obj	c	n p	d	N obj	σ obj	c	n p	d
		< 29	304	∼ 65	∼ 18	80	∼ 17	∼ 11
	XDF e	< 29.5	413	∼ 88	∼ 17	132	∼ 28	∼ 16
		all	513	∼ 110	∼ 15	202	∼ 43	∼ 21
	GOODS-S f	all	2061	∼ 32	∼ 46	303	∼ 4.7	∼ 10

Table A1 .

 A1 Summary of the main instrumental effects affecting simulated NIRSpec spectra computed with the approach adopted in this work, and with the official JWST/NIRSpec ETC.

	Effect	This paper	JWST ETC	Impact on our simulations
	MSA slit-losses	Averaged	Yes	See Appendix B
				Significant impact on the peak
	Line spread function	No	Yes	signal-to-noise ratio of emission lines (optimistic, average effect of 35 per cent,
				linear impact)
		Included in the noise model;	Included in the noise model; adjusts	Limited impact since background shot
	Background	conservative level; perfect	level to on-sky position and	noise is only dominant in a small
		subtraction	observation time; perfect subtraction	wavelength range around 1.2 µm
	Detector noise model	Conservative total noise value (7 instead of nominal 6 electrons)	According to Rauscher et al. (2007) formalism	Significant impact on sensitivity limits (pessimistic, linear impact at long wavelengths)
	Noise correlation	No	No (capability not activated for NIRSpec)	Impact not evaluated. Note that a readout mode (IRS 2 ) optimised to reduce correlated noise is available for NIRSpec. a
				Significant impact on our simulations
	Detector defects	Yes (conservative end-of-life values)	No	compared to ETC ones, as we consider that almost 20 per cent of the exposure time is
				lost due to bad pixels
	Cosmic rays	No	Yes	Only minor impact as we are considering relatively short exposures
	Dark current	Included in the noise model; conservative rate of 0.01 e s -1 ; non-perfect subtraction	Included in the noise model; measured (lower) rates; perfect subtraction	Limited impact as even conservative values of dark-current are very low
				Negligible impact as we assume that
	Flat fielding	Non-ideal	Non-ideal	dithering will be used to mitigate this
				problem
	Spectro-photometric calibration	Perfect	Perfect	Impact on retrieval of galaxy properties not evaluated
	Wavelength calibration	Perfect	Perfect	No impact

a https://jwst-docs.stsci.edu/display/JTI/NIRSpec+IRS2+Detector+Readout+Mode

The Lyα line is theoretically the most luminous emission line at UV and optical wavelengths, but being a resonant line it suffers from radiative transfer effects which affect its visibility and make its physical interpretation extremely challenging.MNRAS 000, 000-000 (0000)

MNRAS 000, 000-000 (0000)

The interstellar metallicity was indicated as Z gas in the original paper describing the beagle tool(Chevallard & Charlot 2016), while here we follow the nomenclature adopted inGutkin et al. (2016).

Our definition of log U S implies a volume-averaged ionization parameter U = 9/4 U S (see equation 1 of Hirschmann et al. 2017).MNRAS 000, 000-000 (0000)

We note that the adopted prior for τV would correspond to the noninformative, scale-invariant Jeffreys prior for a one-dimensional Gaussian likelihood function depending only on τV .

The HUDF/XDF area is larger than the on-sky area of a single MSA quadrant (1.58 × 1.40 arcmin, excluding vignetted shutters), but the separation between quadrants (0.38 arcmin along one axis, 0.62 arcmin along the other) and unknown telescope roll-angle will in practice limit the possibility of selecting HUDF/XDF targets in more than one quadrant.MNRAS 000, 000-000 (0000)

ACKNOWLEDGEMENTS

The authors thank the anonymous referee for a careful reading of the manuscript which significantly improved the final version of the paper. JC, ECL and S. Charlot acknowledge support from the European Research Council (ERC) via an Advanced Grant under grant agreement no. 321323-NEOGAL. S. Carniani acknowledges financial support from the UK Science and Technology Facilities Council (STFC). SA is funded by the Spanish Ministerio de Economía y Competitividad (MINECO) under grant ESP2015-68964-P. RS acknowledges support from a NWO Rubicon program with project number 680-50-1518. RM and RA acknowledge support from the ERC Advanced Grant no. 695671-QUENCH.

APPENDIX A: GENERATING SIMULATED NIRSPEC SPECTRA

In this Appendix, we describe the different steps of the approach that we used to simulate NIRSpec spectra.

A1 Exposures and observing strategy

In this work, we use the standard exposure parameters and observational strategy baselined to monitor the NIRSpec performances during the instrument development. For the individual exposures, we consider the MULTIACCUM readout pattern 22 × 4 listed in Table 2 of Rauscher et al. (2007), which yields an effective integration time of t eff = 934.38 s per-exposure. Each exposure lasts 955.86 s, including the time needed to reset the array at the beginning of the exposure. The simulations are performed for a total of n exp = 108 exposures, corresponding to a cumulated effective integration time of 100923.84 s, slightly in excess of 100 ks.

The building block of a standard NIRSpec MOS observation of faint objects is a minimum set of 3 nodding-pattern exposures in which the object is observed consecutively in three neighbouring micro-shutters forming a 1 × 3 slitlet. For each of the 3 object positions, two measurements of the telescope and zodiacal light backgrounds are available and can be used for an efficient and accurate pixel-level background subtraction. This implies considering the average of n b = 2 observations, over the same pixels used for the object spectrum, to estimate the contribution of the background subtraction to the total variance spectrum.

The effective integration time is further reduced by a factor γ bad pixels = 0.807 to account for bad pixels. The factor γ bad pixels is estimated by computing the probabilities of loosing one, two or three measurements in a 3-slitlet nodding pattern, assuming a fraction of bad pixels of 3.5 per cent. Accounting for bad pixels hence yields a final integration time of t final = t eff ×γ bad pixels = 754.04 s perexposure. We note that the fraction of bad pixels here considered is conservative, as it corresponds to end-of-life values. The official JWST exposure time calculator ('Pandeia' ETC; [START_REF] Pontoppidan | Observatory Operations: Strategies, Processes, and Systems VI[END_REF] does not account for this effect, and this can introduce differences among our simulations and those obtained with the JWST ETC, at fixed observational setting. We do not model the impact of