Spatially Resolved Properties of High Redshift Galaxies in the SMACS0723 JWST ERO Field

Galaxies


INTRODUCTION
Corresponding author: Clara Giménez-Arteaga clara.arteaga@nbi.ku.dk * Hubble Fellow By characterizing the physical properties of galaxies in the redshift range 5 < z < 10, we can study the epoch of reionization, when the Universe experienced its last phase transition (see e.g., Treu et al. 2013;Mason et al. 2018;Robertson 2022, for a review).With well-sampled photometry of high redshift galaxies, we can robustly model their spectral energy distributions (SEDs) and infer the properties of their stellar populations.Up until now, the rest-frame optical emission from galaxies was unavailable at z > 7, having been redshifted to the part of the near-infrared spectrum where our facilities lacked sensitivity and spatial resolution.While the rest-frame UV emission we have had access to is a good tracer of unattenuated star formation, it is a poor tracer of stars older and less massive than O and B-type, that make up the bulk of total stellar mass for populations older than a few Myr.
The latest addition to the space fleet of telescopes, the James Webb Space Telescope (JWST ), has unprecedented sensitivity and spatial resolution in the near infrared.This has opened up a new window into the restframe optical emission of high redshift galaxies, allowing us to understand their stellar populations for the first time.The Near-Infrared Camera (NIRCam; Rieke et al. 2005) on board JWST allows us to reach this spectral range with a unique depth and resolution.The Near Infrared Spectrograph (NIRSpec; Jakobsen et al. 2022) provides high resolution spectroscopy in the nearinfrared, which is key to robustly determine the redshift.This improves the modelling of the SEDs by constraining a free parameter, thus breaking the degeneracies that the redshift has with age and dust (see e.g., Conroy 2013, for a review).
With JWST we can extend for the first time resolved studies beyond redshift ∼ 2, introducing the possibility to study in unique detail the first generations of galaxies (e.g., Chen et al. 2022;Hsiao et al. 2022).These resolved studies will allow us to place unique, new constraints on the formation and evolution of the first galaxies: their mass assembly histories, modes of growth, chemical enrichment, and earliest quenching mechanisms.In order to build a complete picture of galaxy assembly, a resolved view of its components is required, to fully understand the interplay between the stellar population, dust and gas in z > 6 galaxies.
The impact of having resolved observations has so far only been studied at low redshifts (z 3).Various works have compared the inferred physical properties obtained with resolved and unresolved observations (see e.g., Wuyts et al. 2012;Sorba & Sawicki 2015, 2018;Vale Asari et al. 2020;Fetherolf et al. 2020), with diverse conclusions.A resolved approach can have multiple advantages, such as decreasing degeneracies in the stellar population synthesis models and producing more realistic star formation histories (Pérez-González et al. 2022).In highly star forming galaxies, the outshining of old stellar populations by young ones is of particular importance, which a resolved analysis could untangle.Sorba & Sawicki (2018) find that the total stellar mass can be underestimated by a factor of ∼ 5 without taking this outshining effect into account with spatiallyresolved SED modelling.
In this paper we present the observations and first results of a multi-wavelength analysis of five high redshift galaxies (5 < z < 9) observed with JWST, to study the spatially-resolved properties of their stellar populations.These galaxies are the highest spectroscopically confirmed targets in the JWST ERO SMACS0723 field.Using NIRCam imaging in 6 bands spanning the wavelength range 0.8 − 5µm, we perform spatially-resolved SED fitting with Bagpipes (Carnall et al. 2018).There have been multiple works on these targets, albeit always from an integrated perspective (e.g., Carnall et al. 2022;Tacchella et al. 2022;Schaerer et al. 2022;Trump et al. 2022;Rhoads et al. 2022;Curti et al. 2022;Brinchmann 2022;Arellano-Córdova et al. 2022;Heintz et al. 2022;Fujimoto et al. 2022).These papers derive different stellar masses, some of them with extremely early star formation histories.Here we present the first spatially resolved analysis on these high redshift galaxies, which we propose as a more robust approach to accurately calculate their stellar masses.
This paper is structured as follows.In Section 2, we introduce the JWST observations and data reduction procedure.Section 3 describes the methodology we use for the modelling of the SEDs with Bagpipes.In Section 4, we present the main results and inferred properties of our sample, both in integrated and resolved approaches, as well as discussing the implications of our analyses.Finally, in Section 5 we present the summary and conclusions of our work.Throughout this paper, we assume a simplified ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3 and Ω Λ = 0.7.No lensing correction is applied throughout this work.Hence, intrinsic stellar masses and star formation rates can be obtained by dividing by the magnification factor (µ, see Table 1).

