We present the first statistical investigation of spatially resolved emission-line properties in a sample of 63 low-mass galaxies at $4\leq z<10$, using JWST/NIRSpec MSA data from the JWST Advanced Deep Extragalactic (JADES) survey focusing on deep, spatially resolved spectroscopy in the GOODS-S extragalactic field. By performing a stacking of the 2D spectra of the galaxies in our sample, we find an increasing or flat radial trend with increasing radius for [O iii]5007/H and a decreasing one for [Ne iii]3869/[O ii]3727 (3–4 $\sigma$ significance). These results are still valid when stacking the sample in two redshift bins (i.e., $4\leq z<5.5$ and $5.5\leq z<10$). The comparison with star-formation photoionization models suggests that the ionization parameter increases by $\sim 0.5$ dex with redshift. We find a tentative metallicity gradient that increases with radius (i.e., ‘inverted’) in both redshift bins. Moreover, our analysis reveals strong negative gradients for the equivalent width of H (7$\sigma$ significance). This trend persists even after removing known AGN candidates, therefore, it is consistent with a radial gradient primarily in stellar age and secondarily in metallicity. Taken all together, our results suggest that the sample is dominated by active central star formation, with possibly inverted metallicity gradients sustained by recent episodes of accretion of pristine gas or strong radial flows. Deeper observations and larger samples are needed to confirm these preliminary results and to validate our interpretation.

Spatially resolved emission lines in galaxies at $4\leq z<10$ from the JADES survey: evidence for enhanced central star formation

1 Introduction

The distribution of metals within galaxies is a consequence of the processes of metal production, circulation, and dilution, the latter primarily attributed to the phenomenon of accretion. Therefore, studying metallicity gradients is crucial for understanding how these physical processes work. Mostly due to selection effects, the majority of high-redshift galaxies ($z>4$) are likely caught in a complex and turbulent phase of their evolution, characterized by gas accretion, frequent merging events, and by gas outflows due to stars and/or active galactic nuclei (AGNs) feedback (see e.g. Bischetti et al., 2019, 2021; Tripodi et al., 2024b; Neeleman et al., 2021; Shao et al., 2019; Dayal & Ferrara, 2018). These mechanisms clearly affect the content and distribution of heavy elements in the galaxy’s interstellar medium (ISM, Davé et al., 2011). In the local Universe, the spatial distribution of metals in galaxies has been studied precisely, assessing the presence of radial variations and chemical enrichment levels across the galaxy. A large number of local star-forming (SF) galaxies exhibit negative radial gradients in their gas-phase metallicity (or simply metallicity), with the inner regions more chemically enriched with respect to the outskirts (e.g., Magrini et al., 2010; Berg et al., 2015; Belfiore et al., 2017; Bresolin et al., 2016; Li et al., 2018). These gradients are generally interpreted as being indicative of the so-called inside-out growth scenario of galaxy formation (Samland et al., 1997; Portinari & Chiosi, 1999; Gibson et al., 2013; Pezzulli et al., 2017). Deviations from negative gradients provide specific information on the evolutionary phase of galaxies. The observation of flattening gradients beyond a certain radius may result from radial mixing processes (e.g. outflow; see e.g. Choi et al., 2020), or (re)accretion of metal-enriched gas in the outer regions (Bresolin et al., 2012), or an ongoing merger (Kewley et al., 2010; Rupke et al., 2010a, b). Positive (or inverted, i.e. increasing with radius) gradients may indicate a high rate of accretion of pristine gas (Sánchez Almeida et al., 2018). Moving from the local Universe to cosmic noon, it has been shown that galaxies at $1<z<2$ also grew inside-out, similar to their local counterparts (Nelson et al., 2016; Suzuki et al., 2019; Tacchella et al., 2015, 2018; Wang et al., 2019, 2020, 2022) . However, metallicity gradients seem to evolve (see e.g., Li et al., 2022). Curti et al. (2020) investigated the metallicity gradients of a sample of 42 galaxies between $1.2<z<2.5$, finding that $\sim 85\%$ of their galaxies are characterized by metallicity gradients shallower than 0.05 dex kpc${}^{-1}$, and $\sim 89\%$ are consistent with a flat slope within $3\sigma$, suggesting a mild evolution with cosmic time. They also discovered three galaxies with inverted gradients, suggesting recent episodes of pristine gas accretion or strong radial flows. The advent of JWST enabled for the first time the study of metallicity gradients at $z>3\text{--}5$, with the same rest-frame optical diagnostics used for galaxies at $0<z<2$ (for studies of galaxies at $z<5$ see e.g., Troncoso et al., 2014). However, current studies are limited to a couple of bright, extended systems that are not representative of the galaxy population at $z>3$ (see studies performed using integral field unit (IFU) or NIRISS observations or modelling FIR emission lines from ALMA, e.g., Rodríguez Del Pino et al., 2023; Arribas et al., 2023; Venturi et al., 2024; Vallini et al., 2024; Wang et al., 2022). The metallicity gradients of the bulk of galaxies at redshifts higher than $z>4$ are still poorly understood. Being able to measure the metallicity at early epochs gives us the invaluable opportunity to study the history of the baryonic cycle and its influence on the evolution of galaxies. Moreover, at very high-z the interpretation of radial measurements of emission line ratios is more challenging than at low-z, since high-z galaxies show irregular shapes with several (merging) components and/or clumps. The collisionally excited [O iii]5007 and the recombination H optical emission lines are commonly used as diagnostics for the properties of gas in star-forming regions. In particular, H is driven primarily by ionizing radiation, while [O iii]5007 is sensitive to the gas-phase metallicity and the ionization parameter. In combination with other lines (e.g., [Ne iii]3869 and [O ii]3726,3729), these two diagnostics allow us to measure the metallicity, electron temperature, ionization parameter, and hardness of the radiation field in the ISM of a star-forming galaxy (Osterbrock, 1989). Indeed, a powerful way to trace metallicity gradients at high-z is by comparing the [O iii]5007/H emission line ratio with the ratio between [Ne iii]3869 and [O ii]3726,3729 (hereafter [O ii]3727) emission lines (e.g., Nagao et al., 2006; Levesque & Richardson, 2014). Neon is produced during the late evolutionary stages of massive stars and it is expected to closely track oxygen abundance (Thielemann et al., 1994; Henry & Worthey, 1999). Along with oxygen, neon is one of the principal coolants in H ii regions. Moreover, [Ne iii]3869 arises from a broader range of regions in the ionized nebula when compared to [O ii], given that [Ne iii]3869’s high critical density ($\log n[{\rm cm^{-3}}]=7$, see Appenzeller & Oestreicher 1988) makes its flux insensitive to the electron density even in high-density regions. The $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$ ratio has been proven a better diagnostic of ionization parameter than $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$, with a greater sensitivity at shorter wavelengths that accommodates more of the ionizing photons produced by young massive stars. It is also insensitive to reddening effects and usable as an empirical diagnostic of ionization parameter out to higher redshift than $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$, since [O iii]5007 is emitted at longer wavelength than [Ne iii]3869. The degeneracy between metallicity and ionization parameter has indeed proven challenging to disentangle when trying to calibrate abundance diagnostics for SF galaxies and H ii regions. Therefore, when investigating metallicity, it is indispensable to have a tracer that is also able to constrain the ionization parameter (e.g., McGaugh, 1991; Kewley & Dopita, 2002). Additionally, H is particularly powerful in probing ionized regions around young, massive stars. Leaving dust extinction aside, the intensity of the H line is directly linked to the ionizing photons emitted by these stars, providing a valuable metric for quantifying the ongoing star-forming activity. As a SFR diagnostic, H, like all the Balmer lines, inherits the same strengths and weaknesses of H: it is equally sensitive to variations in the IMF and to absorption of Lyman-continuum photons by dust within star-forming regions (Moustakas et al., 2006). More precisely, assuming a star-forming origin, the equivalent width (EW) of H is a tracer of the specific SFR (sSFR) in star-forming galaxies, modulo the absorption of ionizing photons by dust in HII regions. In this work, we aim to investigate the emission-line gradients in a sample of 63 galaxies at $4\leq z<10$ from the JADES survey. Specifically, we focus on spatially resolved analyses of the [O iii]5007/H and $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$ ratios to investigate the distribution of metallicity in high-z galaxies, and of the EW${}_{\rm\text{H\textbeta}}$ as a tracer of the sSFR. We focus on the EW and these specific line ratios as, given that they are very close in wavelength, they do not suffer severely of the wavelength dependent PSF. Although it would be tempting to consider line ratios involving other lines more widely separated in wavelength, the wavelength dependent PSF results in blending different regions at different wavelength and different slit losses. Selecting diagnostics consisting in lines that are close in wavelength also mitigate dust reddening effects. This is the first attempt to perform spatially resolved metallicity studies at $z>4$ in a large sample of galaxies. The paper is structured as follows. In Sect. 2, we present the observations and the data reduction; in Sect. 3, we report the analysis and results on stacked samples; in Sect. 4, we discuss the implications of our findings, comparing our results with observations at different redshifts and with photoionization models. Throughout the paper, we adopt a $\Lambda$CDM cosmology from Planck Collaboration et al. (2020): $H_{0}=67.4\ \rm km\ s^{-1}\ Mpc^{-1}$, $\Omega_{m}=0.315$ and $\Omega_{\Lambda}=0.685$. We assume everywhere a Chabrier initial mass function (Chabrier, 2003). All physical distances are proper distances.

