JWST NIRCam+NIRSpec: Interstellar medium and stellar populations of young galaxies with rising star formation and evolving gas reservoirs

We present an interstellar medium and stellar population analysis of three spectroscopically confirmed $z>7$ galaxies in the ERO JWST NIRCam and JWST NIRSpec data of the SMACS J0723.3-7327 cluster. We use the Bayesian spectral energy distribution (SED) fitting code \texttt{Prospector} with a flexible star-formation history (SFH), a variable dust attenuation law, and a self-consistent model of nebular emission (continuum and emission lines). Importantly, we self-consistently fit both the emission line fluxes from JWST NIRSpec and the broad-band photometry from JWST NIRCam, taking into account slit-loss effects. We find that these three $z=7.6-8.5$ galaxies ($M_{\star}\approx10^{8}~M_{\odot}$) are young with rising SFHs and mass-weighted ages of $3-4$ Myr, though we find indications for underlying older stellar populations. The inferred gas-phase metallicities broadly agree with the direct metallicity estimates from the auroral lines. The galaxy with the lowest gas-phase metallicity ($\mathrm{Z}_{\rm gas}=0.06~\mathrm{Z}_{\odot}$) has a steeply rising SFH, is very compact ($<0.2~\mathrm{kpc}$) and has a high star-formation rate surface density ($\Sigma_{\rm SFR}\approx22~\mathrm{M}_{\odot}~\mathrm{yr}^{-1}~\mathrm{kpc}^{-2}$), consistent with rapid gas accretion. The two other objects with higher gas-phase metallicity show more complex multi-component morphologies on kpc scales, indicating that their recent increase in star-formation rate is driven by mergers or internal, gravitational instabilities. We discuss effects of assuming different SFH priors or only fitting the photometric data. Our analysis highlights the strength and importance of combining JWST imaging and spectroscopy for fully assessing the nature of galaxies at the earliest epochs.


INTRODUCTION
The James Webb Space Telescope (JWST) promises to reveal a new view of galaxy formation in the early Universe. For the first time, astronomers can now combine photometry (NIRCam) and spectroscopy (NIRSpec) in the rest-frame UV/optical to observe directly the growth of stellar populations and galactic structure in the first few hundred million years of cosmic history. The discoveries enabled by JWST will allow for the formation history of stars, gas, metals, and dust to be detailed during the critical Epoch of Reionization (EoR; for a review, see Robertson 2022). Together, NIRCam and NIRSpec provide information on both broad-band stellar light and emission lines powered by Lyman Continuum (LyC) photons emitted by newly ★ E-mail: st578@cam.ac.uk formed stars. This powerful combination allows us to precisely infer the properties of early galaxies without redshift uncertainties. The focus of this paper is to use jointly the initial JWST Early Release Observations (EROs) NIRCam and NIRSpec data in the SMACS0723 cluster (Pontoppidan et al. 2022) to infer the properties of galaxies with spectroscopic redshifts at ∼ 7.6 − 8.5 and thereby gain a foothold on the physics of galaxy formation in the early Universe.
Detailed stellar population analyses of early galaxies advanced significantly during the pre-JWST era thanks to combining Hubble Space Telescope (HST) and Spitzer Space Telescope (Spitzer) 3.6 − 4.5 m data. Several studies used these datasets to constrain the stellar populations of > 8 galaxies, finding populations of young galaxies (mass-weighted stellar ages of a few tens of Myr; e.g., Stefanon et al. 2021), but also more evolved, older objects (with ages of a few hundred Myr; e.g., Hashimoto et al. 2018;Laporte et al. 2021). Significant challenges with these observations are the low signal-to-noise ratio, source blending as a result of the low spatial resolution of Spitzer, and strong emission lines (i.e., H and [OIII]) that can mimic a strong Balmer/4000Å break (Finkelstein et al. 2013;Labbé et al. 2013;Smit et al. 2015;De Barros et al. 2019;Endsley et al. 2021). Because of this, Tacchella et al. (2022c) and Whitler et al. (2023) have investigated how priors might influence the inference of the stellar ages and star-formation histories (SFHs) more generally, finding that the current HST + Spitzer data are not very constraining in age-dating the highest-galaxies. Along these lines, the ultra-violet (UV) colours have been used to study the earliest phases of chemical enrichment (e.g., Wilkins et al. 2011;Finkelstein et al. 2012;Bouwens et al. 2014;Bhatawdekar & Conselice 2021). However, the amount of attenuation, the attenuation law itself and the stellar and gas-phase metallicities remain largely unconstrained at > 7.
JWST will improve the situation tremendously by delivering larger wavelength coverage and higher spectral resolution, both with spectra and medium/narrow band imaging. This will help break the degeneracy between strong emission lines and a strong Balmer/4000Å break, which is needed to age-date these galaxies accurately Tacchella et al. 2022c), enabling not only with learning more about the physical properties of these high-galaxies during and before the EoR, but also with timing cosmic dawn by delivering accurate SFHs. Several JWST programmes will deliver such datasets (i.e. spectra and medium band imaging), including the JWST Advanced Deep Extragalactic Survey (JADES), which is a Guaranteed Time Observer (GTO) programme using about 800 hours of prime time and 800 hours of parallel time to study the formation and evolution of galaxies, combining NIRSpec, NIRCam, and MIRI data in a coordinated observing programme (Rieke et al. 2019).
In this work, we make use of the JWST ERO of the SMACS J0723.3-7327 cluster field (Pontoppidan et al. 2022). The ERO multimode observations of SMACS J0723.3-7327 include NIRCam and MIRI imaging, NIRSpec Multi-Object Spectroscopy, and NIRISS Wide-Field Slitless Spectroscopy of the cluster and surrounding field. The NIRSpec MSA configuration was designed based on a photometric catalogue constructed from the NIRCam imaging, following a workflow similar to that of future science programs such as JADES. We focus on the three galaxies that have a spectroscopic redshift > 7, namely Galaxy ID 04590, 06355 and 10612 (Table 1). These three objects have been discussed in several recent papers (e.g., Arellano-Córdova et al. 2022;Brinchmann 2022;Carnall et al. 2023;Curti et al. 2023;Katz et al. 2023;Rhoads et al. 2023;Schaerer et al. 2022;Trump et al. 2023;Trussler et al. 2022). The different studies rely on different measurement techniques, but also on different reduction and calibration pipelines, making comparisons difficult. Carnall et al. (2023) performed SED fitting of the NIRCam imaging and the existing HST imaging for these three objects, finding low stellar masses (log( ★ / ) = 7.1 − 8.2) and correspondingly young mean stellar ages of only a few Myr. Curti et al. (2023) used the [OIII]4363 auroral line to measure metallicities with the direct e method, finding Galaxy ID 04590 to be extremely metal poor (12 + log(O/H) ≈ 7), Galaxy ID 10612 to be about 1/10 solar and Galaxy ID 06355 to be about one-third solar. The latter two objects are marginally consistent with the Fundamental Metallicity Relation (FMR; Mannucci et al. 2010;Onodera et al. 2016;Curti et al. 2020;Suzuki et al. 2021;Papovich et al. 2022), while Galaxy ID 04590 deviates significantly from it, consistent with being far from the smooth equilibrium between gas flows, star formation and metal enrichment in place at later epochs. This possibly indicates that we enter at > 7 another epoch of galaxy formation, where galactic systems are out of equilibrium and evolve more stochastically. Consistent with this, other studies find high ionization, low metallicity and high pressure in the interstellar medium (ISM; Katz et al. 2023;Schaerer et al. 2022;Rhoads et al. 2023;Trump et al. 2023;Trussler et al. 2022).
In this paper, we contribute to these interesting results by performing a full spectro-photometric analysis by simultaneously fitting the NIRCam photometry and the NIRSpec spectroscopy, taking into account slit-loss effects. We confront the measured gas-phase metallicity with the inferred SFHs and morphology of these three = 7.6−8.5 galaxies, highlighting the strength and importance of combining JWST imaging and spectroscopy for fully assessing the nature of galaxies at the earliest epochs. Specifically, Section 2 presents the NIRCam and NIRSpec data processing and describes the details of the stellar population inference from the spectro-photometric modelling. We present our key results in Section 3 and conclude in Section 4. We assume the latest Planck flat ΛCDM cosmology with H 0 = 67.36, Ω = 0.3153, and Ω Λ = 0.6847 (Planck Collaboration et al. 2020) and Z = 0.0142. All quantities are corrected for magnification if not otherwise indicated (see Section 2.5). When quoting values of derived quantities, we quote the median and 16 − 84% quantile-based errors if not otherwise stated.