DATA AND OBSERVATIONS
We use the public data of the galaxy cluster SMACS J0723.3-7327(SMACS0723), observed by JWST as part of the Early Release Observations (ERO; Programme ID 2736, Pontoppidan et al. 2022).We use the Near-Infrared Camera (NIRCam; Rieke et al. 2005) photometric data from the catalog reduced by Brammer et al. (in prep).The data has been reduced using the public software package grizli (Brammer 2019;Brammer & Matharu 2021;Brammer et al. 2022).The photometry is corrected for Milky Way extinction assuming E(B − V ) = 0.1909 (Schlafly & Finkbeiner 2011) and the Fitzpatrick & Massa (2007) extinction curve.The images are PSF-matched to the F444W band on a common 0. 04/pixel scale.We adopt the PSF models for use with the grizli mosaics 1 , which are based on the WebbPSF models.We compute matching kernels for each of the PSFs to the F444W PSF using a Richardson-Lucy deconvolution algorithm (Richardson 1972;Lucy 1974), and then convolve the images with the resulting kernels to match the PSF resolution in F444W.
In this work we focus on the five galaxies at 5 < z < 9 that have spectroscopic redshifts confirmed from NIR-Spec observations, presented in Carnall et al. (2022).These are detected in the 6 deep NIRCam imaging filters F090W, F150W, F200W, F277W, F356W, and F444W.The targets also have shallower imaging data obtained with NIRISS with the F115W and F200W filters, which we exclude in the analysis presented in this work due to their lower resolution.The sample spans a redshift range from 5.275 for the closest galaxy, up to 8.498 for the highest redshift object (Carnall et al. 2022).The sources have magnification factors between 1.6 up to 10.1 in the Glafic lens models from Oguri (2010).They are not significantly distorted by the gravitational lensing so that it affects the spatially resolved SED fitting.
Table 1 provides the basic information for the targets.We do not apply any lensing correction to our results, although we report on Table 1 the lensing factors for these sources presented in Carnall et al. (2022).Figure 1 displays the cutout images of the five galaxies in all available observed bands, as well as the colour images built combining the F150W, F277W, and F444W filters.

SED fitting with Bagpipes
To model the spectral energy distribution of the individual pixels and derive the physical properties, we use  Carnall et al. (2022), where the lensing factors (µ) are taken from Oguri (2010).
We set the NIRSpec spectroscopic redshifts indicated on Table 1, in order to break the degeneracy with age and dust that a photometric or uncertain redshift would introduce.We use the SPS models by Bruzual & Charlot (2003) and include the nebular emission with Cloudy (Ferland et al. 2017), extending the Bagpipes default grid (that normally reaches up to log 10 (U ) = −2) so that the ionization parameter, U , varies from −3 < log 10 (U ) < −1, since z > 6 galaxies display higher ionization parameters than low-redshift galaxies (e.g., Sugahara et al. 2022;Curti et al. 2022).We assume a Kroupa (2001) initial mass function (IMF), and a Calzetti et al. (2000) attenuation curve, in order to reduce the number of free parameters in our fits.We choose a constant star formation history model, following Carnall et al. (2022), which seems adequate to fit our galaxies (we obtain reduced χ 2 values within the range 0.1 − 7.5, shown in further detail in the following sections).The formation of very young, low-mass, and lowmetallicity galaxies is likely bursty, and a constant SFH accurately resembles this on short timescales.We let the maximum age grid to vary from 1 Myr to 1 Gyr, to allow for the presence of more evolved stellar populations, and further limited by the age at the given redshift.
We set the visual extinction to vary from A V = 0 to A V = 2, and the metallicity from 0 to Z , with uniform priors.Even though the metallicity has been calculated with integrated NIRSpec spectra for these targets (see e.g., Schaerer et al. 2022;Curti et al. 2022;Brinchmann 2022), they obtain varying results within our allowed range.Moreover, we want to allow for spatial variation across the galaxy, thus we do not set Z as a fixed parameter in the fit.Finally, we set the lifetime of birth clouds to 10 Myr.