2 Observations and data reduction

We use publicly available JWST/NIRSpec data from the JWST Advanced Deep Extragalactic Survey (JADES, Eisenstein et al., 2023b; Bunker et al., 2023; Rieke & Jades Collaboration, 2023), a collaboration between the JWST/NIRCam and NIRSpec GTO teams. In this work, we focus on deep, spatially resolved spectroscopy in the GOODS-S extragalactic field (Giavalisco et al., 2004). Specifically, we analyse 63 galaxies at $4\leq z<10$. These data were obtained from Programme ID 1210 (PI: N. Lützgendorf; henceforth: PID), and used the NIRSpec micro-shutter assembly (MSA; Jakobsen et al., 2022; Ferruit et al., 2022). For the galaxy images, we used publicly available NIRCam data from JADES itself (PID 1180; PI: D. Eisenstein), from the JWST Extragalactic Medium-band Survey (JEMS, PID 1963, PIs: C. C. Williams, S. Tacchella and M. Maseda; Williams et al., 2023), and from the First Reionization Epoch Spectroscopic COmplete Survey (FRESCO, PID 1895; PI: P. Oesch; Oesch et al., 2023). Even though JADES includes observations with a range of dispersers (Bunker et al., 2023; Carniani et al., 2023), in this work we use only the PRISM/CLEAR observations, because they are the deepest available while still providing sufficient spectral resolution for our purpose of separating [O ii]3727 from [Ne iii]3869, and H from [O iii]4959,5007 at redshifts $z>4$. These data were observed using a 3-shutter slit, with 1-shutter nodding to provide accurate background subtraction. The observations also included dithering to explore different regions of the detector and to safeguard against ‘disobedient’ shutters, i.e., shutters that are in a different open or closed state than the requested configuration. The observations consist of three pointings, where each target is allocated to up to three pointings according to its priority (as described in Bunker et al., 2023). The optimisation of the MSA allocation was performed using the EMPT software (Bonaventura et al., 2023). For our targets, the exposure times range from 9.3 to 28 hours. The data reduction is described in Bunker et al. (2023); we report here only the most relevant steps. The data reduction pipeline is based on the ESA NIRSpec Science Operations Team pipeline (Alves de Oliveira et al., 2018; Ferruit et al., 2022), and will be described in Carniani et al. (in prep.). Background subtraction was performed by using the local background from adjacent shutters, but extended sources and sources with contaminants were pre-identified by visual inspection and were not self subtracted, by using only empty shutters for the background. Some degree of self subtraction may still be present for the most extended sources (i.e., those extending more than two shutters), but these are only relevant at the lowest redshifts and do not affect the sample we use in this paper. The wavelength calibration includes a correction for intra-shutter target position. The pipeline also applies a wavelength-dependent path-loss correction, to account for flux falling outside the micro-shutters. This takes into account the wavelength-dependent size of the NIRSpec PSF as well as the intra-shutter target position, assuming a point-source light distribution. While this assumption is incorrect at lower redshifts, it is an excellent approximation for the redshift range used in this paper ($z>4$). Because the nominal spectral resolution of the prism varies by an order of magnitude ($R=30\text{--}330$ over the wavelength range, Jakobsen et al., 2022), the spectra are binned on an irregular grid, with varying pixel size matching the spectral resolution and ensuring Nyquist sampling of the resolution element.

3 Analysis and Results

3.1 Line fitting procedure

The main goal of this work is to study the line properties of $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}},\text{H\textbeta},\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}$ and [O ii]3727. Therefore, in this section, we briefly explain the procedure adopted for fitting the 2D spectra of our targets. We fit separately in pairs [O iii]5007 and H, and [O ii]3727 and [Ne iii]3869, with their corresponding underlying continua. We model the emission lines with Gaussian functions and the two underlying continua with 1^st-order polynomials. The doublet [O iii]4959,5007 has been fitted fixing the ratio between the peak fluxes (peak${}_{\text{[{O}\,{\sc{iii}}]{\textlambda}{4959}{}{}}}/{\rm peak}_{\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}}=0.335$) and the wavelength separation ($\Delta\lambda=47.94$ Å) of the two emission lines, and using the same FWHM for both emission lines. Similarly, we fit simultaneously the [Ne iii]3869 and [Ne iii]3967 doublet, adopting a ratio of 0.301 between the latter and the former, and a rest-frame wavelength separation of $98.73$ Å, and considering the same FWHM for both emission lines. The doublet [O ii]3727 is always blended given the resolution of our spectra, so we fitted it with a single Gaussian. This reduces to 3 the number of free parameters for both the [O iii] and the [Ne iii]3869 doublets. Therefore, we have 8 free parameters for each pair of lines (i.e., peak flux, peak wavelength, FWHM for [O iii]5007 and same for H, slope and intercept for their underlying continuum; same for [O ii]3727 and [Ne iii]3869). We explore the 8-dimensional parameter space for each pair of lines using a Markov chain Monte Carlo (MCMC) algorithm implemented in the EMCEE package (Foreman-Mackey et al., 2013), assuming uniform priors for the fitting parameters, considering 10 walkers per parameter and 1500 trials (the typical burn-in phase is $\sim$ 500 trials). Analogously, we fit the [O iii]5007, H, [O ii]3727 and [Ne iii]3869 emission lines and corresponding continua in the 5-pixel boxcar extraction of the 1D spectrum. Finally, we compute the integrated fluxes by integrating the best-fitting functions for each emission line. Moreover, regarding the 2D spectra, we measure the EW of H for each trial in the chain. Then, we derive the best-fitting value from the 50^th percentile of the EW${}_{\rm\text{H\textbeta}}$’s chain, and the error on EW${}_{\rm\text{H\textbeta}}$ from the 16^th and 84^th percentiles.

3.2 Line blending and contamination correction