DATA AND METHOD
All JWST observations used in this work are taken as part of the SMACS0723 ERO (Programme #2736; Pontoppidan et al. 2022). We focus in this analysis on three galaxies that are both covered with NIRCam and NIRSpec and lie at > 7 (Tab. 1).

NIRCam imaging and photometry
We use the deep NIRCam imaging data of SMACS0723 in the F090W, F150W, F200W, F277W, F356W, and F444W filters, which cover an observed wavelength range of obs = 0.8 − 5μm. The images reach (and surpass) ≈ 29 mag. We reduce the raw level-1 data products with the public JWST pipeline (v1.9.2; https: //github.com/spacetelescope/jwst), using the latest available calibration files (CRDS_CTX=jwst_1039.pmap). An additional, local background subtraction is performed on the final mosaiced images and applied to the individual exposures, see details in Robertson et al. (2023) and Tacchella et al. (2023).
We measure object fluxes using the Bayesian model-fitting program forcepho (Johnson et al., in prep.). This program simultaneously fits multiple PSF-convolved Sérsic profiles, each having 11 parameters, to all of the background-subtracted stage 2 cal exposures for every NIRCam filter. The joint posterior distribution for the parameters of all of the profiles is sampled via Markov Chain Monte Carlo (MCMC), allowing estimation of covariances between the parameters of a single object (such as Sérsic index and half-light radius) as well as covariance between objects (due to blending), and identification of any non-Gaussian posterior distributions. The simultaneous fit to all bands and the MCMC sampling allow us to obtain self-consistent estimates of object fluxes in each band. We represent Galaxy ID 06355 and 10612 with several sub-components to account for clumpy structure obvious in the higher resolution F150W imaging, and we simultaneously model all sources within 1 to account for the contribution of extended structures. To determine total source fluxes we then combine the fluxes of the sub-components while taking into account covariances among them, which we accomplish by summing the sub-component fluxes for each posterior sample and then computing the mean and standard deviation of the summed fluxes. In Figures 1, 2  and stacked residuals for one draw from the MCMC chain. We have varied the number of sub-components used for the sources and find similar total fluxes. Due to uncertainties in photometric calibration and other factors affecting these very early data, we adopt an error floor of 10% on the photometry.
The forcepho fits to the images also provide size estimates. For the single component of Galaxy ID 04590 this is taken directly from the half-light radius of the fitted Sérsic profile. For the other two Galaxy IDs that have multiple components, we estimate a combined effective half-light radius by generating a model image with the PSF removed using the F150W fluxes, compute the barycenter, and determine the radius containing half the combined model flux. This yields half-light radii of 0.039 ± 0.004 (0.18 ± 0.02 kpc), 0.077 ± 0.008 (0.39 ± 0.04 kpc) and 0.100 ± 0.010 (0.51 ± 0.05 kpc) respectively for Galaxy ID 04590, 06355, and 10612. The size for Galaxy ID 10612 is driven by the separation between the two dominant components, which are each 0.03 in half-light radius. We note that we have not corrected these sizes for lensing distortions (Section 2.5). We expect this effect to be small for Galaxy ID 06355 and 10612, because they are only weakly magnified. This effect could however be more substantial for Galaxy ID 04590, indicating that our size is possibly an upper limit.

NIRSpec spectroscopy
NIRSpec observations for SMACS0723 were carried out using two disperser/filter combinations, but given the redshift of our targets, in this work we use only data from the G395/F290LP combination, covering the wavelength range 2.9-5.2 μm with a spectral resolution ∼ 1000. The observations consist of two pointings, with a total integration of 8,840 seconds. We show the slit position in In Figures  1, 2, and 3. For all three galaxies, the slits cover the main part of the galaxies, which supports our assumption of modelling the photometry together with the spectroscopy (i.e., treating the spectroscopic observations as representative of the galaxy). Nevertheless, slit losses are important (see below). We estimate the slit loss to be 12%, 29% and 39% for Galaxy ID 4590, 6355, and 10612, respectively, using the F356W band, which probes the wavelength range of the spectrum.
In this paper, we use the fully reduced spectra from Curti et al. (2023) 1 , which are based on level-2 data from the MAST archive and were processed using the GTO pipeline (NIRSpec/GTO collaboration, in prep.; Ferruit et al. 2022). The data and data reduction procedure are already described in Curti et al. (2023), so here we report only a summary.  The GTO pipeline uses a custom flagging and masking algorithm for bad pixels and cosmic rays, which effectively removes all artefacts from the final spectra. The extraction is performed with a custom aperture, using a boxcar extraction. After inspecting the exposures of the individual nods, it was found that one of the shutters on which Galaxy ID 04590 was nodded did not open in Obs 7 (Rawle et al. 2018); the affected nod positions were discarded before stacking the data. Stacking was performed with the GTO pipeline, taking into account the variance and quality arrays of both pointings. We refined the flux calibration using an empirical response function, derived from a spectro-photometric star observed during commissioning (Programme #1128, Gordon et al. 2022). Because our three targets are only marginally resolved, we apply the path-loss correction derived for a point-like source. The path-loss correction appropriate for an extended uniform source gives a spectrum that has higher overall flux, but approximately the same shape as our default choice 2 . As we discuss below (see Section 2.3), we introduce a nuisance parameter in the spectro-photometric modelling to account for any mismatch between the photometry and spectroscopy (related to slit loss, calibration issues as well as physical effects). Therefore, it does not matter in this analysis whether we use the point-like or extended source approximation. The resulting spectra are shown in the lower-left corner of Figures 1-3.
Emission-line fluxes were measured using (Cappellari 2017), using the configuration described in Curti et al. (2023). Briefly, we fit the stellar continuum using a linear combination of simple stellar-population spectra from the C3K library ), using the MIST isochrones (Choi et al. 2016). The continuum is scaled using a 10 th -order Legendre polynomial. Emission-lines are modelled as Gaussians. All the model spectra (both continuum and emission-line) are constrained to have the same velocity and velocity dispersion; note that the continuum is only marginally resolved; replacing the SSP template spectra with a constant does not change the emission-line fluxes (Curti et al. 2023). The resulting fluxes have already been published (Curti et al. 2023, their Table 1) and are not repeated here. We remark that fitting each line individually does not change our conclusions, as reported in Katz et al. (2023). The spectroscopic redshifts in this work are based on the emission-line velocities measured by . In the fitting, we mask the following emission lines: [OIII]4364,H ,H ,and [NeIII]3870. We mask [OIII]4364 because we cannot reproduce very high [OIII]4364 emission line fluxes (as measured here) with our current cloudy-based modelling for the nebular emission (see next section). As discussed in Dors et al. (2011), photoionization models have trouble predicting [OIII]4364 at a necessary accuracy (see also Brinchmann 2022, who has also excluded this line from their modelling). Directly related, Katz et al. (2023) show that in cloudy one can add an extra, hard ionization source, which then, for fixed ionization parameter and gas-phase metallicity, increases [OIII] Curti et al. (2023), this might be because of possible background subtraction issues. Finally, we mask both H and [NeIII]3870 because they are blended and the current decomposition method seems to be unable to provide reliable fluxes for each of the lines (see also Section 2.4).

Modelling of the spectro-photometric data
We use the SED fitting code Prospector (Johnson et al. 2021) to simultaneously fit the photometric and spectroscopic data. We adopt a similar, 13-parameter physical model and priors as described in Tacchella et al. (2022c). Specifically, we fit for the stellar mass, stellar and gas-phase metallicities, dust attenuation (using a twocomponent dust model, which includes a diffuse component for the entire galaxy, extra attenuation around the youngest stars (< 10 Myr), and a flexible slope; 3 free parameter) and ionization parameter for the nebular emission. We adopt the flexible SFH prescription  with 6 time bins with the bursty-continuity prior (Tacchella et al. 2022c). This model includes 5 free parameters which control the ratio of SFR in six adjacent time bins; the first two bins are spaced at 0 − 5 Myr and 5 − 10 Myr of lookback time, and the remaining four bins are log-spaced to = 20. This implies that we are not able to infer ages below 2.5 Myr, which we believe is a realistic lower limit for galaxies. We also fit for a nuisance parameter that scales all the emission line fluxes by a constant factor to account for flux calibration offsets or, to first order, slit losses or physical effects (see below for details). For all the fits, we assume the MESA Isochrones and Stellar Tracks (MIST) stellar models (Choi et al. 2017) and a Chabrier (2003)

initial mass function (IMF).
Prospector has the ability to model (and fit) nebular emission line fluxes using the nebular line and continuum emission predictions provided by the Flexible Stellar Population Synthesis (FSPS) code, as described in Byler et al. (2017). These predictions are based on interpolated cloudy (version 13.03; Ferland et al. 2013) photo-ionization models that were constructed on a grid of ionization parameter ( ), gas-phase metallicity, and Single Stellar Population (SSP) age, using stellar spectra of the same age and metallicity as the ionizing source. The solar nebular abundances are taken from Anders & Grevesse (1989). The cloudy modelling used a constant electron density and the highest ionization parameter in the grid is log( ) = −1, making it difficult to model extreme [OIII] ratios (e.g. Ferland et al. 2013;Katz et al. 2023). Critically, the nebular emission is scaled by the number of ionizing photons predicted at each SSP age for the given SFH.
Crucial assumptions in this model for the emission lines are that all emission lines are powered by star formation, the LyC escape fraction is zero, and there is no dust absorption of LyC. Under these assumptions the predicted nebular emission line fluxes are consistent with the ionizing flux of the modelled stellar population. As discussed in Tacchella et al. (2022a, see also Smith et al. 2022, these assumptions do not necessarily hold at these early epochs: LyC absorption by dust can be important (absorbing up to 50% of the LyC photons during a starburst phase), about 5 − 10% of the LyC photons are absorbed by helium, and a small fraction of the LyC can also escape the galaxy. However, we add additional flexibility in two ways. First, we allow the gas-phase and stellar metallicity to take on different values, which in detail means that the shape of the model ionizing spectrum may be different than that used to predict the nebular lines ratios, though the normalisation (and hence Balmer line flux) is kept consistent. Second, we rescale all the model emission-line fluxes multiplying them by f scale , a fittable nuisance parameter (we do not rescale the nebular continuum). This constant factor can account for flux cali-    bration offsets with NIRSpec or, to first order, slit losses, differential magnification effects 3 or physical effects like LyC escape or dust absorption.

Different runs
In addition to the fiducial Prospector setup described above, we run three further configurations in order to investigate how our results depend on the modelling assumptions as well as the data present. Specifically, our fiducial approach simultaneously fits the NIRCam photometry and NIRSpec emission lines with the bursty-continuity prior ("EL & bursty" fits). Since the inferred stellar population parameters can significantly depend on the assumed SFH prior Leja et al. 2019;Suess et al. 2022;Tacchella et al. 2022b,c;Whitler et al. 2023), we vary the SFH prior. Specifically, the bursty-continuity SFH prior assumes that Δ log(SFR) between adjacent time bins is described by Student's t-distribution with = 1.0 and = 2. We also run fits with the standard continuity prior ("EL & continuous" fits), which assumes = 0.3, which weights against sharp transitions between the different time bins of the SFH ). In addition, to study the effect of including emission lines in the fitting, we fit those two SFH prescriptions to only the photometry, i.e. ignoring the emission line constraints ("phot-only & bursty" and "phot-only & continuous").
The key results of our fiducial fits are tabulated in Table 1, while  Table 2 shows the effects on the inferred stellar masses, SFRs and stellar ages when varying the SFH prior and including/excluding emission lines. We present and discuss these results in Section 3. We focus on the impact of including emission lines in the fitting in Sections 3.3 and discuss the SFH prior in Section 3.4. Figure 4 shows the goodness of our fiducial fits (i.e. fitting both photometry and emission lines with the bursty SFH prior): the top panels compare the inferred model and observed SEDs, while the bottom panels compare the individual emission lines. In the top panels, the observed NIRCam photometry obtained with forcepho is indicated with blue circles, while the inferred model photometry is shown as orange squares. The absolute values of for the individual bands are typically < 1 and the total 2 tot values are 1.7, 1.6 and 0.8 for Galaxy ID 04590, 06355 and 10612, respectively.
The emission lines are overall well reproduced (bottom panels of Figure 4)  Importantly, the bottom panels of Figure 4 include the nuisance parameter that rescales all the emission fluxes. We find that the emission lines predicted from the model need to be suppressed by 40 − 70% in order to be consistent with the observed fluxes. Not including this nuisance parameter would lead to nonphysical high stellar metallicities (i.e. > Z ) in order to suppress the ionizing photon production. Specifically, the rescale factors scale are 0.57 +0.07 −0.07 , 0.50 +0.11 −0.07 , and 0.31 +0.05 −0.04 for Galaxy ID 04590, 06355 and 10612, respectively. Our inferred rescale factors are physically plausible. They are larger than the slit loss corrections (12%, 29% and 39% for Galaxy ID 4590, 6355, and 10612), implying that the slit loss can only partially explain our inferred scale , and that the remaining 31%, 21% and 30% need be explained via other means, for example Lyman continuum absorption by dust or escape. Overall, this highlights that understandsource, emission lines from more compact regions could be more magnified, leading to a bias in the line ratios ("magnification bias"). This effect is probably negligible for Galaxy ID 06355 and ID 10612, which both have a very low magnification (Table 1).
ing slit loss corrections is fundamentally important to shed light onto both the stellar populations but also radiative transfer effects.
Finally, we compare the goodness of the resulting fits with the four different setups introduced above. The current uncertainties and this small sample do not allow us to infer whether any of the models is preferred (i.e. Bayes factor is inconclusive). However, we find that the bursty SFH prior consistently produces comparable or lower 2 values both when including or excluding the emission lines in the fitting than the standard continuity prior: we find -for Galaxy ID 04590, 06355 and 10612 -2 tot = 1.2, 3.0 and 1.2 and 2 tot = 1.3, 3.0 and 1.6 for the standard continuity prior with and without emission line constraints. Including or excluding emission lines has only little effect on the resulting 2 values.

Magnification factors
We use the magnification factors derived by Curti et al. (2023) for each galaxy (Table 1). Briefly, Curti et al. (2023) adopted the publicly available lens model of Mahler et al. (2023) to derive the magnification corrections, which combine ancillary HST with novel JWST/NIRCam data to better constrain the cluster mass distribution. For each object, Curti et al. (2023) derived the magnification maps for each target redshift and then obtained the fiducial magnification factors by averaging its value within a 1 -wide box around the central coordinates of each galaxy.
Other lens models have been published in the literature (e.g., Pascale et al. 2022; Caminha et al. 2022). The two weakly magnified objects (Galaxy ID 06355 and 10612) are robust, but we note that for the highest magnification galaxy ID4590, the different values obtained by the available models can impact the stellar mass and SFR determination for this source by up to 30 − 50%.

RESULTS
After showing in the previous section that our fits are able to reproduce the observational data, we present and discuss here the key results for our fiducial Prospector run regarding the early stellar mass build-up (Section 3.1) and enrichment (Section 3.2). We then discuss why self-consistent emission line modelling is important (Section 3.3) and what the effects are of changing the SFH prior (Section 3.4). Finally, we close by highlighting several caveats and future improvements (Section 3.5). Figure 5 plots the inferred SFHs and also gives the stellar masses and mass-weighted ages. The SFHs are increasing for all three galaxies and have similar stellar masses of about 10 8 , with Galaxy ID 04590 and 10612 being the least massive ( ★ ≈ 10 7.8 ) and and Galaxy ID 06355 being the most massive galaxy ( ★ ≈ 10 8.6

Stellar mass build-up
). Galaxy ID 04590 has increased its SFR from < 0.1 M yr −1 100 Myr ago to > 5 M yr −1 in the past 5 Myr, i.e. this galaxy is undergoing a recent burst. This is consistent with its young, massweighted age of 50 = 4 +177 −3 Myr. We note that older ages are still possible and difficult to rule out because of the strong outshining effect due to the ongoing starburst. This young age is consistent with the morphological appearance, which indicates that this source is compact with a size of less than 0.2 kpc.
Galaxy ID 06355 is going through an even larger burst of star   (Asplund et al. 2009). We find a good agreement between the gas-phase metallicity and the directly inferred oxygen abundance. We find a Z gas Z ★ for Galaxy ID 04590, while Galaxy ID 06355 and 10612 have Z gas Z ★ -but we note that the stellar metallicity is uncertain.
formation. In this case, older ages can largely be ruled out 4 and we inferred the youngest age of this sample with 3 +186 −1 Myr. It is interesting to note that this galaxies is more extended than Galaxy ID 4590 and is a clumpy object (Figure 2), consistent with the idea that mergers or an internal, gravitational instability triggered this intense starburst. Galaxy ID 10612 has an age of 50 = 3 +59 −1 Myr. The SFH shows a hint of elevated star formation 30 − 60 Myr ago. Together with insight that this source is made-up of two components (Figure 3), it is consistent with a picture where we are witnessing a merger in progress.
For comparison, Carnall et al. (2023) fitted the HST and JWST photometry with the SED fitting code Bagpipes (Carnall et al. 2018), inferring stellar masses (corrected for the different assumptions related to magnification) of log( ★ / ) = 7.5 +0.  Carnall et al. (2023) are at lower end of our quantiles, but still consistent with our estimates. However, Carnall et al. (2023) seem to rule out any underlying older component for those galaxies, which contrasts with our estimates. It is not straightforward to pinpoint the cause for this difference, but it most probably has to do with the different SFH prescriptions. Carnall et al. (2023) assume a simple constant SFH, while we assume a more flexible SFH. In addition, different extraction methods of the photometry could also lead to some differences. Finally, Carnall et al. (2023) only fit the photometry, while we fit both photometry and emission lines. As shown in Section 3.3, we would however expect that this leads to the opposite trend: excluding emission lines leads typically to older ages (Table 2). Figure 6 investigates our inferred gas-phase and stellar metallicities. In particular, the top panel compares our gas-phase metallicities (Z gas ) from the Prospector SED fitting with the ones from Curti et al. (2023), which used the direct e method enabled by the detection There is no clear trend between Z gas and sSFR 10 (and also stellar age). Intriguingly, Galaxy ID 04590, which has to lowest Z gas , has an extremely high Σ sSFR , i.e. is vigorously forming stars in a very compact configuration, indicating that this galaxy is likely experiencing a rapid rising phase of enrichment while being fuelled by accretion of pristine gas.

Dust & enrichment
of [OIII]4363. Encouragingly, we find very good agreement between the two approaches, confirming the rather low gas-phase metallicity of Galaxy ID 04590 with Z gas = 0.06 Z .
We find that the gas-phase and stellar metallicities are marginally consistent, though we emphasise that the inferred stellar metallicities are uncertain and that the elements that the gas-phase and stellar metallicity trace are different. Specifically, gas-phase metallicity traces the Oxygen abundance, while stellar metallicity traces mostly Iron. Therefore, one would expect an enhancement in the gas-phase metallicity of ∼ 2 − 5 times compared with the stellar metallicity (Strom et al. 2022;Arellano-Córdova et al. 2022). We find that Galaxy ID 04590 has Z gas Z ★ , while the other two have Z gas Z ★ . Figure 7 correlates the metallicities with the inferred stellar population and morphological parameters, including sSFR 10 (averaged over past 10 Myr), dust attenuation A V , SFR surface density Σ SFR , and sSFR surface density Σ sSFR . We do not find any convincing trend for Z gas with either sSFR 10 (as well as 50 ). Interestingly, the two galaxies (Galaxy ID 04590 and ID 06355) with an elevated V and Σ SFR have very different gas-phase metallicities. The possible cause from this can be found in the right panel of Figure 7: these two galaxies have very different sSFR surface densities with Σ sSFR = 328 +156 −94 Gyr −1 kpc −2 (Galaxy ID 04590) and 79 +42 −21 Gyr −1 kpc −2 (Galaxy ID 06355), given by the fact that Galaxy ID 04590 is a factor of 2 − 3 times more compact and lower in stellar mass than Galaxy ID 06355 (Figures 1 and 2). Galaxy ID 04590 is vigorously forming stars in a very compact configuration. This is consistent with a scenario in which Galaxy ID 04590 is undergoing rapid accretion of pristine gas, which leads to a steeply rising SFH and a high (s)SFR surface density. Despite having a low gas-phase metallicity, the high V can be explained by its compactness (i.e. high gas surface density 5 ). Consistent with this rapid inflow of pristine gas is the tentative evidence of Z gas Z ★ . Furthermore, the source just left of Galaxy ID 04590 in Fig. 1 might be part of this accretion stream. The two other galaxies (Galaxy ID 06355 and 10612) with higher gas-phase metallicity show more complex multi-component morphologies on kpc scale, indicating that their recent increase in SFR is driven by mergers or internal, gravitational instabilities. This 5 The dust-star geometry plays an important role in setting the attenuation (e.g., Zuckerman et al. 2021).
is consistent with the FMR analysis of Curti et al. (2023), where both Galaxy ID 06355 and 10612 are roughly consistent with the FMR at lower redshifts, while Galaxy ID 04590 deviates significantly, implying it being far from the smooth equilibrium between gas flows, star formation and metal enrichment in place at later epochs.

The constraining power of emission lines
From a pure Bayesian modelling approach, we want to include information to constrain the posterior distribution as tightly as possible. We therefore want to fit both the broad-band photometry and the emission lines. The emission lines themselves include crucial information about the gas properties (e.g., metallicity, temperature, and density) and ionizing sources (stars and other sources), including the most recent SFR. Figure 8 shows the posterior distribution of the fits including both photometry and emission lines (red; fiducial) and only photometry (blue) for Galaxy ID 06355. The results are summarised for all objects in Table 2. An obvious, straightforward difference is that the gas-phase metallicity is basically unconstrained when fitting only the photometry. In addition, since the relative line ratios are sensitive to the dust attenuation (in particular Balmer line series), the dust attenuation is more tightly constrained when including emission lines, which also affects the stellar age constraint via the rest-UV (although some degeneracy remains due to the flexible dust attenuation law). Directly related, another large difference can be seen for the inferred SFHs: without emission line constraints, the inferred stellar ages are older since the most recent SFH is less well constrained. This then also leads to larger stellar masses (up to ∼ 0.4 dex in the case of Galaxy ID 04590), indicating that most inferred properties of the galaxies are affected by whether emission lines are included or not in the fits.

Effects of the SFH prior
Since all these galaxies seem to be undergoing a recent burst in star formation (Figure 5), young stellar populations dominate the SED and outshine older stellar populations. Therefore, a significant amount of stellar mass could be hidden and is difficult to rule out, which implies that the prior on the SFH is expected to play an important role. As mentioned in Section 2.4, we run Prospector with two SFH priors: our fiducial "bursty" prior, which allows for large variations between adjacent time bins, and the standard "continuity" prior, which down-weights extreme bursts. Figure 9 shows the effect of the SFH prior on key inferred quantities, including the SFH, stellar mass, sSFR, dust attenuation, age, stellar metallicity and gas-phase metallicity. The results from the bursty and standard continuity prior are shown in red and blue, respectively. The same observational data, i.e. both photometry and spectroscopy, are fitted ("EL & bursty" and "EL & continuous" runs in Table 2).
As expected, the SFH prior affects the inferred SFH: the continuous SFH prior leads to more stellar mass being formed early and a less steep rise in recent time. This leads to overall older ages, larger stellar masses and lower sSFRs than what is inferred with the bursty prior. Specifically, as tabulated in Table 2, the stellar ages increase from 10 Myr to 100 Myr for Galaxy ID 4590 and 10612, highlighting that it is difficult to rule out early star formation. Importantly, the 50 posterior distributions of the two runs are overlapping, i.e. they are consistent with each other. There is also a big impact on the inferred stellar masses: they increase by up to 0.6 dex for those two galaxies. bursty prior continuity prior Figure 9. Effect of changing the SFH prior for Galaxy ID 06355. The results adopting the bursty prior and the continuity prior for the SFH are shown in red and blue, respectively (the red histogram is the same as in Figure 8). The continuity prior results in a more constant SFH, which leads to a higher stellar mass, lower sSFR and older age.
The stellar age constraints come from the rest-UV and 4000Åbreak (e.g., Conroy 2013). Since we are using a flexible dust attenuation law (Section 2.3), the rest-UV has only little constraining power. The 4000Å-break for our three galaxies at ∼ 7 − 9 is probed by the broad-band photometry, but -as shown in Figure 4 -there is a degeneracy with the strength of the emission lines. Specifically, both a strong continuum break and strong emission lines can lead to a boost of the longer-wavelength filter relative to the shorter-wavelength one. But why are we not able to do better when including emission lines in the fitting? The problem is that our nuisance parameter, which rescales the emission line fluxes, is degenerate with the strength of the 4000Å-break, i.e. our emission line measurements mainly constrain the gas-phase metallicity and the dust attenuation, which in second order leads to a tighter constraint on the stellar ages via the rest-UV.
Further improvements (and solving the degeneracy with slit losses) can be achieved by exactly assessing the slit position (and performing spatially resolved SED modelling) and by directly probing emission lines with the photometry, which can be achieved with narrow and medium bands (e.g., Roberts-Borsani et al. 2021;Tacchella et al. 2022c). The medium bands on/off the lines would allow the measurement of the 4000Å-break more directly.

Caveats and outlook
We list here the most important caveats in our analysis, which we hope to improve upon in the future. From the Bayesian SED modelling perspective, our goal is to include as much information as possible that can be described by our physical model. This means we want to fit both the NIRCam photometry and the NIRSpec emission line constraints. One fundamental challenge with this is that the photometry and spectroscopy could probe different regions of the galaxy. Specifically, the photometry attempts to include the total emission of the galaxy, while the spectroscopy only captures the light that is within the shutter of the NIRSpec MSA. Slit losses affect most spectroscopic observations, but because not every object can be aligned perfectly with the NIRSpec MSA, slit losses may be a larger issue for these observations. Besides this technical challenge 6 , other observational and physical effects can also lead to a mismatch between the photometry and spectroscopy. On the observational side, differential magnification across the source can lead to boosting of flux of certain galactic sub-regions (typically of more compact regions), leading to a bias between the photometry and spectroscopy (as well as of certain emission lines). On the physical side, we assume in this work that the emission lines (and the nebular continuum) are powered by ionizing photons of the stellar populations, without taking into account LyC absorption by dust of LyC escape, which could be important, in particular for objects in the starburst phase (e.g., Tacchella et al. 2022a).
In this work, we addressed these technical, observational and astrophysical issues by fitting for a nuisance parameter, scale , which scales all the emission lines by a constant factor, mainly intended to account for slit loss. We find scale to be 0.57 +0.07 −0.07 , 0.50 +0.11 −0.07 , and 0.31 +0.05 −0.04 for Galaxy ID 04590, 06355 and 10612, respectively. If this rescaling is ignored, we find extremely high values for the stellar metallicities (Z ★ > Z ). This highlights that this factor is important when performing the fitting. Our inferred scale are larger than calculated from slit loss, implying that 31%, 21% and 30% need be explained via other means, for example Lyman continuum absorption by dust or escape. Furthermore, working with unlensed galaxies will resolve the issue with the differential magnification bias.
Related to this is the assumption in our approach that stellar populations are driving the emission lines. Firstly, active galactic nuclei (AGN) and other non-stellar sources could power the emission lines as well as contribute to the rest-UV emission. Within Prospector, only the rest-frame mid-infrared emission of AGN can be modelled (i.e. dusty torus emission). There are on-going efforts to expand upon this. Secondly, a recent analysis by Katz et al. (2023) shows that the inclusion of high-mass X-ray binaries or a high cosmic ray background in addition to a young, low-metallicity stellar population can provide the additional heating necessary to explain the observed high [OIII]4364/[OIII]5007 ratio while remaining consistent with other observed line ratios. However, because it is challenging photoionization models to predict [OIII]4364 (e.g., Dors et al. 2011), we mask this line. More generally, it is notoriously difficult to distinguish between AGN and star formation as an ionization source at low metallicity. Interestingly, Galaxy ID 6355 appears to have an AGN contribution (narrow-line AGN given [NeIV]2422,2424 is detected; Brinchmann 2022). We believe that our results are robust even for this object given that such high ionization lines can be obtained by adding an Xray source without strongly affecting the other lines (Katz et al. 2023).
For completeness, we also mention the caveats related to stars. In particular, in this work we assume a solar abundance pattern scaled by the metallicity. However, there is an indication that at higher redshifts (for ≈ 2 results see Strom et al. 2022) galaxies are more -enhanced (but see Arellano-Córdova et al. 2022 for a recent analysis of this for those three galaxies, finding no evolution in this -element ratio). This is also theoretically expected given that some of the highgalaxies are old enough to have seen enrichment from intermediatemass stars, but are still young enough that Type Ia supernovae have not had time to contribute significantly to their enrichment (e.g., Kriek et al. 2016). We further assume the MIST stellar models (Choi et al. 2017), which include rotation, but not binarity. Investigating the effects the binary-based stellar model (e.g., Eldridge et al. 2017) or varying the IMF is beyond the scope of this work.
Finally, this work investigates three galaxies, the only three > 7 galaxies that have both NIRCam and NIRSpec observations. Obviously, we cannot draw strong conclusions regarding the whole galaxy population from three objects. Although we think that galaxies have mostly increasing SFHs at these early times (e.g., Tacchella et al. 2018;Endsley et al. 2021), the substantial bursts found in these galaxies might be a selection effect and not representative of all galaxies at ∼ 7 − 9 (e.g., Looser et al. 2023). Therefore, we stress that the main purpose of this work is to deliver interesting insights into the properties of those three galaxies by combining NIRSpec and NIRCam data, discuss technical but important details for doing this, and discuss how priors play an important role in driving some of the results.

CONCLUSIONS
We present a careful ISM and stellar population analysis of the ERO NIRCam and NIRSpec data of three = 7.6 − 8.5 galaxies in the SMACS cluster field. These three galaxies have diverse morphologies (Figures 1-3), from being compact (Galaxy ID 04590) to being made up of several components (at least two components in Galaxy ID 10612 and four components in Galaxy ID 06355). We perform the photometry with the new photometry tool forcepho, which is a Bayesian model-fitting program that simultaneously fits multiple PSF-convolved Sérsic profiles to all filters.
Within the Prospector framework, we fit a 13-parameter model to NIRCam photometry and the NIRSpec emission lines. This physical SED model includes a flexible SFH, a multi-component dust model including a variable dust attenuation law, different gas-phase and stellar metallicities, a free ionization parameter for the nebular emission and a nuisance parameter that scales the emission line fluxes. The latter parameter is important to account for possible slitloss effects and physical effects such as LyC dust absorption and escape. We find that this factor is important: model-based emission line fluxes need to be rescaled by factors of 0.3 − 0.6, caused both by slit loss and astrophysical reasons (for example escape of ionising radiation). Overall, we find that we are able to reproduce the pho-tometry and emission line measurements (Figure 4). An exception is the emission line [OIII]4364, which is brighter than predicted by our modelling, which could be attributed to non-stellar sources powering this line, to a conservative cloudy grid (too low ionization parameter values).
We infer for all three galaxies rising SFHs and stellar masses of ★ ≈ 10 8 ( Figure 5 and Table 1). These galaxies are all young with mass-weighted ages of 50 = 3 − 4 Myr. However, we find indications for underlying older stellar populations, implying the SFHs extend at least over several tens of Myr. We emphasise that the SFHs, stellar masses and stellar ages depend on the adopted SFH prior (Figure 9 and Table 2): assuming a bursty SFH prior leads to younger, lower-mass galaxies, which is consistent with previous studies (Tacchella et al. 2022c;Whitler et al. 2023). Emission lines are helpful to constrain mainly the gas-phase metallicity and ionization parameter, and only have a second order effect on the inferred SFHs and ages ( Figure 8 and Table 2). This is because of the aforementioned rescaling of the emission lines, i.e. we are fitting for relative emission line strength and not their absolute strength. However, we still find that emission lines help with constraining the SFHs because they pin down the dust attenuation parameters, which then in turn constrains the rest-UV emission -another important age indicator.
Focusing on metallicity, our SED fitting approach (with masking the [OIII]4364 line) delivers gas-phase metallicities consistent with the direct method from the [OIII]4364 line from Curti et al. (2023, Figure 6). We find no convincing trend between the gas-phase metallicity and the total sSFR (and stellar age) of the galaxies (Figure 7). However, Galaxy ID 04590 with the lowest gas-phase metallicity shows the highest star-formation intensity as measured by the sSFR surface density Σ sSFR (Σ sSFR ≈ 330 Gyr −1 kpc −2 ), i.e. this galaxy has very high SFR surface (Σ SFR ≈ 22 M yr −1 kpc −2 ) density for its stellar mass ( ★ ≈ 10 8 M ), indicating that this galaxy is vigorously forming stars in a very compact configuration. This is consistent with a scenario in which Galaxy ID 04590 is undergoing rapid accretion of pristine gas, which leads to a steeply rising SFH and a high (s)SFR surface density. Despite having a low gas-phase metallicity, the high V can be explained by its compactness (i.e. high gas surface density). Consistent with this rapid inflow of pristine gas is the tentative evidence of Z gas Z ★ . The two other galaxies (Galaxy ID 06355 and 10612) with higher gas-phase metallicity show more complex multi-component morphologies on kpc scales, indicating that their recent increase in SFR is driven by mergers or internal, gravitational instabilities.
In summary, our work highlights the great potential for combining photometric with spectroscopic JWST data to study early galaxy formation, such as the JADES survey (Rieke et al. 2019;Bunker et al. 2023;Cameron et al. 2023;Curtis-Lake et al. 2023;Robertson et al. 2023;Saxena et al. 2023;Tacchella et al. 2023).