Pixel-based modelling
In this work we perform SED fitting on a pixel-bypixel basis.With the setup described in the previous subsection, we fit the spectral energy distribution and infer the physical parameters of each individual pixel.1.The RGB images are built combining the F150W (B), F277W (G) and F444W (R) PSF-matched images (to the F444W band).The physical scale is calculated in the lens plane.
This allows us to recover the 2D distribution of properties such as the stellar mass and star formation rate (SFR).The pixels in our maps correspond to physical sizes between 180 and 240 parsecs.In order to fit the SED of individual pixels, we impose a signal-tonoise ratio (S/N) threshold of 2 on both the F150W and F200W bands, which are the noisiest.We find that this threshold is enough to produce trustworthy fits, obtaining good reduced χ 2 values in the fits of individual pixels, as we show in more detail in the following section.To produce the maps of the physical properties and study their spatial distribution, we display the 50th percentile of the inferred parameter, calculated with the posterior distribution that Bagpipes provides.We can also present the uncertainties for each pixel extracted from the 16th and 84th percentiles of the posterior distribution (see Carnall et al. 2018 for details).

RESULTS AND DISCUSSION
In this section we present and discuss the results of our study.We provide both an integrated and a spatially resolved analysis of the physical properties that we infer for the five targets that comprise this work.

Spatially Resolved Physical Properties
For all galaxies studied, we have produced maps of the physical parameters inferred with Bagpipes.These include the star-formation rate surface density (SFRD), the stellar mass surface density (SMD), the visual extinction A V , the mass weighted age, the UV slope β, and the equivalent width (EW) of the [O III] and Hβ emission lines.The last two are measured from the Bagpipes posterior SEDs.We also display maps of the empirical colours F150W-F277W, which is a proxy for the UV slope, and F356W-F444W, which generally traces the strength of the inferred [O III]+Hβ emission (except for the lowest and highest redshift galaxies, where no available couple of bands capture the lines and continuum accordingly).
Figures 2-6 display the resulting maps for the five galaxies presented in this work, from the lowest redshift z = 5.275, to the highest z = 8.498.Firstly, we see that all galaxies are resolved and display strong empirical colour gradients, both in the blue bands as well as the red bands.These gradients appear on larger scales than the FWHM of the F444W PSF (0. 145, equivalent to ∼ 3.6 pixels), confirming that we can resolve trends and structures in our sample.In general, we find that even at these early times, most of the galaxies display multiple star forming clumps, traced by the F200W contours, as is found also by other recent works (e.g., Claeyssens et al. 2022;Treu et al. 2022;Chen et al. 2022).These are regions of very high inferred equivalent widths of the [O III]+Hβ emission (in the range ∼ 300 − 4000 Å restframe), embedded within larger structures that are not undergoing a burst of star formation.In these regions with extreme line EWs, the inferred ages are extremely young (<10 Myr), corresponding to a bursty clump of young stars that is resolved in targets ID10612 (Figure 4) and ID6355 (Figure 5), and marginally unresolved in ID5144 (Figure 3) and ID4590 (Figure 6), given the scale of the F444W PSF FWHM.Around these high-EW bursty clumps, we also find underlying older stellar populations (∼100 Myr), which would be missed in an integrated analysis, as we will discuss in more detail in the next subsection.
Figure 2 displays the maps for galaxy ID8140, at redshift z = 5.275.We obtain smooth maps for all physical properties and colours.The SFRD and SMD appear entirely co-spatial, and the UV slope traces perfectly the dust obscuration map.This galaxy has ∼ 2 times older stellar populations compared with the average of the rest of the sample, in line with being the galaxy at lowest redshift.It shows strong colour, A V and UV slope gradients, with two distinct clumps, one red and one blue.This could indicate that this galaxy is undergoing a merger, even though the EWs (the lowest of all targets) and ages show little variation across the object.
Figure 5 shows the resulting maps for the galaxy ID6355 at z = 7.665.It is the galaxy with most extreme EWs (reaching ∼ 4000 Å).We see that the region with very high EW is very extended for this galaxy, leaving barely a shell where we find underlying older stellar pop-ulations, which are otherwise outshined (or not present) by the younger stars in these strong line emission regions.The clear gradients in the empirical colours give us confidence that this shell is real and not an artifact of the age-dust degeneracy in the SED fitting process, which we also discuss in further detail in §4.4.On top of this, the shell is larger than the PSF scale.In Figure 7 we present and analyse the fits for three individual pixels within this source, so that we can study further whether the "shell" of older stars in this particular galaxy is real or an artifact.We select the pixels A, B, and C that are indicated in the age and EW maps in Figure 5, since they appear to be very distinct regions within this galaxy.Albeit being towards the edge of the galaxy, all three pixels fulfill our S/N threshold, so that we can produce robust fits (with reduced χ 2 values within 0.10 − 0.53).These pixels are also far enough from each other so that the PSF is not blending the information they encode, and we can thus resolve their different stellar populations.
We can clearly see that the SEDs look different, reflecting the gradients that we already see in Figure 5, both on the maps of the inferred physical parameters, as well as the empirical colour maps.The greatest difference is observed in the strength of the inferred [O III]+Hβ emission lines, since the SED for pixel C has extreme EW, reaching 4264 ± 533 Å rest-frame.This yields a considerable difference in the inferred ages, with pixel A having a mass weighted age of 159 +115 −108 Myr, and pixels B and C displaying very young stellar ages under 10 Myr.
Figure 6 shows the results for the galaxy ID4590, which is the highest redshift z = 8.498 in this work.It is a very compact source, with only 31 pixels where S/N>2 for both F150W and F200W.We see a star forming clump, which also corresponds to the highest inferred stellar mass, and towards the dustiest zone in the A V map.We see a marginally unresolved centrally located clump of young stellar population, which is not entirely co-spatial with the star-forming burst, but traces perfectly the higher equivalent width of the inferred [O III] and Hβ emission lines.On top of this, the empirical colour maps display a clear gradient, which follows the ones observed for the SFRD, SMD, A V and β maps.The rest of targets exhibit similar trends for all physical properties, with the main characteristic being this region with extremely high EWs and therefore very young stellar populations.