When measuring [Ne iii]3869, we also include flux from H ($\lambda=3890.17$ Å) and He i 3889, both of which are blended with [Ne iii]3869 at the resolution of the prism. To estimate the contamination, we proceed as follows. For H, we calculate the ratio H/([Ne iii]3869 + H + He i 3889) from the stacked 1D spectrum (see Sect. 3.3), finding a value of $0.35\pm 0.01$. This means that the ratio H/([Ne iii]3869 + H + He i 3889) will be less than 0.14, assuming the Balmer ratios from Case B recombination, $T_{\mathrm{e}}=10,000$ K and $n_{\mathrm{e}}=100\leavevmode\nobreak\ \mathrm{cm}^{-3}$ (Storey & Hummer, 1995). We apply no reddening correction to H/[Ne iii]3869, therefore this contamination fraction is strictly an upper limit. For He i 3889, the estimate is much more uncertain, because this line can be optically thick (e.g., Robbins, 1968). As a tentative estimate, we use the prism spectra to measure the He i 5877 line and then use models to infer the He i 3889 contamination to [Ne iii]3869. For $n_{\mathrm{e}}=100\leavevmode\nobreak\ \mathrm{cm}^{-3}$ and temperatures $T_{\mathrm{e}}=5,000$, 10,000 and 20,000 K, Benjamin et al. (1999) gives He i 3889/He i 5877 of 0.6, 0.8 and 1.1, respectively (these are intrinsic values, without reddening correction). From the stacked 1D spectrum, we also measure a ratio of $\text{{He}\,{\sc{i}}{\,\textlambda}{5877}{}{}}/\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}=0.23\pm 0.01$, thus even assuming the highest He i 3889/He i 5877 value of 1.1, the overall contamination should be of order 0.25, similar to what we estimated for H. Thus we have both $\text{H\textzeta}/(\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}+\text{H\textzeta}+\text{{He}\,{\sc{i}}{\,\textlambda}{3889}{}{}})<0.14$ and $\text{{He}\,{\sc{i}}{\,\textlambda}{3889}{}{}}/(\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}+\text{H\textzeta}+\text{{He}\,{\sc{i}}{\,\textlambda}{3889}{}{}})<0.25$. Adding these two constraints, we obtain $(\text{H\textzeta}+\text{{He}\,{\sc{i}}{\,\textlambda}{3889}{}{}})/(\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}+\text{H\textzeta}+\text{{He}\,{\sc{i}}{\,\textlambda}{3889}{}{}})<0.39$, implying $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/(\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}+\text{H\textzeta}+\text{{He}\,{\sc{i}}{\,\textlambda}{3889}{}{}})>0.61$. This is an upper limit on the contamination, because, as we noted, we did not apply a reddening correction for H/H or He i 3889/He i 5877. In addition, it has to be noted that He i models are still quite uncertain, even for optically thin lines (Benjamin et al., 1999) – let alone the optically-thick He i 3889. In light of these difficulties, we do not apply the estimated correction. However, we show its magnitude in the relevant figures.

Figure 1: (a): 3-arcsec false-colour cutout. The nominal location of the shutters is overlaid; the MSA acquisition accuracy is better than one NIRCam pixel (0.03 arcsec); the gray circle represents the FWHM NIRSpec PSF from webbpsf, at the observed wavelength of [O iii]5007. (b): 2D spectrum of 4270 at $z=4.023$. The y axis has been shifted so that the centre of the galaxy is at pixel 0. The physical scale of one pixel is reported on the plot (0.1 arcsec/pixel). (c): 1D spectrum of 4270. Note the complex morphology of this galaxy: the synergy between spectroscopy and imaging is vital for interpreting correctly spatially resolved slit spectroscopy. The simultaneous presence of high-ionisation and low-ionisation species, including auroral lines (cf. He ii 4686 vs. [N ii]5755, [O i]6300,6363, [O ii]7319–7332, possibly [N i]5198,5200) may indicate a complex interplay between AGN photoionization and shocks.

Figure 2: Same as Fig. 1 for 16745 at $z=5.566$.

3.3 Line ratios: stacking analysis

We aim to perform a spatially resolved study of the [O iii]5007/H and [Ne iii]3869/[O ii]3727 line ratios in a sample of galaxies at $4\leq z<10$. Therefore, we initially select all the publicly available galaxies belonging to the JADES survey with S/N$\geqslant 5$ on the flux of each of the four spectral lines at $4\leq z<10$: we find that 19 objects satisfy these criteria. In Figs. 1 and 2, we show 3-arcsec false-colour cutouts and spectra of the two galaxies with the most extended and bright [O iii]5007, H, [Ne iii]3869, and [O ii]3727 emissions. We extract the spectra from 7 pixels (one spectrum for each pixel’s row from -3 to +3, with pixel scale = 0.1 arcsec/pixel; see e.g., 2D spectrum in Fig. 1) centred on the peak flux of the emission lines of each galaxy. We do not find any galaxy showing extended emission beyond pixel +3/-3. Considering that the spatial resolution of our data is about 2 pixels (see Appendix A for estimate of the PSF size), we are not able to resolve emission within 1 pixel from the galaxy centre, that is within pixel -1 and 1. We obtain detections or upper/lower limits on the [O iii]5007/H and/or [Ne iii]3869/[O ii]3727 line ratios with a separation of more than 1 pixel from the galaxy centre only for 5 sources. They show flat radial trends of the line ratios. However, the complex morphology shown in the NIRCam images makes the interpretation of these tentative gradients arduous. Details on the analysis and results of individual sources can be found in Appendix B.

Given the tentative results on the 19 individual sources, to increase the S/N of the emission lines, we stack the 2D and 1D spectra of all the sources at $4\leq z<10$ that have S/N$\geqslant 5$ on [O iii]5007. This sample consists of 63 galaxies. To avoid being biased by the brightest objects in the sample, we start by dividing each spectrum by the integrated [O iii]5007 flux (for both the 1D and 2D analyses). We then weight the normalized spectra by the inverse of the variance. In addition to the full stack at $4<z<10$, we also consider two separate subsets, at $4\leq z<5.5$ ($z_{\rm med}=4.7$) and at $5.5\leq z<10$ ($z_{\rm med}=6$). These redshift bins are chosen to ensure a comparable number of objects in the two bins (30 and 33, respectively), and to cover time intervals of similar duration (0.50 and 0.57 Gyr, respectively). As an example, in Fig. 3, we present the stacked spectra extracted from pixel 0 (top panel) and pixel 2 (bottom panel, corresponding to an average radial distance of $\sim 1.2$ kpc) of the 2D full stacked spectrum. The fit to the emission lines of interest has been performed as explained in Sect. 3.1.

3.4 Radial gradients of emission-line ratios

The averaged radial profiles of line ratios in the three samples are shown in Fig. 4: the errors on the averaged results account for both the statistical error of each measurement and the standard deviation between the measurements at each radius. We use a simple linear regression to quantify the gradient of these line ratios. The slopes for each sample and line ratio are reported in Tab. 1. While $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$ shows a decreasing trend towards the outskirts in all the samples, the $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{H\textbeta}$ ratio is almost constant with radius. The measured gradients are only tentative, with a statistical significance between 3 and 4 $\sigma$; to confirm a real trend we would need deeper observations that allow to reduce the uncertainties and possibly to extend the detections to larger radii. We also note that the average value of both line ratios varies with redshift comparing the two redshift-bin sub-samples.

To better understand the implication of our results, in Fig. 5 we compare them to previous measurements at different redshifts found in the literature, and to photo-ionization models. In each of the three panels, the line ratios from the 1D spectrum are intermediate between the spatially resolved line ratios. Overall, we see radial trends in the three samples, which may be caused by a combination of varying ionization parameter and/or metallicity. In particular, for the $4\leq z<10$ stack, the radial trend appears to be due mostly to variations in the ionization parameter. However, when we divide the sample at $z=5.5$, we find that the two redshift bins have different average ionization parameter, and the radial variation also shows a component consistent with positive metallicity gradients. This is a direct consequence of the different average value of line ratios seen in Fig. 4. This will be extensively discussed in Sect. 4.