Integrated Analysis
Besides providing an invaluable insight into the internal structure of galaxies, we want to test whether spatially resolved observations yield other consequences, such as inferring different physical properties, compared to only integrated measurements.
To perform this test, we sum the photometry in each observed band for the pixels that fulfil our S/N criteria, so that we only consider the same pixels that we fit in the spatially resolved analysis shown in Figures 2-6.With the sum of the photometry in each filter, we then use Bagpipes to find the best fit SED and infer the integrated physical parameters.We use the exact same set up as for the spatially resolved run, described in §3.1.We present the integrated fits with the best fit SEDs in Figure 8.The inferred integrated physical properties for each galaxy can be found in Table 3, in Appendix B. In Figure 8, we can see that both the resolved and integrated models (red and black curves, respectively) fit adequately the photometry (turquoise points), with reduced χ 2 values within the range 0.1 − 7.5.The surprising finding is that both best fit SEDs are considerably different.For all galaxies, we find that the high equivalent widths that we could spatially locate in the resolved analysis within a clump, now completely dominate the overall fit.This results in inferring extremely young ages in the integrated light, and potentially too low stellar masses as a result.
Figure 9 shows the comparison between the stellar mass estimates that we infer in the spatially resolved analysis and the integrated fit.Table 2 provides these values, as well as the mass weighted ages.We obtain the  2. The one-to-one line is indicated, as well as the 0.5 and 1 dex offset lines.