Table 1: Radial gradient of $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{H\textbeta}$ and $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$ ratios, EW${}_{\rm\text{H\textbeta}}$ and metallicity for the three stacked samples.

Sample	Median redshift	% of galaxies	$\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{H\textbeta}$	$\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$	EW${}_{\rm\text{H\textbeta}}$	$\log(Z/Z_{\odot})$
			[1/pixel]	[1/pixel]	[Å/pixel]	[dex/pixel]
$4\leq z<10$	5.5	63	0.024 $\pm$ 0.008	-0.083 $\pm$ 0.021	-23 $\pm$ 2	0.08 $\pm$ 0.05
$4\leq z<5.5$	4.7	30	0.013 $\pm$ 0.014	-0.085 $\pm$ 0.030	-20 $\pm$ 2	0.07 $\pm$ 0.06
$5.5\leq z<10$	6.0	33	0.023 $\pm$ 0.007	-0.120 $\pm$ 0.033	-35 $\pm$ 5	0.04 $\pm$ 0.05

Figure 3: Stacked spectra using the NIRSpec/MSA prism spectra of galaxies at $4\leq z<10$ belonging to the JADES survey. Top panel: spectrum extracted from the central pixel of the 2D stacked spectrum. The cyan (violet) line is the best-fit function for the continuum and [O iii]5007, H ($\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}},\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}$) spectral lines. Bottom panel: spectrum extracted from two pixel off central (i.e., $\sim 1.2$ kpc radius) of the 2D stacked spectrum. Colour coding and labels are the same as top panel. Even from a visual inspection of these two panels, one can notice the decrease of the [Ne iii]3869/[O ii]3727 ratio towards larger radii, i.e., from pixel 0 ([Ne iii]3869/[O ii]3727 $\sim 1$) to pixel 2 ([Ne iii]3869/[O ii]3727 $<1$).

Figure 4: Averaged radial profiles of line ratios from the 2D stacked spectra at different redshift ranges. Cyan (violet) dots mark the logarithm of [O iii]5007/H ($\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$) ratio as a function of radius expressed in pixels (pixel 0 is the centre of the galaxy). The left panel show the results stacking all galaxies at $4\leq z<10$; the central panel stacking all galaxies at $4\leq z<5.5$; the right panel stacking all galaxies at $5.5\leq z<10$. The fit for the [O iii]5007/H ($\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$) profile is shown as a dashed (dash-dotted) line with shadowed region ($3\sigma$ confidence interval). The orange vertical lines mark the spatial resolution of the observations given by the 50th (solid), 16th (dashed left) and 84th (dashed right) percentiles of the distribution of the PSF size, computed from the redshift distribution of galaxies in each sample. For each panel, top x-axis is computed considering the median redshift of each sample that is: $z=5.5$ for the total stack, $z=4.7$ for the low-z stack, and $z=6$ for the high-z stack. These correspond to a scale of 6.113 kpc/”, 6.614 kpc/” and 5.832 kpc/”, respectively.

Figure 5: [O iii]5007/H–[Ne iii]3869/[O ii]3727 line ratio diagram for the stacked samples. For comparison, we show $z\sim 0$ SDSS AGNs (galaxies) as blue (pink) colormap with contours, $z\sim 2$ MOSDEF galaxies and AGNs (black contours), SMACS 06355, 10612 and 04590 red squares; Trussler et al., 2023; the left-most square of the three is 06355, the type-II AGN identified by Brinchmann, 2023, the type-I AGN host GS 3073 at $z=5.55$ (filled and hollow pentagon, the latter estimating the flux of [Ne iii]3869 based on the Case B assumption modulated by the median dust attenuation; Ji X. in prep.), and the $z=7.15$ AGN type-I’s host galaxy ZS7 yellow cross and plus, depending on whether line fluxes are computed from the BLR location or [O iii] centroid, respectively; see Übler et al., 2023 for details. Spatially resolved ratios for our stacked samples of galaxies at $4\leq z<10$, at $4\leq z<5.5$, and at $5.5\leq z<10$ are plotted as solid circles in all the three panels and are colour and size coded based on the distance (in pixel) from the centre of the galaxy (i.e., pixel 0). Ratios for 1D stacked spectra of our samples are plotted as green rhombs. Overlaid are the star-formation photoionization models of Gutkin et al. (2016) at hydrogen densities $\log n{\rm[cm^{-3}]}=2.0$. The dotted and solid coloured lines show the variation of the ionization parameter at fixed metallicity and the variation of metallicity at fixed ionization parameter, respectively (color scales on the right-hand side of the figure). With the spectral resolution of the prism, we measure a blend of [Ne iii]3869, He i 3889 and H; the red horizontal arrow at the bottom of each panel represents the maximum correction for [Ne iii]3869/([Ne iii]3869 +He i 3889 +H) of $\approx$0.2 dex (see Sect. 3.2). This correction applies only to our data, and to GS 3073; all other data in the figure have sufficient spectral resolution to deblend [Ne iii]3869.

3.5 Analysis of H equivalent width

Fig. 6 presents the results on the radial profiles of EW${}_{\rm\text{H\textbeta}}$ for the three samples (full stacked and two z-bins). Overall EW${}_{\rm\text{H\textbeta}}$ shows a decreasing trend with radius, implying that the emission of H is centrally concentrated relative to the continuum, and it is highest in the high-z sample. We use a linear regression to quantify the gradient of EW${}_{\rm\text{H\textbeta}}$, and the results for the slopes derived for each sample are reported in Tab. 1. A similar trend is also found in the radial profile of 4270 at $z=4.023$, shown in Fig. 13 (image and spectra of 4270 are shown in Fig. 2). In this case, we observe a negative radial gradient, steepest in the central regions, followed by a flattening at the outermost radii.

This rising trend of EW${}_{\rm\text{H\textbeta}}$ towards the center has some interesting implications on the nature of the innermost regions of these high-z galaxies, depending on the source of ionization. Assuming a star-forming origin, EW${}_{\rm\text{H\textbeta}}$ is a tracer of the sSFR in SF galaxies, which is directly connected to the age of the stellar population (and indirectly connected to its metallicity, given that a low-metallicity population shows higher EW). The right y-axis of Fig. 6 shows the scaling between the EW${}_{\rm\text{H\textbeta}}$ and the age at two fixed different metallicities ($\log(Z/Z_{\odot})=-1.0,-0.5$). This scaling has been computed using synthetic spectra calculated with Prospector (Johnson et al., 2021), with nebular emission treatments based on Cloudy (see Byler et al. 2017, for details). We assumed a delayed exponential star-formation history and constant metallicity, and defined the age as the look-back time when the star formation started (we tested different values of the e-folding time $\tau$ between 0.1 and 1 Gyr, but show the results only for a fiducial value of $\tau=0.25$ Gyr). The observed rising trend of EW${}_{\rm\text{H\textbeta}}$ indicates a negative radial gradient either in metallicity or in age. In principle, the same trends could be due to central emission boosted by an AGN or by a radial trend in the escape fraction of ionizing photons. These alternative scenarios will be better discussed in the following section.

Figure 6: Radial profile of EW${}_{\text{H\textbeta}}$ in the stacked samples. From left to right: total stacked sample, stacked sample of $z=4-5.5$ galaxies, stacked sample of $z=5.5-10$ galaxies. The fit for each sample is shown as a dashed line with shadowed region ($1\sigma$ confidence interval). At pixels 1 and 2, different markers are used to discriminate between pixels above (squares) and below (circles) the central pixel 0 in the slit (see also Fig. 2). Top x-axis as in Fig. 4.

Figure 7: Radial profile of the metallicity in the stacked samples. From left to right: total stacked sample, stacked sample of $z=4-5.5$ galaxies, stacked sample of $z=5.5-10$ galaxies. The fit is shown as a dashed line with shadowed region (1$\sigma$ confidence interval). Orange vertical lines and top x-axis as in Fig. 4.