. 5 d e x 1 d e x
spatially resolved mass estimate by summing the stellar mass inferred in each individual pixel.The resolved average age is inferred as the mean mass weighted age of all pixels.In both cases, the resolved uncertainties reported in Table 2 are calculated with the 16th and 84th percentile of the summed posterior distribution over all pixels.For the integrated run, the output from Bagpipes is directly reported, and the uncertainties are calculated with the 16th and 84th percentiles from the posterior distribution.We find that the integrated run estimates systematically lower stellar masses than the spatially resolved one, from ∼0.5 up to 1 dex lower, seemingly without any trend with the redshift.This, as explained above, is a consequence of being forced to choose extremely young stellar populations to fit the integrated light, due to the strong emission lines dominating it completely.We see this in Table 2, since all galaxies have an integrated age of under 10 Myr, except the one at lowest redshift, with a best fit age of 14.9 +1.7 −1.0 Myr.We demonstrate this by fixing the age in Bagpipes to the average stellar age inferred in the spatially resolved analysis.By doing this, the estimate of the stellar mass in the integrated run increases, retrieving closer values to the spatially resolved masses.Moreover, we see that the average mass weighted ages in the resolved analysis are all > 10 Myr, and significantly older than the integrated ages, even when considering the large uncertainties associated with the age estimate.
Figure 10 shows how the SFH affects the inferred stellar mass.We plot the sum of the SFH inferred for the spatial pixels, as well as the SFH estimated in the unresolved analysis, for the galaxy ID10612 at z = 7.663.The integrated SFH consists of a single burst with very young age (∼ 2 Myr), whereas the spatially resolved SFH is a distribution that covers a wider age range, reaching up to ∼ 300 Myr.For this galaxy, this would mean a formation redshift of z ∼ 12.We see that, whereas the integrated analysis forms all stellar mass within less than 10 Myr, this only corresponds to ∼ 6% of the spatially resolved stellar mass, directly proving where the mass discrepancy is coming from.We obtain the same results in the SFH comparison for all galaxies studied in this work.This could considerably change our current picture of mass assembly in the early Universe, particularly while our samples are limited to the brightest candidates in small numbers.These systematics would affect from the stellar mass functions that we have derived so far at high redshifts, to our cosmological models of galaxy formation and mass build-up, since all our observations and mass estimates at high redshift until now have been based on integrated measurements.Overall, the nature of these galaxies can completely change by having a resolved picture.
On top of this, we adopt a parametric constant SFH model.Some works have shown that using instead a non-parametric SFH model can lead to inferring larger stellar masses by up to 1 dex (all in an integrated approach), in particular for galaxies like the ones studied here, where a burst with very young stellar ages dominates the integrated light (e.g., Leja et al. 2019;Lower et al. 2020;Tacchella et al. 2022;Whitler et al. 2022).Non-parametric SFHs therefore seem like the better option when considering integrated SEDs, since we have shown that there can be significant spatial variation in star formation (see Figure 10).Moreover, works like Bisigello et al. (2019) argue that we need to include the two MIRI bands in order to constrain better the stellar mass estimates in high redshift galaxies, particularly in young galaxies with nebular emission lines, such as the ones we study here.