4 Discussion

In the former sections, we presented the study of spatially resolved emission line properties in a sample of galaxies at $4\leq z<10$, by stacking the 2D and 1D spectra to get the radial profiles of the line ratios. We also divided this larger sample in two redshift bins ($4\leq z<5.5$ and $5.5\leq z<10$), and we performed the stacking for these two sub-samples. Fig. 4 and Tab. 1 show a slightly increasing or flat (for the $4\leq z<5.5$ sub-sample) trend for the $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{H\textbeta}$ ratio as a function of radius, while the $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$ ratio decreases with radius at more than $3\sigma$ confidence level for both the larger sample and the two sub-samples.

In Fig. 5, we compare our results with low-z data in the literature and with star-forming photoionization models. In particular, we show the line ratios measured in the local Universe, both in star-forming galaxies and AGNs (see SDSS DR7 data in pink and blue, Kauffmann et al., 2003; Brinchmann et al., 2004; Abazajian et al., 2004), in a sample of $\sim 1500$ galaxies and AGNs at $1.37<z<3.8$ (from the MOSDEF survey, black contours Kriek et al., 2015; Reddy et al., 2015), in 3 galaxies at $z>7$ (from the SMACS 0723 JWST Early Release Observation field, red sqaures, Pontoppidan et al., 2022; Trussler et al., 2023), in the type-I AGN host GS 3073 at $z=5.55$ (filled and hollow pentagon; Ji X. in prep.), and in the $z=7.15$ AGN type I’s host galaxy ZS7 (yellow cross and plus Übler et al., 2023). For SDSS, we include all the sources at $z>0$ and with S/N$>2$ on the emission lines of interest ([O iii]5007, H, [O ii]3727, [Ne iii]3869, H, [N ii]6584). It is worth mentioning that SDSS spectra investigate the nuclear region of the galaxies, since they are extracted from the central 3 arcsec of each galaxy. We discriminate between galaxies and AGNs using the BPT diagram (Baldwin et al., 1981) and considering conservatively as AGNs all the objects lying 0.2 dex above the Kewley et al. (2001) and Schawinski et al. (2007) lines, while as SF galaxies all the objects 0.2 dex below the Kauffmann et al. (2003) line. Overplotted are the photoionization models presented in Gutkin et al. (2016), at fixed hydrogen density $\log(n)=2\leavevmode\nobreak\ \mathrm{cm}^{-3}$. Increasing or decreasing the density by a factor of 10 does not alter the predicted line ratios. The variation of the ionization parameter and the metallicity along the models are highlighted according to the colormaps in the right legend.

We see two different behaviours for the full stack at $4\leq z<10$ and for the two redshift bins. For the former, the radial trend appears to be due mostly to variations in ionization parameter. However, when we divide the sample at $z=5.5$, the picture changes. First, we find that the two redshift bins have different average ionization parameter ($\Delta\sim 0.5$ dex) and metallicity ($\Delta\sim 0.4$ dex), with the highest-redshift bin having the highest average value of $\log\,U$ and the lowest value of $Z$, consistent with other studies (Cameron et al., 2023; Sanders et al., 2023; Curti et al., 2023b; Nakajima et al., 2023). Second, the direction of the radial variation of the observed line ratios is intermediate between grid lines of constant $\log\,U$ and constant metallicity. Note that, in the model grid, we fixed the dust-to-metal mass ratio $\xi_{\mathrm{d}}=0.3$. Radial variations in $\xi_{\mathrm{d}}$ would affect both the emission-line ratios we consider, because higher dust depletion would preferentially remove oxygen from the ISM gas (relative to hydrogen and neon). Specifically, the diagnostic diagram suggests a radial increase in gas metallicity, which could be explained by a radially decreasing $\xi_{\mathrm{d}}$ at constant metallicity of the ISM (i.e., gas and dust). However, a decreasing radial gradient of $\xi_{\mathrm{d}}$ would also increase EW${}_{\text{H\textbeta}}$, because a lower dust fraction would increase the flux ratio between the UV ionising photons (causing H emission) and the rest-optical photons in the continuum near H. Thus, while we cannot rule out a role for radial trends in $\xi_{\mathrm{d}}$, if such a trend was dominant, it would produce an increasing trend of EW${}_{\text{H\textbeta}}$ with radius, which is contrary to our observations.

For all the three samples, the radial gradients are in the direction of decreasing ionisation parameter and increasing metallicity, although the large measurement uncertainties prevent us from drawing strong conclusions. For the full stacked sample, the radial trend may be mostly due to the variation of ionization parameter. For the two z-bins samples, we find tentative inverted metallicity gradients, and an almost constant ionization parameter as a function of radius, with the only exception being the outermost radius of the higher-z sample for which the ionization parameter drops by $\sim 0.5$ dex. Moreover, the metallicity is on average higher ($\sim 0.4$ below solar) in the $z=4-5.5$ sub-sample than in the higher-z one ($\sim 0.8$ dex below solar), as one would expect for more evolved objects, and is intermediate between the two in the total stacked sample. We note, however, that this strong metallicity evolution is much larger than what is measured using e.g. strong-line methods (e.g. Curti et al., 2023b). It is therefore unclear whether the observed increase reflects a physical evolution, or arises from a sample-selection bias. In fact, the highest-z galaxies may have higher average SFR than the targets in the lowest-z bin (see e.g., Scoville et al., 2023; Tripodi et al., 2024a; Speagle et al., 2014; Calabrò et al., 2024; Genzel et al., 2010). The fundamental metallicity relation (FMR; Mannucci et al., 2010; Lara-López et al., 2010; Cresci et al., 2019; Baker et al., 2023a) predicts that at fixed stellar mass, metallicity decreases with increasing SFR, which could explain the observed trend. The recent observation of deviations from the FMR at $z>6$ (Curti et al., 2022, 2023b) would only reinforce this trend.

To investigate these inverted metallicity gradients further, we convert both the $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{H\textbeta}$ and $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$ ratios in O/H, a proxy for gas-phase metallicity, using the calibrations presented in Curti et al. (2020), and updated in Curti et al. (2023a, b). The derived metallicity radial profiles for the three stacked samples are shown in Fig. 7. Also in this case, the metallicity in the higher-z sub-sample is -on average- $\sim 0.3$ dex lower than in the $z=4-5.5$ sub-sample. We perform a fit to these profiles using a linear regression, and the results are reported in Tab. 1. The increasing trend with radius in all three subsets is suggestive of inverted metallicity gradients, but the large uncertainties on the individual measurements prevent us from drawing definitive conclusions (all three fits are statistically consistent with no gradient).

In Appendix B, we consider the radial gradients of five spatially resolved galaxies. However, interpreting the nature of these individual gradients is typically challenging. This complexity arises from the compact nature of galaxies within the considered mass range and redshifts. de Graaff et al. (2023) also found small sizes and complex kinematic features in six $5.5<z<7.4$ galaxies belonging to JADES, three of them showing significant spatial velocity gradient. As a result, individual gradients can be measured only for complex and possibly interacting systems. In this context, the synergy between photometry and spectroscopy is crucial to drawing the correct view about the nature of the gradients (see Fig. 8). In contrast, when stacking, we ‘marginalize’ over the precision and accuracy of the shutters alignment, making the radial gradients easier to interpret. Whether a sample of 63 galaxies is sufficient for marginalising over these effects, and enables us to infer the true average gradients of the population, would require simulated observations of realistic galaxies, but this is beyond the scope of the present paper.