Comparison with Other Works
To put our results into context, we compare the physical parameters that we infer with other published works on these targets.As stated before, only integrated analyses are available in the literature.We expect our integrated measurements and estimates of the physical properties to be lower than any other work, since aperture photometry yields larger fluxes on all bands (∼ 20% larger for the galaxy ID10612 at z = 7.663 if we use a 0. 5 aperture), given that we only consider pixels that fulfil our S/N criteria.On top of this, previous works have used varying data reduction, zero points, SED fitting codes, SFH models, attenuation curves, magnification corrections, amongst other differences.Therefore, here we mostly focus on comparing the physical nature of the sources, as inferred by different works.
The mass-weighted ages that we infer in our integrated analysis are all consistent within the uncertainties with those found for these same five targets by Carnall et al. (2022), from which we reproduce the SED fitting assumptions and parameters, except using a differently reduced photometry, a Calzetti et al. (2000) attenuation curve, as well as extending our nebular grid as explained in §3.Tacchella et al. (2022) studies the stellar populations of three of our targets, the ones at highest redshifts z = 7.663, z = 7.665 and z = 8.498 (IDs 10612, 6355 and 4590, respectively).They use Prospector with a flexible non-parametric SFH prescription instead to model the SEDs.They find that the highest redshift galaxy (ID4590) is undergoing a recent burst, inferring a young stellar age under 10 Myr, like we obtain in our integrated analysis.They cannot rule out older stellar ages in their analysis.In our spatially resolved maps (see the age map in Fig. 6), we find that most pixels have a mass weighted age above ∼ 100 Myr, except the centrally located young burst, which dominates the integrated light.The particularly striking case is the z = 7.665 galaxy (ID6355).Tacchella et al. (2022) also infer an extremely young age of 3 +29 −1 Myr, and they rule out the presence of older stellar populations, since the extreme [O III]+Hβ lines dominate the emission.In this work we find a shell of older stars, confirmed by the empirical colour gradients as well as the other tests discussed in §4.4 and §A.Finally, for the target ID10612 at z = 7.663, they argue that its SFH together with its morphology, could indicate that this source is undergoing a merger.This is consistent with what we find, both in terms of the empirical colour gradients, as well as the clumpiness of the blue F150W and F200W bands.Moreover, we find two distinct populations within the galaxy in the spatially resolved maps (see Figure 4).Consistent with our results here, they find that inferring older stellar ages (in their case by adding emission line constraints with NIR-Spec spectra) leads to larger stellar masses of up to 1 dex for the galaxy ID4590 at z = 8.495, which is exactly what we find in our integrated versus spatially resolved comparison.
Our work confirms the issue discussed in Tacchella et al. (2022); Topping et al. (2022) and Whitler et al. (2022), where young stars outshine and dominate the emission when compared to older stars, the presence of which is difficult to rule out via integrated measurements.This leads to inferring lower stellar masses.In Tacchella et al. (2022) they conclude that the SFH prior is of vital importance, and they only infer stellar ages older than 10 Myr in their fits when using a non-parametric SFH model with a continuous prior and older populations present.This leads to an increase in the stellar masses of up to 0.6 dex.Whitler et al. (2022) also use various SFH models to explore the potential presence of old stellar populations in seemingly young galaxies.They find stellar masses larger by up to an order of magnitude with non-parametric SFH versus constant SFH models.The impact of the SFH model in the inferred stellar mass has been studied by other works (see e.g., Leja et al. 2019;Suess et al. 2022), with similar results.
Outshining and its effects on stellar mass estimates have been studied at lower redshifts (see e.g., Maraston et al. 2010;Pforr et al. 2012).Our results agree with previous works such as Sorba & Sawicki (2018), where they find a discrepancy in the inferred stellar masses of up to a factor of ∼ 5, when having resolved SED fitting, albeit their study only reaches z = 2.5.Moreover, they propose that unresolved studies should apply corrections to their mass estimates.This resolves the mass missing problem, in which a tension is found between the observed stellar mass density of the Universe and the star formation rate density (see e.g., Leja et al. 2015).By correcting stellar mass functions with resolved estimates, they find that these agree better with the observed star formation densities collected by Madau & Dickinson (2014).Leja et al. (2020) also solve this discrepancy by using flexible non-parametric SFH priors, which produce older ages, thus inferring ∼ 50% higher stellar mass density.Endsley et al. (2022) study a population of UV-faint galaxies at a similar redshift range z ∼ 6.5 − 8, also finding that the SEDs are dominated by young stellar populations, exhibiting low masses.They find the majority of their objects to appear very blue (β ∼ −2), with some dusty galaxies (β ∼ −1).With our integrated analysis, we find three galaxies with β ∼ −2, and two targets with a value closer to −1.In the spatially resolved maps, we find values −2 < β < −1, with some very dusty regions reaching values around β ∼ −0.5.
At this redshift range, the majority of targets with high EW([O III]+Hβ) have an inferred young stellar age when considering a constant SFH, just as we find in our integrated analysis (see Fig. 9 in Endsley et al. 2022).On the other hand, there are works that find evolved stellar populations (>100 Myr) at even higher redshifts, such as Furtak et al. (2022) with z ∼ 10 − 16 candidates in the SMACS0723 field, and Leethochawalit et al. (2022) with 7 < z < 9 photometrically-selected galaxies in the GLASS-JWST ERS program, which infer a median mass-weighted age of 140 Myr.
From a resolved point of view, in a sample of z ∼ 6−8 galaxies in the Extended Groth Strip (EGS) field, Chen et al. (2022) find multiple clumps dominated by young stellar populations, as well as significant variations in the equivalent width of the [O III]+Hβ lines.They find EWs with extreme values such as the ones we find for most of our targets (of the order ∼ 300 − 3000 Å), which also yield young ages in their fits, confirming once more what we are finding for these high redshift targets.Moreover, Pérez-González et al. ( 2022) also find strong [O III]+Hβ emission in a spatially resolved analysis using HST and JWST data from the CEERS survey in the EGS.They also link these findings with very young starburst with possibly an underlying older stellar population.
In summary, our results are consistent with the works that have been published so far studying these same galaxies, or targets at a similar redshift range with JWST.By integrating our maps, we can produce similar results and draw equivalent conclusions to the integrated works performed so far in these sources.By producing a spatially resolved analysis, we can demonstrate the presence of underlying older stellar populations that are otherwise outshined in the integrated analyses, inferring larger stellar masses and considerably affecting our picture of the nature of these high redshift galaxies.