As we have noted, the analysis of the EW of H reveals a steep radial profile increasing towards the center in the stacked spectra. Additionally, while analysing the individual galaxy 4270, which shows the clearest gradients among the brightest 5 galaxies in Appendix B, we found a similar trend for EW${}_{\text{H\textbeta}}$ vs radius. In SF galaxies, the EW of Balmer lines like H and H tracks the cosmic evolution of the sSFR (Fumagalli et al., 2012; Sobral et al., 2014). Indeed, for a fixed age (and, therefore, for a fixed stellar mass-to-light ratio), EW${}_{\text{H\textbeta}}$ can be considered as a rough proxy for sSFR. Therefore, the higher EW${}_{\text{H\textbeta}}$ in the center could be due to our sample being dominated by galaxies undergoing centrally concentrated star formation, leading to the build-up of stellar mass in the core of galaxies (e.g., Baker et al. 2023b; see Zhang et al. 2012 for a local example in a dwarf galaxy). This phase of ‘core formation’ is indeed expected from both numerical simulations and analytical models (Dekel et al., 2009; Krumholz et al., 2018; Tacchella et al., 2016; Tripodi et al., 2023; Zolotov et al., 2015) and observations of high-z galaxies (see e.g., Ji & Giavalisco, 2023; Ji et al., 2023), and is fuelled by rates of pristine/low-metallicity gas accretion that are much faster than in the local Universe. This interpretation would agree with the inverted metallicity gradients, for which the most natural explanation is the continuous dilution of otherwise metal-rich gas in the central regions due to the same gas inflows powering the central star formation. In principle, the presence of an AGN could also cause a central enhancement of EW${}_{\text{H\textbeta}}$, given that AGNs tend to have high-EW nebular emission (e.g., Caccianiga & Severgnini, 2011, for EW${}_{\rm\text{[{O}\,{\sc{iii}}]{}{}{}{}}}$). This may be a plausible scenario for 4270, because the overlap of its line ratios with those of local SDSS AGNs already suggested that 4270 may itself be an AGN host (see Fig. 12).

To ensure that the presence of AGNs did not bias our results on the stacked samples, we excluded from the stacking the 15 AGNs identified in our sample (Scholtz et al., 2023), and performed the analysis on the line ratios and EW${}_{\text{H\textbeta}}$ again. The results are unchanged both for the total and the two z-bin stacked samples. That is, even excluding the AGNs, we find flat radial profiles for [O iii]5007/H, and declining for [Ne iii]3869/[O ii]3727, and a steep increase of EW${}_{\text{H\textbeta}}$ towards the center. This suggests that AGNs do not explain the observed trends in EW${}_{\text{H\textbeta}}$, which, therefore, must be due to stellar-population trends.

It should also be noted that the redshift range we are probing is relatively uncharted for this topic, therefore our interpretation is uncertain. In particular, star-forming galaxies in burst phases are much more common in the young Universe (Ciesla et al., 2023; Dressler et al., 2023; Looser et al., 2023; Endsley et al., 2023), as also expected from models (see e.g., Faucher-Giguère, 2018; Tacchella et al., 2020). These objects have much higher EW${}_{\text{H\textbeta}}$ and EW${}_{\rm\text{[{O}\,{\sc{iii}}]{}{}{}{}}}$ than star-forming galaxies in the local Universe (e.g., Smit et al., 2014; van der Wel et al., 2011; Maseda et al., 2014; Boyett et al., 2024), therefore they may equal or even exceed the EWs of AGNs (for type-2 AGN, EW${}_{\text{[{O}\,{\sc{iii}}]{}{}{}{}}}<500$ Å; e.g., Caccianiga & Severgnini, 2011). In principle, the observed trends could be caused by the fact that, on average, EW${}_{\text{H\textbeta}}$ increases with z (e.g., Smit et al., 2014), whereas the typical size of galaxies decreases (van der Wel et al., 2014; Ji et al., 2024). However, we still observe the same trend between EW${}_{\text{H\textbeta}}$ and radius even after splitting our sample in two redshift bins, which suggests that the trend is real and not due to the observational conspiracy between EW${}_{\text{H\textbeta}}$, size and redshift. Alternatively, if starburst galaxies had both higher EW${}_{\text{H\textbeta}}$ and smaller size than non-starburst galaxies, this could also produce the decreasing radial gradients of EW${}_{\text{H\textbeta}}$ in the stack. This is because, under this hypothesis, starburst galaxies (with their high EW${}_{\text{H\textbeta}}$) would dominate the central pixels of the stack, while more extended, less bursty galaxies (with lower EW${}_{\text{H\textbeta}}$) would dominate the outer pixels of the stack, thus creating a spurious gradient. However, distinguishing between inverted gradients and correlations between galaxy structure and emission-line ratios (and EW${}_{\text{H\textbeta}}$) requires further dividing the sample; for this technique to be effective, we need a larger sample than the currently available one. NIRCam could help evaluate this possibility, combining SED fitting (to identify starbursts) and NIRCam medium- vs wide-band imaging (to measure EW gradients).

Alternatively, IFU observations with NIRSpec can directly measure spectral gradients in galaxies (Rodríguez Del Pino et al., 2023; Arribas et al., 2023; Venturi et al., 2024). However, NIRSpec/integral field spectroscopy (IFS) is best suited for large, bright galaxies, primarily due to its low multiplexity. Using this instrument configuration for targeting a large sample of low-mass, high-z galaxies with a depth comparable to JADES, would require an inordinate amount of time. This is compounded by the lower instrument sensitivity and, possibly, by a larger PSF at 1–3 µm (at least in the current cube reconstruction; e.g., D’Eugenio et al., 2023). For all these reasons, we deem it unlikely that NIRSpec/IFS can be competitive for building large samples – at least in the redshift and mass range explored in this article.

Figure 8: 3-arcsec false-colour cutouts of all the galaxies with high S/N and with available JWST/NIRCam photometry, sorted by increasing redshift $z$ (targets 4270 and 16745 are reported already in Figs. 1 and 2 and are not repeated here). The nominal location of the shutters is overlaid; the MSA acquisition accuracy is better than one NIRCam pixel (0.03 arcsec). This figure illustrates the key synergy between NIRCam and NIRSpec/MSA; imaging clarifies that for individual galaxies gradients along the slit may result from complex morphologies. However, because the intra-shutter positions are effectively random between different targets, a stacking analysis over a sufficiently large sample may ‘marginalise’ over peculiar morphologies.

5 Summary and Conclusions

We study the radial profiles of the $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{H\textbeta}$ and $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{\textlambda}{}{3726,3729}}$ line ratios, and of the EW${}_{\rm\text{H\textbeta}}$ in a sample of 63 galaxies at $4\leq z<10$, using publicly available JWST/NIRSpec data from JADES. We focus on deep, spatially resolved spectroscopy in the GOODS-S extragalactic field, since we want to study spatially resolved emission line gradients. Our findings are summarized as follows.

• By performing an initial selection of 19 galaxies with S/N$>5$ on the lines of interest, we find that 5 of them show flat radial trends of the line ratios. Specifically, for the two galaxies 4270 at $z=4.023$ and 16745 at $z=5.5666$, these trends are detected or marginally detected up to $\sim 2$ kpc from the center of the galaxy. These can be caused by merging events or re-accretion of metal enriched gas. However, the complex morphology shown in the NIRCam images makes the interpretation of gradients arduous. This highlights the importance of the synergy between imaging and spectroscopy for capturing the physical properties of these objects.

• To increase the S/N ratio on the lines of interest and to marginalize over the complex galaxy morphologies, we stack the whole sample of 63 galaxies and two sub-samples, created by dividing the larger sample in two redshift bins ($4\leq z<5.5$, and $5.5\leq z<10$). We found a slightly increasing or flat radial trend for $\text{[{O}\,{\sc{iii}}]{\textlambda}{5007}{}{}}/\text{H\textbeta}$, while a decreasing radial trend for $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{\textlambda}{}{3726,3729}}$ in all the three stacked samples.

• By comparing with photoionization models, we find that the two redshift bins have different average ionization parameters ($\Delta\sim 0.5$ dex, increasing with redshift), with tentative evidence of an inverted metallicity gradient. The metallicity is on average $\sim 0.3$ dex higher in the $z=4-5.5$ sub-sample, implying an evolution stronger than measured using e.g., strong-line methods. This is still valid if adopting the calibrations presented in Curti et al. (2020) to convert line ratios to metallicity. It remains difficult to discriminate whether this reflects a physical evolution or a selection bias, since galaxies with higher SFR at high-z are expected to have lower metallicity at fixed stellar mass.

• We also find negative radial gradients for EW${}_{\text{H\textbeta}}$, even when removing known AGNs from the sample. This may indicate that – on average – our sample of high-z galaxies is dominated by central starbursts, fuelled by high mass accretion rates. The latter could also explain the inverted metallicity gradients.

This work led to the first tentative detection of metallicity gradients and of negative gradient of EW${}_{\text{H\textbeta}}$ in a stacked sample of $4\leq z<10$ galaxies. There are several possible avenues to confirm and refine this study. For galaxies with high-EW emission lines, we could exploit the spatial resolution of NIRCam and the excellent wide- and medium-filter coverage of the JADES Origins Field (Eisenstein et al., 2023a). NIRCam slitless spectroscopy also offers promising possibilities e.g. FRESCO, Oesch et al., 2023; CONGRESS, (GO-3577), Sun et al. in prep. For fainter and/or low-EW targets, and especially regarding the study of metallicity gradients, we need larger samples of comparable depth. Finally, for the highest-mass galaxies, which are typically the most extended, we would require a different approach to the background subtraction (to avoid self subtraction) or to the observations e.g., the NIRSpec/IFS, Böker et al., 2022; Rodríguez Del Pino et al., 2023; Arribas et al., 2023.

Acknowledgment. RT and MB acknowledge support from the ERC Grant FIRSTLIGHT and Slovenian national research agency ARIS through grants N1-0238 and P1-0188. RT acknowledges financial support from the University of Trieste. RT acknowledges support from PRIN MIUR project “Black Hole winds and the Baryon Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper”, contract #2017PH3WAT. FDE, RM, JS, WB, XJ and JW acknowledge support by the Science and Technology Facilities Council (STFC), by the ERC through Advanced Grant 695671 “QUENCH”, and by the UKRI Frontier Research grant RISEandFALL. RM also acknowledges funding from a research professorship from the Royal Society. MC acknowledges the support from the ESO Fellowship Programme. AJB and AJC acknowledge funding from the “First Galaxies” Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 789056). JT acknowledges support by the ERC Advanced Investigator Grant EPOCHS (788113) from the European Research Council (ERC) (PI Conselice). SA acknowledges support from Grant PID2021-127718NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI/ 10.13039/501100011033). SC and GV acknowledge support by European Union’s HE ERC Starting Grant No. 101040227 - WINGS. ZJ and CNAW acknowledge funding from JWST/NIRCam contract to the University of Arizona NAS5-02015. BER acknowledges support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015, and JWST Program 3215. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. The authors acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with programme PID 1210 (NIRCam-NIRSpec galaxy assembly survey - GOODS-S - part #1b; PI N. Lützgendorf). Reduced and calibrated images and spectra were obtained from the JADES Collaboration, via the high-level science products page on MAST.

Facilities: JWST. Software: astropy (Astropy Collaboration et al., 2022), Matplotlib (Hunter, 2007), SciPy (Gommers et al., 2023), emcee (Foreman-Mackey et al., 2013).

Appendix A Assessing the spatial resolution of the data

A key requirement for making meaningful gradient measurements is that the data do contain spatial information. In this section, we show that the data do contain marginally resolved spatial information, by comparing the spatial profile of the [O iii]5007 nebular emission with the model point-spread function (PSF).

Figure 9: Example of model fit to the galaxy spatial profile, for galaxy ID 7892. We show the data (black), the model background (dotted red), the background plus the negative traces (blue), and the best-fit model (green). The dark green colour is due to overlapping blue and green. We report the observed wavelength of [O iii]5007 (bottom left) and the best-fit $\sigma$ (bottom right).

Figure 10: Comparison of the radial extent of our data (using [O iii]5007, circles) to the PSF model frow webbpsf (dashed line) and to the model from the data reduction pipeline (dotted line). The galaxies in our sample extend beyond the PSF $\sigma$.

For the galaxies, we measure the light profile along the slit by summing the 2D spectrum in a window of three wavelength pixels, centred on the [O iii]5007 wavelength. As a model, we used a linear background (two free parameters), a Gaussian representing the source profile (three free parameters), and four more Gaussians representing the negative source traces, which arise from the nod sky subtraction Sect. 2 and Bunker et al. 2023. These four negative Gaussians add nine free parameters in total, as we discuss below. The centroid of each of the four negative Gaussians is free (four out of nine parameters), but is subject to a Gaussian probability prior placing it one or two shutters away from the source (prior dispersion 0.2 arcsec). The flux of each of the four negative Gaussians is also free (four out of nine parameters), but again is subject to being smaller (in absolute value) than the flux of the source; in addition, we place a Gaussian probability prior such that the negative traces have flux which is within a fraction of the source flux (with prior dispersion 0.5). The final free parameter for the negative traces is their common dispersion, which is constrained to be identical to that of the source within a Gaussian probability prior with fractional dispersion 0.5. The simultaneous fitting of the negative traces is crucial for a reliable estimate of the galaxy’s profile in the regime where the PSF is mildly under-sampled. This model is optimised using again a MCMC integrator, and our fiducial value of the source extent is the $\sigma$ of the source Gaussian, for which we quote the 50^th percentile as fiducial value, and half the 16^th–84^th inter-percentile range as uncertainty. An example galaxy profile is reported in Fig. 9. The $\sigma$ of each galaxy in our sample is shown as a dot in Fig 10.

For the MSA PSF, we use the model from webbpsf, for which we report the Gaussian-equivalent $\sigma$, defined by $\sigma\equiv\mathrm{FWHM}/\sqrt{\ln{256}}$ (dashed line in Fig. 10). We tested that following the approach of de Graaff et al. (2023) gives consistent results.

We note that the process of co-adding the 2D spectra in the data reduction pipeline does not align the traces to the sub-pixel level. This is done to avoid resampling the data, which increases the correlation between the pixels. The pipeline currently aligns the 2D spectra from each visit by finding their peak. To model the effect of this procedure, we create mock 2D spectra with the wavelength-dependent $\sigma$ of webbpsf, but for each realisation, we sum 20 exposures with an added random offset of $\pm 0.5$ pixel, or 0.05 arcsec. We then fit the spatial profile of this distribution, and find that indeed its profile is broader than the model PSF (dotted line in Fig. 10), but still below the values measured for our data. By comparing the dotted line with the dashed line representing webbpsf, we estimate the loss in resolution from the alignment procedure to be equivalent to convolving the native PSF with a Gaussian PSF having $\sigma=0.03$ arcsec.

Appendix B Analysis of individual sources

Figure 11: Five best cases of radial profiles of line ratios. Cyan (violet) dots mark the logarithm of [O iii]5007/H ($\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$) ratio as a function of the radius expressed in pixels (pixel 0 is the centre of the galaxy). Different markers are used to discriminate between pixels above (squares) and below (circles) the central pixel 0 in the slit (see also Fig. 2). The ID of each object and the redshift are specified in the top left corner of each plot.

In order to perform a spatially resolved study of the line ratios [O iii]5007/H and [Ne iii]3869/[O ii]3727 in a sample of galaxies at $4\leq z<10$, we initially select all the publicly available galaxies belonging to the JADES survey with S/N$>=5$ on the flux of each spectral line at $4\leq z<10$: we find that 19 objects satisfy these criteria.