Caveats
As briefly mentioned before, one could argue that the "shells" of older stellar populations where the young stars with high EWs are embedded, could be instead a result of the dust-age degeneracy present in SED fitting softwares.To test if the gradient in age is real, we perform a test in which we fix the extinction to the value given by Tacchella et al. (2022).We find a similar gradient, where there is a shell of older stellar populations surrounding the bursty young star forming region.The effects of dust and age are now blended into the age map, since we fix the A V , but the gradient persists.On top of that, the individual fits that we obtain fixing the dust obscuration are very poor, compared to leaving A V as a free parameter in the Bagpipes fit.This gives us confidence that our fits with A V as a free parameter sample better the galaxy properties, and the gradient observed is real.In Appendix A, we discuss in more details the age uncertainties, focusing on the target ID6355 at z = 7.665.
Another caveat that our spatially resolved analysis could have is whether the process of PSF-matching affects our inferred maps.One could argue that the mass weighted age map could result as an artifact of the PSFmatching procedure, where flux is re-distributed radially to match the resolution of the F444W band.We only observe this radial distribution on the age and EW maps.We observe horizontal gradients on all the rest of physical properties, as well as the empirical colours.The gradients that we observe extend across spatial scales larger than the FWHM of the F444W PSF.We therefore conclude that this is not an effect of PSF-matching, but a true young stellar population clump centrally located in most galaxies.
Finally, besides the S/N threshold that we impose in the noisiest bands, one could still doubt whether there is enough S/N per pixel to be able to infer robust physical parameters.To test this, we apply a Voronoi tesselation binning method on the targets, in order to achieve bins with a constant minimum S/N across the image and filters.Imposing a minimum S/N of 5 or even up to 10 in all bands, we find the same gradients and trends that we observe in the maps of the various inferred physical parameters in all galaxies.This, combined with the fact that our fits display good reduced χ 2 values, gives us confidence that the S/N in each native pixel is sufficient to provide trustworthy estimates.

SUMMARY AND CONCLUSIONS
We present the first spatially resolved analysis of spectroscopically confirmed 5 < z < 9 galaxies in the SMACS0723 ERO field.We use images in 6 bands obtained with NIRCam onboard JWST, spanning the wavelength range 0.8 − 5µm.With the SED fitting software Bagpipes, we model the spectral energy distributions on a pixel-by-pixel basis, being able to infer the physical parameters on a 180 − 240 parsec scale.Our main findings and conclusions are the following: • All galaxies are resolved and display strong empirical colour gradients.Even at these early times, these galaxies display multiple star forming clumps.
• We find regions that exhibit high EW of the [O III]+Hβ emission (up to ∼ 3000 − 4000 Å).These extreme starbursts are embedded within regions with less specific star formation, which points to very bursty star formation happening on small scales (< 1 kpc), not galaxy-wide.
• The strong line emission regions dominate the integrated light, biasing the fits towards very young inferred ages of the stellar population (< 10 Myr).
Only a resolved analysis demonstrates the presence of older stellar populations, which can be seen in the spatial maps.
• Resolving the stellar populations on a pixel-bypixel basis leads to inferring from 0.5 up to ∼ 1 dex larger stellar masses, when compared to an integrated analysis.Our analysis extends previous findings on the problem of outshining and its effects on stellar mass estimates, which so far has only been studied at lower redshifts (up to z ∼ 3).
Current and upcoming observations with JWST will allow us to characterise the early Universe and first galaxies in a new and more complete way.The combination of having confirmed redshifts with NIRSpec, and the unprecedented resolution and depth of NIRCam imaging, will transform how we study galaxies, changing our current views on their internal structure and mass assembly, amongst others.The systematics in stellar mass estimates found in this work would have strong implications in the shape and evolution of the stellar mass function at high redshift, particularly while samples are limited to small numbers of the brightest candidates.Furthermore, the process of galaxy formation could be more extended and earlier than previously thought, as is implied by the presence of evolved older stellar populations being outshone by the youngest stars.Only with a spatially resolved analysis, we can begin to untangle the complexity of the internal structure of galaxies at this epoch.3 shows the resulting physical parameters inferred with Bagpipes on the integrated fit for each galaxy.As a reminder, we expect these values to be lower than the ones inferred by other works, since we only consider the pixels here that fulfil a certain S/N criteria, instead of performing aperture photometry, which would yield larger fluxes and different physical estimates.The integrated run is performed like this to be able to do a one-to-one comparison with the spatially resolved analysis, as discussed in §4.2.

Figure 1 .
Figure1.Cutout images of the five SMACS0723 galaxies in all available NIRCam bands.The cutouts are 1. 2 across and centered in the coordinates provided in Table1.The RGB images are built combining the F150W (B), F277W (G) and F444W (R) PSF-matched images (to the F444W band).The physical scale is calculated in the lens plane.

ID8140Figure 2 .
Figure 2. Maps of the galaxy ID8140 at z = 5.275.Left: Maps of the empirical colours in AB mag, inferred as the difference of the bands F150 and F277W (top) and the images F356W and F444W (bottom).The contours correspond to the non PSFmatched F200W image.The physical scale is calculated in the lens plane.The FWHM of the F444W PSF is indicated.Right: Maps of the physical properties inferred with Bagpipes.The top row, from left to right, shows the resulting maps for the star formation rate surface density (SFRD), the logarithm of the stellar mass surface density (SMD), and the visual extinction AV .The bottom row shows the maps for the mass-weighted stellar age, the UV slope (β) and the equivalent width of the inferred [O III]+Hβ emission lines.No lensing correction has been applied.

ID5144Figure 3 .
Figure 3. Maps of the empirical colours and the physical properties inferred with Bagpipes on the galaxy ID5144 at z = 6.383.See Figure 2 for more details.

Figure 5 .Figure 6 .
Figure 5. Maps of the empirical colours and the physical properties inferred with Bagpipes on the galaxy ID6355 at z = 7.665.See Figure 2 for more details.The three boxes A, B, C indicate the pixels that are analysed in more detail in the text, as well as in Figure 7, where the best fit SEDs are shown for each individual pixel.

Figure 7 .
Figure7.Best fit SEDs for the three pixels A, B and C indicated in Figure5, from the galaxy ID6355 at z = 7.665.The turquoise points and errorbars correspond to the NIRCam photometry, the orange points to the best fit model, and the orange curve and shaded region is the best fit SED inferred with Bagpipes, and corresponding 16th and 84th percentile uncertainty interval.The inferred physical parameters are indicated for each pixel, as well as the reduced χ 2 of the fit.

Figure 8 .
Figure 8. Best fit SEDs and models for the integrated (black curve and circles) and resolved (red curve and squares) modelling of the five galaxies studied in this work.The turquoise points and errorbars correspond to the integrated NIRCam photometry.The red curve is inferred by summing the posterior distributions in all pixels, and calculating the 50th percentile of the resulting one.The shaded regions correspond to the 16th and 84th percentile of the summed posterior distribution.The inset cutouts correspond to the same RGB images as in Figure 1.The reduced χ 2 values of each fit are indicated.

Figure 9 .
Figure 9.Comparison between the stellar mass that we infer in the spatially resolved analysis, versus the integrated fit.No lensing correction is applied in any of the two estimates.The plotted values are shown Table2.The one-to-one line is indicated, as well as the 0.5 and 1 dex offset lines.

The
Figure 11.Mass weighted age 16th, 50th and 84th percentiles of the posterior distribution inferred with Bagpipes on the galaxy ID6355 at z = 7.665.

Table 2 .
Values for the stellar mass and mass weighted stellar age that we infer with Bagpipes, both in the integrated run and the spatially resolved analysis presented in this work.We plot the mass values in Figure9.
*Figure10.Comparison between the star formation history that we infer in the spatially resolved analysis (red curve), versus the integrated fit (black curve), for the galaxy ID10612 at z = 7.663.The shaded regions correspond to the 16-84th percentile range in each case.