Fig. 1 shows the 2D spectrum of one of those objects (ID 4270 at $z=4.023$), where the y-axis has been shifted with respect to the brightest pixel along the slit, which represents the centre of the galaxy. We shift all the 2D spectra according to the centre of each galaxy. The [O iii]5007 emission of 4270 is very bright and clearly extended over 7 pixels, that is $\sim 5$ kpc at the galaxy redshift. We extract the spectra from these 7 pixels (from pixel -3 to +3). For this object, [Ne iii]3869 is blended with [Ne iii]3967, and we take this into account when fitting the spectrum (see Sect. 3.1)

Another interesting object is 16745, which also has a very bright [O iii]5007 emission, extended over 7 pixels that corresponds to $\sim 4$ kpc at the galaxy redshift (see Figs. 2 and 2 for 2D and 1D spectra). de Graaff et al. (2023) found disc-like kinematics in 16745, with $v/\sigma\sim 2$. We fitted both pairs of emission line as above, considering that in this case there is no strong blending between the [Ne iii]3869,3967 doublet.

Top left and top right panels of Fig. 11 present the radial profile of the ratios between the integrated fluxes of [O iii]5007 and H (cyan dots), and between [Ne iii]3869 and [O ii]3727 (violet dots) for 4270 and 16745 along with other three bright sources. Since the choice of the pixel sign is arbitrary, we show the fluxes as a function of the absolute value of the radius. The interpretation of these profiles is challenging, since the morphology of these sources is complex (see Figs. 1 and 2), probably caused by the presence of a companion or satellite galaxy.

We apply the same method for each object in this sample, always taking into account the blending of the [O iii]4959,5007 doublet, and the blending of the [Ne iii]3869,3967 doublet. The results for the other 4 best cases are shown in the other panels of Fig. 11, where the object ID and redshift are specified in the top left corner of each plot. We select the cases in which the [O iii]5007/H ratio is detected up to pixel 2 at least, and [Ne iii]3869/[O ii]3727 up to pixel 1 at least. The other sources present upper/lower-limits or non-detections in the off-center pixels, especially for the $\text{[{Ne}\,{\sc{iii}}]{\textlambda}{3869}{}{}}/\text{[{O}\,{\sc{ii}}]{\textlambda}{3727}{}{}}$ ratio.

Overall, the most extended target is found to be 4270, which shows emission of [Ne iii]3869 and [O ii]3727 up to pixel 3 (i.e., $\sim 2$ kpc from the galaxy centre). Fig. 12 compares our emission-line ratios for the 2D and 1D spectra to the literature. When comparing with low-z data in the literature and with the star-forming photoionization models of Gutkin et al. (2016), we see that 4270 is consistent with the redshift evolution of the [O iii]5007/H and [Ne iii]3869/[O ii]3727 ratios, which both increase from $z=0$ to $z\sim 2$ and $z\sim 7$. In particular, the measurements for this target overlap with the high-[Ne iii]3869/[O ii]3727 envelope of the SDSS AGNs. The spatially resolved measurements have large uncertainties, yet they show some tentative dependencies on metallicity and ionization parameter. We can see that the ratios of the 1D spectrum are averages between the resolved measurements, which is a reassuring consistency check, implying that the 1D spectrum captures the average characteristics of the galaxy.

Figure 12: [O iii]5007/H–[Ne iii]3869/[O ii]3727 line ratio diagram for ID 4270. Ratios for 4270 at $z=4.023$ are plotted as solid circles, that are colour and size coded based on the distance from the centre of the galaxy (i.e., pixel 0 is the lightest and biggest). For comparison, we show $z\sim 0$ SDSS AGNs (galaxies) as blue (pink) colormap with contours, $z\sim 2$ MOSDEF galaxies and AGNs (black contours), SMACS 06355, 10612 and 04590 red squares; Trussler et al. 2023; the left-most square of the three is 06355, the type-II AGN identified by Brinchmann 2023, the type-I AGN host GS 3073 at $z=5.55$ (filled and hollow pentagon; Ji X. in prep.), and the $z=7.15$ AGN type-I’s host galaxy ZS7 yellow cross and plus, depending on whether line fluxes are computed from the BLR location or [O iii] centroid, respectively; see Übler et al. 2023 for details. Spatially resolved ratios for our stacked samples of galaxies at $4\leq z<10$, at $4\leq z<5.5$, and at $5.5\leq z<10$ are plotted as solid circles in all the three panels and are colour and size coded based on the distance (in pixel) from the centre of the galaxy (i.e., pixel 0). Ratios for 1D stacked spectra of our samples are plotted as green rhombs. Overlaid are the star-formation photoionization models of Gutkin et al. (2016) at hydrogen densities $\log n{\rm[cm^{-3}]}=2.0$. The dotted and solid coloured lines show the variation of the ionization parameter at fixed metallicity and the variation of metallicity at fixed ionization parameter, respectively (color scales on the right-hand side of the figure). With the spectral resolution of the prism, we measure a blend of [Ne iii]3869, He i 3889 and H; the red horizontal arrow at the bottom of each panel represents the maximum correction for [Ne iii]3869/([Ne iii]3869 +He i 3889 +H) of $\approx$0.2 dex (see Sect. 3.2). This correction applies only to our data, and to GS 3073; all other data in the figure have sufficient spectral resolution to deblend [Ne iii]3869.

The flat shape of these radial profiles can be caused by merging events or re-accretion of metal-enriched gas (Oppenheimer & Davé 2008; Übler et al. 2014). Specifically for the cases of 4270 at $z=4.023$ and 16745 at $z=5.566$, both these scenarios may be likely. This is because both sources have a complex morphology (Figs. 1 and 2) and it is unclear whether we are observing satellite galaxies or simply a clumpy morphology. The position of the shutters is such that we may trace some part of the companion emission, and the central pixel (pixel 0) is actually tracing part of the central emission, more than the actual centre of the source. Moreover, note that the position of the shutters is known with finite precision. As a word of caution, complex morphology and slit position make the interpretation of gradients in individual sources really challenging, especially when the position of the shutters is not aligned with the source. As we have seen, the synergy between slit spectroscopy and imaging is of critical importance to grasp the complex morphology and physical properties of high-redshift galaxies; 2D maps of the stellar mass and/or gas kinematics are critical to understanding the structure of individual sources, and to correctly interpreting the observed spatial variations. In Fig. 12, we compare the results for 4270 with other galaxies and AGNs at different redshifts from the literature. Interestingly, the line ratios of 4270 overlap with the SDSS AGNs, suggesting that this object may host an AGN. Applying the H +He i 3889 contamination correction enhances the overlap with SDSS AGN, and would not change our conclusions. When comparing with photoionization models, at fixed hydrogen density ($\log n{\rm[cm^{-3}]}=2.0$), the metallicity is almost solar at the galaxy centre and increases up to 0.3 dex with radius. This inverted gradient can be explained by recent episodes of pristine gas accretion to the galaxy centre or strong radial flows, and it has been previously found in galaxies at lower redshift, $1.2<z<2.5$ (Curti et al. 2020). There also seems to be a trend with ionization parameter, with higher ionization towards larger radii, although the measurement uncertainties are large. Indeed, since the major photoionization sources (e.g., stars, AGN) are typically concentrated in the centre, this scenario of higher ionization in the outskirts appears unlikely, even if considering the presence of a companion. Therefore, we conservatively conclude that the ionization parameter radial profile is consistent with flat, given the large errorbars along the x-axis. For this target, we also find that EW${}_{\text{H\textbeta}}$ decreases with radius (Fig. 13), which would be in agreement with both the AGN scenario as well as a central starburst, fuelled by inflow of relatively metal-poor gas.

Figure 13: Radial profile of H EW in 4270 at $z=4.023$. At pixels 1 and 2, different markers are used to discriminate between pixels above (squares) and below (circles) the central pixel 0 in the slit (see also Fig. 2).