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

Roberta Tripodi, Francesco d'Eugenio, Roberto Maiolino, Mirko Curti, Jan

Scholtz, Sandro Tacchella, Cosimo Marconcini, Andrew J. Bunker, James A.

A. Trussler, Alex J. Cameron, et al.

► To cite this version:

Roberta Tripodi, Francesco d'Eugenio, Roberto Maiolino, Mirko Curti, Jan Scholtz, et al. Spatially resolved emission lines in galaxies at $4 \le z < 10$ from the JADES survey: Evidence for enhanced central star formation. Astronomy & Astrophysics - A&A, 2024, 692, 10.1051/0004-6361/202449980. insu-04851484

HAL Id: insu-04851484

https://insu.hal.science/insu-04851484v1

Submitted on 21 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spatially resolved emission lines in galaxies at $4 \le z < 10$ from the JADES survey: Evidence for enhanced central star formation

Roberta Tripodi<sup>1,2,3,4,5,6,*</sup> (D), Francesco D'Eugenio<sup>5,6,D</sup>, Roberto Maiolino<sup>5,6,7</sup>, Mirko Curti<sup>8</sup>, Jan Scholtz<sup>5,6,D</sup>, Sandro Tacchella<sup>5,6</sup>, Cosimo Marconcini<sup>9,10</sup>, Andrew J. Bunker<sup>11</sup>, James A. A. Trussler<sup>12</sup>, Alex J. Cameron<sup>11</sup>, Santiago Arribas<sup>13,D</sup>, William M. Baker<sup>5,6</sup>, Maruša Bradač<sup>1,D</sup>, Stefano Carniani<sup>14,D</sup>, Stéphane Charlot<sup>15,D</sup>, Xihan Ji<sup>5,6</sup>, Zhiyuan Ji<sup>16</sup>, Brant Robertson<sup>17</sup>, Hannah Übler<sup>5,6</sup>, Giacomo Venturi<sup>14,D</sup>, Christopher N. A. Willmer<sup>16,D</sup>, and Joris Witstok<sup>5,6</sup>

Received 14 March 2024 / Accepted 4 November 2024

## ABSTRACT

We present the first statistical investigation of spatially resolved emission-line properties in a sample of 63 low-mass galaxies at $4 \le z < 10$ using James Webb Space Telescope (JWST)/NIRSpec Micro Shutter Assembly (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 \text{III}]\lambda5007/H\beta$ and a decreasing trend for the blended spectral complex ($H\eta + [Ne \text{III}]\lambda3869 + He I\lambda3889 + H\zeta$)/$[O \text{II}]\lambda3727$ ($3-4\sigma$ significance). These results are still valid when stacking the sample in two redshift bins (i.e. $4 \le z < 5.5$ and $5.5 \le z < 10$). The comparison with star-formation photoionisation models suggests that the ionisation parameter increases by ~0.5 dex with redshift. Under the hypothesis that radial variations in ($H\eta + [Ne \text{III}]\lambda3869 + He I\lambda3889 + H\zeta$)/$[O \text{II}]\lambda3727$ are dominated by trends in $[Ne \text{III}]\lambda3869/[O \text{II}]\lambda3727$, 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\beta$ ($7\sigma$ significance). This trend persists even after removing known active galactic nucleus candidates, and is therefore consistent with a radial gradient primarily in stellar age and secondarily in metallicity. Taken 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.

**Key words.** techniques: imaging spectroscopy – galaxies: high-redshift – galaxies: 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

characterised by gas accretion, frequent merging events, and by gas outflows due to stars and/or feedback from active galactic nuclei (AGNs; see e.g. Bischetti et al. 2019, 2021; Tripodi et al. 2024a; 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 in detail, and the presence of radial variations and chemical enrichment levels across the galaxy has been assessed. A large number of local star-forming (SF) galaxies exhibit negative

<footnote n="1">¹ University of Ljubljana FMF, Jadranska 19, 1000 Ljubljana, Slovenia</footnote>

<footnote n="2">² Dipartimento di Fisica, Università di Trieste, Sezione di Astronomia, Via G. B. Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="3">³ INAF – Osservatorio Astronomico di Trieste, Via G. Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="4">⁴ IFPU – Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34151 Trieste, Italy</footnote>

<footnote n="5">⁵ Cavendish Laboratory – Astrophysics Group, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="6">⁶ Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="7">⁷ Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="8">⁸ European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany</footnote>

<footnote n="9">⁹ Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Firenze, Italy</footnote>

<footnote n="10">¹⁰ INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="11">¹¹ Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="12">¹² Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="13">¹³ Centro de Astrobiología (CAB), CSIC-INTA, Cra. de Ajalvir Km 4, 28850 Torrejón de Ardoz, Madrid, Spain</footnote>

<footnote n="14">¹⁴ Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="15">¹⁵ Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="16">¹⁶ Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA</footnote>

<footnote n="17">¹⁷ Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote>* Corresponding author; roberta.tripodi@inaf.it</footnote>
radial gradients in their gas-phase metallicity (or simply metal-
licity), with the inner regions being found to be more chemically
enriched than 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 forma-
tion (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 galax-
ies. 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). Recently, metallicity gradients in the Large Magellanic
cloud indicated that its stellar structure is effectively outside
in (Frankel et al. 2024). 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 gra-
dients seem to evolve (see e.g. Li et al. 2022). Curti et al. (2020)
investigated the metallicity gradients of a sample of 42 galax-
ies within 1.2 < z < 2.5, finding that ~85% of these galax-
ies are characterised by metallicity gradients shallower than
0.05 dex kpc⁻<fnref n="1" />, and ~89% are consistent with a flat slope within
3σ, suggesting a mild evolution with cosmic time. These authors
also discovered three galaxies with inverted gradients, suggest-
ing recent episodes of pristine gas accretion or strong radial
flows.

The advent of the *James Webb* Space Telescope (JWST) enabled for the first time the study of metallicity gradients at *z* > 3–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 far-infrared emission lines from ALMA; e.g. Rodríguez Del Pino et al. 2024; Arribas et al. 2024; Venturi et al. 2024; Vallini et al. 2024; Wang et al. 2022). The metallicity gradients of most galaxies at redshifts of 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*, because high-*z* galaxies show irregular shapes with several (merging) components and/or clumps.

The collisionally excited [O III]λ5007 and the recombina-
tion Hβ optical emission lines are commonly used as diag-
nostics for the properties of gas in star-forming regions. In
particular, Hβ is driven primarily by ionising radiation, while
[O III]λ5007 is sensitive to the gas-phase metallicity and the
ionisation 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, ion-
isation parameter, and hardness of the radiation field in the ISM
of a star-forming galaxy (Osterbrock 1989). Indeed, a power-
ful 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<fnref n="1" /> (e.g. Nagao et al. 2006; Levesque & Richardson 2014). Neon is produced during the late evolutionary stages of massive stars and 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 ionised nebula when compared to [O II], given that the high critical density of [Ne III]λ3869 (log n[cm⁻³] = 7; see Appenzeller & Oestreicher 1988) makes its flux insensitive to the electron density even in high-density regions. The [Ne III]λ3869/[O II]λ3727 ratio has been shown to be a better diagnostic of the ionisation parameter than [O III]λ5007/[O II]λ3727, with a greater sensitivity at shorter wavelength that accommodates more of the ionising photons produced by young massive stars. It is also insensitive to reddening effects and usable as an empirical diagnostic of the ionisation parameter out to higher redshift than [O III]λ5007/[O II]λ3727, as [O III]λ5007 is emitted at longer wavelength than [Ne III]λ3869. The degeneracy between metallicity and the ionisation 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 ionisation parameter (e.g. McGaugh 1991; Kewley & Dopita 2002).

Additionally, Hβ is particularly powerful in probing ionised
regions around young, massive stars. Leaving dust extinction
aside, the intensity of the Hβ line is directly linked to the ionising
photons emitted by these stars, providing a valuable metric for
quantifying the ongoing SF activity. As a SFR diagnostic, Hβ,
like all the Balmer lines, inherits the same strengths and weak-
nesses 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, assum-
ing a star-forming origin, the equivalent width (EW) of Hβ is a
tracer of the specific SFR (sSFR) in SF galaxies, though its accu-
racy is influenced by the absorption of ionising photons by dust
within HII regions.

In this work, we aim to investigate the emission-line gradi-
ents in a sample of 63 galaxies at 4 ≤ z < 10 from the JWST
Advanced Deep Extragalactic Survey (JADES). Specifically, we
focus on spatially resolved analyses of the [O III]λ5007/Hβ
and [Ne III]λ3869/[O II]λ3727 ratios to investigate the distribution
of metallicity in high-z galaxies, and of EW_Hβ as a tracer of
the sSFR. We focus on the EW and these specific line ratios
because they are not severely affected by the fact that the PSF
is wavelength dependent, given that they are very close in wave-
length. 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 wavelengths and with different slit losses. Selecting
diagnostics based on lines that are close in wavelength also mit-
giates dust reddening effects.

This is the first attempt to perform spatially resolved metal-
licity studies at *z* > 4 in a large sample of galaxies. The paper
is structured as follows. In Sect. 2, we present the observa-
tions and the data reduction; in Sect. 3, we report our analy-
sis and results on stacked samples; in Sect. 4, we discuss the
implications of our findings, comparing our results with obser-
vations at different redshifts and with photoionisation models.
Throughout the paper, we adopt a ΛCDM cosmology from

<footnote n="1">¹ Ionisation potentials for O I, O II, and Ne II are 13.61 eV, 35.12 eV, and 40.96 eV, respectively.</footnote>
the Planck Collaboration VI (2020): $H_0 = 67.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$,
$\Omega_m = 0.315$, and $\Omega_\Lambda = 0.685$. We assume a Chabrier initial mass
function everywhere (Chabrier 2003). All physical distances are
proper distances.

**2. Observations and data reduction**

We use publicly available JWST/NIRSpec data from JADES (Eisenstein et al. 2023a; Bunker et al. 2024; 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 Great Observatories Origins Deep Survey-South (GOODS-S) extragalactic field (Giavalisco et al. 2004). Specifically, we analyse 63 galaxies at $4 \le 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. 2024; Carniani et al. 2024), 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 [$\text{O II}$]$\lambda$3727 from [Ne III]$\lambda$3869, and $H\beta$ from [$\text{O III}$]$\lambda$4959, 5007 at redshifts of $z > 4$. These data were observed using a three-shutter slit, with one-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, that is, 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 up to three pointings according to its priority (as described in Bunker et al. 2024). 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. (2024); 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 using the local background from adjacent shutters, but extended sources and sources with contaminants were pre-identified by visual inspection and self-subtraction was avoided 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 the 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–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 sam-
pling 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 III}$]$\lambda$5007, $H\beta$, [$\text{Ne III}$]$\lambda$3869, and [$\text{O II}$]$\lambda$3727. Therefore, in
this section, we briefly explain the procedure adopted for fitting
the 1D spectra along the slit of the 2D spectrum of our targets.

We separately fit [$\text{O III}$]$\lambda$5007 and $H\beta$, and [$\text{O II}$]$\lambda$3727 and
[$\text{Ne III}$]$\lambda$3869 in pairs with their corresponding underlying con-
tinua. We model the emission lines with Gaussian functions
and the two underlying continua with first-order polynomi-
als<fnref n="2" />. The doublet [$\text{O III}$]$\lambda$4959, 5007 was fitted fixing the ratio
between the peak fluxes ($\text{peak}_{[\text{O III}]\lambda4959}/\text{peak}_{[\text{O III}]\lambda5007} = 0.335$)
and the wavelength separation ($\Delta\lambda = 47.94\AA$) of the two
emission lines, and using the same FWHM for both emission
lines. Similarly, we fit the [$\text{Ne III}$]$\lambda$3869 and [$\text{Ne III}$]$\lambda$3967 dou-
blet simultaneously, adopting a ratio of 0.301 between the lat-
ter and the former, and a rest-frame wavelength separation of
98.73Å, and considering the same FWHM for both emission
lines. The doublet [$\text{O II}$]$\lambda$3727 is always blended given the res-
olution of our spectra, and so we fitted it with a single Gaus-
sian. This reduces the number of free parameters for both the
[$\text{O III}$] and the [$\text{Ne III}$]$\lambda$3869 doublets to three. Therefore, we
have eight free parameters for each pair of lines (i.e. peak
flux, peak wavelength, FWHM for [$\text{O III}$]$\lambda$5007 and the same
for $H\beta$, slope and intercept for their underlying continuum;
the same for [$\text{O II}$]$\lambda$3727 and [$\text{Ne III}$]$\lambda$3869). We explore the
eight-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 ten walkers
per parameter and 1500 trials (the typical burn-in phase is ~500
trials). Analogously, we fit the [$\text{O III}$]$\lambda$5007, $H\beta$, [$\text{O II}$]$\lambda$3727,
and [$\text{Ne III}$]$\lambda$3869 emission lines and corresponding continua in
the five-pixel boxcar extraction of the 1D spectrum<fnref n="3" />. Finally, we
compute the integrated fluxes by integrating the best-fitting func-
tions for each emission line.

Moreover, regarding the 2D spectra, we measure the EW of $H\beta$ for each trial in the chain. Then, we derive the best-fitting value from the 50th percentile of the $EW_{H\beta}$ chain, and the error on $EW_{H\beta}$ from the 16th and 84th percentiles.

**3.2. Line blending and contamination correction**

When measuring [$\text{Ne III}$]$\lambda$3869, we also include flux from $H\eta$, $H\zeta$ ($\lambda = 3890.17\AA$) and $H\i \lambda3889$, both of which are blended with [$\text{Ne III}$]$\lambda$3869 at the resolution of the prism. To estimate the contamination, we proceed as follows. For $H\zeta$ and $H\eta$, we calculate the ratio $H\delta/([\text{Ne III}]\lambda3869 + H\zeta + H\i \lambda3889 + H\eta)$ from the stacked 1D spectrum (see Sect. 3.3), find-

<footnote n="2">² As the continuum is only fitted in proximity to each pair of lines, a first-order polynomial is enough to capture the shape of the continuum, and it leads to the same results as if considering a power-law functional form.</footnote>

<footnote n="3">³ As described in Bunker et al. (2024), the five-pixel boxcar extraction is not performed on the combined 2D spectra, but on each individual exposure; the resulting 1D spectra are then combined in the final 1D spectrum.</footnote>
ing a value of 0.35 ± 0.01. This means that the ratio $H\eta/([Ne\,III]\,\lambda3869 + H\zeta + He\,I\,\lambda3889 + H\eta)$ will be less than 0.1 and the ratio $H\zeta/([Ne\,III]\,\lambda3869 + H\zeta + He\,I\,\lambda3889 + H\eta)$ will be less than 0.14, assuming the Balmer ratios from Case B recombination, $T_c = 10\,000$ K, and $n_e = 100\,\text{cm}^{-3}$ (Storey & Hummer 1995). Considering the Anomalous Balmer Emitters in JADES at high-$z$ ($z > 5$), McClymont et al. (2024) found 14 galaxies with Balmer ratios below Case B values and 12 galaxies above Case B values. Therefore, we may say that, on average, Balmer ratios in our stacked sample are aligned with Case B. We note that the contamination caused by line blending diminishes with redshift as the spectral resolution increases. We apply no reddening correction to $H\delta/([Ne\,III]\,\lambda3869$, and therefore this contamination fraction is strictly an upper limit. For $He\,I\,\lambda3889$, 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\,\lambda5877$ line and then use models to infer the $He\,I\,\lambda3889$ contamination to $[Ne\,III]\,\lambda3869$. For $n_e = 100\,\text{cm}^{-3}$ and temperatures $T_c = 5000, 10\,000$, and $20\,000$ K, Benjamin et al. (1999) give $He\,I\,\lambda3889/He\,I\,\lambda5877$ 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 $He\,I\,\lambda3889[5877]/[Ne\,III]\,\lambda3869 = 0.23 \pm 0.01$, and therefore even when assuming the highest $He\,I\,\lambda3889/He\,I\,\lambda5877$ value of 1.1, the overall contamination should be of the order of 0.25, similar to what we estimated for $H\zeta$. Thus, we have $(H\zeta + H\eta)/([Ne\,III]\,\lambda3869 + H\zeta + H\eta + He\,I\,\lambda3889) < 0.24$, and $He\,I\,\lambda3889/([Ne\,III]\,\lambda3869 + H\zeta + H\eta + He\,I\,\lambda3889) < 0.25$.

Adding these two constraints, we obtain $(H\zeta + He\,I\,\lambda3889)/([Ne\,III]\,\lambda3869 + H\zeta + H\eta + He\,I\,\lambda3889) < 0.49$, implying $[Ne\,III]\,\lambda3869/([Ne\,III]\,\lambda3869 + H\zeta + H\eta + He\,I\,\lambda3889) > 0.51$. This is an upper limit on the contamination, because, as noted above, we did not apply a reddening correction for $H\zeta/H\delta$ or $He\,I\,\lambda3889/He\,I\,\lambda5877$. In addition, it has to be noted that $He\,I$ models are still very uncertain for the optically thick $He\,I\,\lambda3889$, and are even uncertain for optically thin lines (Benjamin et al. 1999). In light of these difficulties, we do not apply the estimated correction. However, we show its magnitude in the relevant figures. To assess the validity of our estimates of the contamination, we compare to other spectrally resolved measurements of $H\eta$, $[Ne\,III]\,\lambda3869$, $He\,I\,\lambda3889$, and $H\zeta$ at high redshift. In MACS-1149-JD1 at $z = 9.11$ (Marconcini et al. 2024), in GNz-11 (Maiolino et al. 2024), and in ZS7 at $z = 7.15$ (Übler et al. 2024), $H\eta$ is a weak contaminant. In GNz-11, $H\eta$ is not even detected and $(H\zeta + He\,I\,\lambda3889)/[Ne\,III]\,\lambda3869$ is 0.4 (Maiolino et al. 2024, and at most <0.6; Bunker et al. 2024). In ZS7, $H\eta$ is also undetected and $(H\zeta + He\,I\,\lambda3889)/[Ne\,III]\,\lambda3869 = 0.38$ (Übler et al. 2024).

### 3.3. Line ratios: Stacking analysis

We aim to perform a spatially resolved study of the $[O\,III]\,\lambda5007/H\beta$ and $[Ne\,III]\,\lambda3869/[O\,II]\,\lambda3727$ line ratios in a sample of galaxies at $4 \le z < 10$. Therefore, we initially selected all the publicly available galaxies belonging to the JADES survey with S/N ≥ 5 on the flux of each of the four spectral lines<fnref n="4" /> at $4 \le z < 10$: we find that 19 objects satisfy these criteria. In Figs. 1 and 2, we show 3-arcsec false-colour cut-outs and spectra of the two galaxies with the most extended and bright

$[O\,III]\,\lambda5007, H\beta$, $[Ne\,III]\,\lambda3869$, and $[O\,II]\,\lambda3727$ emissions. We extracted the spectra from 7 pixels (one spectrum for each pixel's row from -3 to +3, with pixel scale = 0.1 arcsec/pixel<fnref n="5" />; see e.g. the 2D spectrum in Fig. 1) centred on the highest 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 an 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 and/or lower limits on the $[O\,III]\,\lambda5007/H\beta$ and/or $[Ne\,III]\,\lambda3869/[O\,II]\,\lambda3727$ line ratios with a separation of more than 1 pixel from the galaxy centre for only five sources. These show flat radial trends in their line ratios. However, the complex morphology shown in the NIRCam images makes the interpretation of these tentative gradients arduous. Details of the analysis and the 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 stacked the 2D and 1D spectra of all the sources at $4 \le z < 10$ that have S/N ≥ 5 on $[O\,III]\,\lambda5007$. This sample consists of 63 galaxies. To avoid bias towards the brightest objects in the sample, we start by dividing each spectrum by the integrated $[O\,III]\,\lambda5007$ flux (for both the 1D and 2D analyses). We then weight the normalised spectra by the inverse of the variance<fnref n="6" />. Moreover, before stacking, we aligned all the galaxies with respect to their centre. This ensures more accurate estimates of emission-line radial properties. In addition to the full stack at $4 < z < 10$, we also considered two separate subsets, at $4 \le z < 5.5$ ($z_{\text{med}} = 4.7$) and at $5.5 \le z < 10$ ($z_{\text{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 ~1.2 kpc) of the 2D full stacked spectrum. The fit to the emission lines of interest was performed as explained in Sect. 3.1. For completeness, the stacked spectra extracted from the other pixels (-2, -1, 2) are shown in Appendix C with their best-fitting models.

### 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 of 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 Table 1. While $[Ne\,III]\,\lambda3869/[O\,II]\,\lambda3727$ shows a decreasing trend towards the outskirts in all the samples, the $[O\,III]\,\lambda5007/H\beta$ ratio is almost constant with radius. The measured gradients are only tentative, with a statistical significance of between 3 and $4\sigma$; to confirm a real trend we would need

<footnote n="4">⁴ These fluxes are also publicly available, and are computed from the line fitting of the 1D spectrum of each source (Bunker et al. 2024).</footnote>

<footnote n="5">⁵ At the median redshift of the three samples studied in this work (i.e. $z = 5.5$, $z = 4.7$, $z = 6$), the pixel physical scale is 0.61 kpc, 0.66 kpc and 0.58 kpc, respectively.</footnote>

<footnote n="6">⁶ We checked that, for each galaxy spectrum, considering the variance either at a fixed wavelength or varying with wavelength did not alter our results for determining line ratios of $[O\,III]\,\lambda5007/H\beta$ and $[Ne\,III]\,\lambda3869/[O\,II]\,\lambda3727$, and $EW_{H\beta}$ on the stacked spectra. Not only the results obtained in these two ways are in agreement within the uncertainties, but also gradients and radial profiles are obtained at the same statistical significance.</footnote>
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Image and spectra of 4270. (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 grey 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. We note the complex morphology of this galaxy: the synergy between spectroscopy and imaging is vital for correctly interpreting 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, and possibly [N I]λ5198, 5200), may indicate a complex interplay between AGN photoionisation and shocks.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Image and spectra of 16 745 at z = 5.566, same as Fig. 1.</figcaption></figure>

deeper observations that allow us to reduce the uncertainties and possibly to extend the detections to larger radii. We also note that the average value of the two line ratios show a variation with redshift when comparing the two redshift-bin subsamples. The drop in [O III]λ5007/Hβ at large radii in the 5.5 ≤ z < 10 sample is not due to physical reasons but to our data reduction. In Appendix D we test an alternative data-reduction method, whereby the emission at pixel #3 is more accurate, but at the expense of lower S/N everywhere compared to our default data reduction. Our main results are unchanged even with the alternative reduction.

An important caveat of our analysis is the spectral blend between Hη, [Ne III]λ3869, He I λ3889, and Hζ in our prism data. In Sect. 3.2, we establish that Hη is low, but Hζ+He I λ3889 can make up to 0.40 of the blended spectral flux. Therefore, radial variations in (He I λ3889+Hζ)/[Ne III]λ3869 could in prin-

ciple explain our observed trends. An empirical assessment of this possibility is impossible with our data. However, we note that in the high-resolution observations of the lensed galaxy MACS-1149-JD1 (Hashimoto et al. 2018; Stiavelli et al. 2023; Bradač et al. 2024; Marconcini et al. 2024), there is a tentative negative radial gradient in [Ne III]λ3869/[O II]λ3727 (see Appendix E), with a magnitude of 0.3 dex, which is in agreement with our findings. Therefore, in the following, we study the hypothesis that the observed radial variations are dominated by true gradients in [Ne III]λ3869/[O II], which we consider the most likely hypothesis.

To better understand the implications of our results, in Fig. 5 we compare them to previous measurements at different redshifts found in the literature, and to photo-ionisation 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 ionisation parameter and/or metallicity. In particular, for the $4 \le z < 10$ stack, the radial trend appears to be due mostly to variations in the ionisation parameter. However, when we divide the sample at $z = 5.5$, we find that the two redshift bins have different average ionisation parameters, and the radial variation also shows a component consistent with positive metallicity gradients. This is a direct consequence of the different average value of the line ratios seen in Fig. 4. This topic is extensively discussed in Sect. 4.

### 3.5. Analysis of the $EW_{H\beta}$ equivalent width

Figure 6 presents the results regarding the radial profiles of $EW_{H\beta}$ for the three samples (full stacked and two z-bins). Overall, $EW_{H\beta}$ shows a decreasing trend with radius, implying that the emission of $H\beta$ 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_{H\beta}$, and the results for the slopes derived for each sample are reported in Table 1. A similar trend is also found in the radial profile of 4270 at $z = 4.023$, as shown in Fig. B.4 (image and spectra of 4270 are shown in Fig. 2). In this case, we observe a negative radial gradient, which is steepest in the central regions, followed by a flattening at the outermost radii.

This rising trend of $EW_{H\beta}$ towards the centre has some interesting implications for the nature of the innermost regions of these high-z galaxies, depending on the source of ionisation. Assuming a star-forming origin, $EW_{H\beta}$ 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_{H\beta}$ and the age at two fixed different metallicities ($\log(Z/Z_\odot) = -1.0, -0.5$). This scaling was 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_{H\beta}$ 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 ionising photons. These alternative scenarios are discussed in more detail in the following section.

## 4. Discussion

In the above sections, we present our study of the spatially resolved emission line properties of a sample of galaxies at $4 \le z < 10$, which involves stacking the 2D and 1D spectra to get the radial profiles of the line ratios. We also divided this larger sample into two redshift bins ($4 \le z < 5.5$ and $5.5 \le z < 10$), and we performed the stacking for these two subsamples. Figure 4 and Table 1 show a slightly increasing or flat (for the $4 \le z < 5.5$ subsample) trend for the $[O III]\lambda5007/H\beta$ ratio as a function of radius, while the $[Ne III]\lambda3869/[O II]\lambda3727$ ratio decreases with radius at a more than $3\sigma$ confidence level for both the larger sample and the two subsamples.

In Fig. 5, we compare our results with low-z data in the literature and with star-forming photoionisation models. In particular,

<figure><img src="image_4.png" /></figure>

**Fig. 3.** Stacked spectra using the NIRSpec/MSA prism spectra of galaxies at $4 \le 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 the $[O III]\lambda5007$ and $H\beta$ ($[O II]\lambda3727$, $[Ne III]\lambda3869$) spectral lines. Bottom panel: Spectrum extracted from the second pixel off central (i.e. ~1.2 kpc radius) of the 2D stacked spectrum. The colour coding and labels are the same as in the top panel. Even from a visual inspection of these two panels, we note a decrease in the $[Ne III]\lambda3869/[O II]\lambda3727$ ratio towards larger radii, i.e. from pixel 0 ($[Ne III]\lambda3869/[O II]\lambda3727 \sim 1$) to pixel 2 ($[Ne III]\lambda3869/[O II]\lambda3727 < 1$).

we show the line ratios measured in the local Universe, both in SF galaxies and AGNs (see the Sloan Digital Sky Survey Data Release 7, SDSS DR7, and specifically the data in pink and blue; Kauffmann et al. 2003; Brinchmann et al. 2004; Abazajian et al. 2004), in a sample of about 1500 galaxies and AGNs at $1.37 < z < 3.8$ (from the MOSFIRE Deep Evolution Field, 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 squares, 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, in prep.), and in the $z = 7.15$
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Averaged radial profiles of line ratios from the 2D stacked spectra at different redshift ranges. Cyan (violet) dots mark the logarithm of the [O III]λ5007/Hβ ([Ne III]λ3869/[O II]λ3727) ratio as a function of radius expressed in pixels (pixel 0 is the centre of the galaxy). The left panel shows the results stacking all galaxies at $4 \le z < 10$; the central panel stacking all galaxies at $4 \le z < 5.5$; the right panel stacking all galaxies at $5.5 \le z < 10$. The fit for the [O III]λ5007/Hβ ([Ne III]λ3869/[O II]λ3727) profile is shown as a dashed (dash-dotted) line with a shadowed region (3σ confidence interval). The orange vertical lines, along with the hatched region, 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, the 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 scales of 6.113 kpc/", 6.614 kpc/", and 5.832 kpc/", respectively.</figcaption></figure>

<table><caption>**Table 1.** Radial gradient of [O III]λ5007/Hβ and [Ne III]λ3869/[O II]λ3727 ratios, $EW_{H\beta}$, and metallicity for the three stacked samples.</caption><thead><tr><th>Sample</th><th>Median redshift</th><th>% of galaxies</th><th>[O III]λ5007/Hβ<br>[1/pixel]</th><th>[Ne III]λ3869/[O II]λ3727<br>[1/pixel]</th><th>EW<sub>Hβ</sub><br>[Å/pixel]</th><th>log(Z/Z<sub>☉</sub>)<br>[dex/pixel]</th></tr></thead><tbody><tr><td>4 ≤ z &lt; 10</td><td>5.5</td><td>63</td><td>0.024 ± 0.008</td><td>-0.083 ± 0.021</td><td>-23 ± 2</td><td>0.08 ± 0.05</td></tr><tr><td>4 ≤ z &lt; 5.5</td><td>4.7</td><td>30</td><td>0.013 ± 0.014</td><td>-0.085 ± 0.030</td><td>-20 ± 2</td><td>0.07 ± 0.06</td></tr><tr><td>5.5 ≤ z &lt; 10</td><td>6.0</td><td>33</td><td>0.023 ± 0.007</td><td>-0.120 ± 0.033</td><td>-35 ± 5</td><td>0.04 ± 0.05</td></tr></tbody></table>

**Notes.** The values reported in this table are the slope of the linear fit of the radial profiles shown in Fig. 4, Figs. 6, and 7. We recall that the pixel size is 0.1", and that the conversion scales (kpc/") used for each sample are 6.113 kpc/", 6.614 kpc/", and 5.832 kpc/", respectively.

host galaxy of the type-I AGN, ZS7 (yellow cross and plus Übler et al. 2024). 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, as 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 conservatively consider all the objects lying 0.2 dex above the Kewley et al. (2001) and Schawinski et al. (2007) lines as AGNs, while considering all the objects 0.2 dex below the Kauffmann et al. (2003) line as SF galaxies. Over-plotted are the photoionisation models presented in Gutkin et al. (2016), at a fixed hydrogen density of $log(n) = 2 \text{ cm}^{-3}$. Increasing or decreasing the density by a factor of ten does not alter the predicted line ratios. The variation of the ionisation parameter and the metallicity along the models is shown by the colour maps on the right.

We see two different behaviours for the full stack at $4 \le z < 10$ and for the two redshift bins. For the former, the radial trend appears to be due mostly to variations in ionisation parameter. However, the picture changes when we divide the sample at $z = 5.5$. First, we find that the two redshift bins have different average ionisation parameters ($\Delta \sim 0.5$ dex) and metallicities ($\Delta \sim 0.4$ dex), with the highest-redshift bin having the highest average value of $log U$ and the lowest value of $Z$, which is consistent with other studies (Cameron et al. 2023; Sanders et al.

2024; Curti et al. 2024; 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. We note that, in the model grid, we fixed the dust-to-metal mass ratio to $\xi_d = 0.3$. Radial variations in $\xi_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_d$ at constant metallicity of the ISM (i.e. gas and dust). However, a decreasing radial gradient of $\xi_d$ would also increase $EW_{H\beta}$, 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_d$, if such a trend were dominant, it would produce an increasing trend in $EW_{H\beta}$ with radius, which is contrary to our observations.

For all 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 ionisation parameter. For the two z-bin samples, we find tentative inverted metallicity gradients, and an almost constant ionisation parameter as a function of radius, with the only exception being the
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** [O III]λ5007/Hβ–[Ne III]λ3869/[O II]λ3727 line ratio diagram for the stacked samples. For comparison, we show z ~ 0 SDSS AGNs (galaxies) as a blue (pink) colour map with contours, z ~ 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; J1, in prep.), and the host galaxy of the z = 7.15 AGN type-I, ZS7 (yellow cross and plus symbol, depending on whether line fluxes are computed from the BLR location or the [O III] centroid, respectively; see Übler et al. 2024 for details). Contours represent the percentiles (from 15% to 95% for SDSS with a step of 10%, and 10%, 50%, and 90% for MOSDEF) of the number of objects in a sample. Spatially resolved ratios for our stacked samples of galaxies at 4 ≤ z < 10, 4 ≤ z < 5.5, and 5.5 ≤ z < 10 are plotted as solid circles in all 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 diamonds. Overlaid are the star-formation photoionisation models of Gutkin et al. (2016) at hydrogen densities of log n[cm⁻³] = 2.0. The dotted and solid coloured lines show the variation of the ionisation parameter at fixed metallicity and the variation of metallicity at fixed ionisation parameter, respectively (colour 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, Hη, 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ζ+Hη) of ≈0.3 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 de-blend [Ne III]λ3869.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Radial profile of EW(Hβ in the stacked samples. From left to right: Total stacked sample, stacked sample of z = 4–5.5 galaxies, and stacked sample of z = 5.5–10 galaxies. The fit for each sample is shown as a dashed line with a shadowed region (1σ 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. 2b). The top x-axis is as in Fig. 4.</figcaption></figure>

outermost radius of the higher-z sample, for which the ionisation
parameter drops by ~0.5 dex. Moreover, the metallicity is higher
on average (~0.4 below solar) in the z = 4–5.5 subsample than
in the higher-z subsample (~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 for example strong-line methods (e.g. Curti et al. 2024).
It is therefore unclear whether the observed increase reflects a
physical evolution, or arises from a sample-selection bias. More-
over, 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. 2024b; 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. 2023) predicts that at fixed stellar mass, metallic-
ity 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, 2024) would only reinforce this
trend.

To investigate these inverted gradients further, we also
derived the gas-phase metallicity (as probed by the oxygen
<figure><img src="image_8.png" /><figcaption>**Fig. 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 a shadowed region ($1\sigma$ confidence interval). The orange vertical lines and top x-axis are as in Fig. 4.</figcaption></figure>

abundance, O/H) by means of empirical diagnostics cali-
brated over the electron temperature ($T_e$) abundance scale,
as described in Curti et al. (2020, 2023). Specifically, we
adopted the $R_3 = \log([\text{O III}]\lambda5007/H\beta)$ and $\text{Ne}_3\text{O}_2 =
\log([\text{Ne III}]\lambda3869/[\text{O II}]\lambda3727)$ emission-line ratios, where both
diagnostics (their respective uncertainties, and the intrinsic dis-
persion of each calibration) are used simultaneously in a $\chi^2$-
minimisation procedure to break the double-branch degeneracy
of the $R_3$ indicator. 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$ subsample is, on average, ~0.3 dex
lower than in the $z = 4–5.5$ subsample. We performed a fit
to these profiles using a linear regression, and the results are
reported in Table 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 spa-
tially 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. (2024) also found small
sizes and complex kinematic features in six $5.5 < z < 7.4$ galax-
ies belonging to JADES, with three of them showing a significant
spatial velocity gradient. As a result, individual gradients can be
measured only for complex and possibly interacting systems. In
this context, understanding the synergy between photometry and
spectroscopy is crucial in order to draw the correct conclusions
regarding the nature of the gradients (see Fig. 8). In contrast,
when stacking, we ‘marginalise’ over the precision and accuracy
of the shutter alignment, making the radial gradients easier to
interpret. Assessing whether or not a sample of 63 galaxies is
sufficient to marginalise over these effects and enable us to infer
the true average gradients of the population, would require sim-
ulated observations of realistic galaxies, but this is beyond the
scope of the present paper.

As noted above, the analysis of the EW of Hβ reveals a steep
radial profile increasing towards the centre in the stacked spectra.
Additionally, while analysing the individual galaxy 4270, which
shows the clearest gradients among the five brightest galax-
ies in Appendix B, we found a similar trend in EW<sub>Hβ</sub> versus
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<sub>Hβ</sub> can be considered as a
rough proxy for sSFR. Therefore, the higher EW<sub>Hβ</sub> in the cen-
tre could be due to our sample being dominated by galaxies
undergoing centrally concentrated SF, leading to the build-up
of stellar mass in the core of galaxies (e.g. Baker et al. 2024;
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.
Shen et al. 2024; Ji & Giavalisco 2023; Ji et al. 2024a), and is
fuelled by rates of pristine or 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 oth-
erwise metal-rich gas in the central regions due to the same
gas inflows as those powering the central SF. In principle, the
presence of an AGN could also cause a central enhancement of
EW<sub>Hβ</sub>, given that AGNs tend to have high-EW nebular emission
(e.g. Caccianiga & Severgnini 2011, for EW<sub>[O III]</sub>). 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. B.3).

To ensure that the presence of AGNs does not bias our
results on the stacked samples, we excluded the 15 AGNs iden-
tified in our sample (Scholtz et al. 2023) from the stacking,
and performed the analysis on the line ratios and EW<sub>Hβ</sub> again.
The results are unchanged both for the total and the two z-bin
stacked samples. That is, even when excluding the AGNs, we
find flat radial profiles for [O III]λ5007/Hβ, and declining pro-
files for [Ne III]λ3869/[O II]λ3727, and a steep increase in EW<sub>Hβ</sub>
towards the centre. This suggests that AGNs do not explain the
observed trends in EW<sub>Hβ</sub>, which must therefore be due to stellar-
population trends.

It should also be noted that the redshift range we are probing
is relatively uncharted for this topic, and therefore our interpre-
tation is uncertain. In particular, SF galaxies in burst phases are
much more common in the young Universe (Ciesla et al. 2024;
Dressler et al. 2024; Looser et al. 2023; Endsley et al. 2024), as
also expected from models (see e.g. Faucher-Giguère 2018;
Tacchella et al. 2020). These objects have much higher EW<sub>Hβ</sub>
and EW<sub>[O III]</sub> than star-forming galaxies in the local Universe
(e.g. Smit et al. 2014; van der Wel et al. 2011; Maseda et al.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Three-arcsecond false-colour cut-outs 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; this 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.</figcaption></figure>

2014; Boyett et al. 2024), and therefore they may equal or even exceed the EWs of AGNs (for type-2 AGN, $EW_{[OIII]} < 500$ Å; e.g. Caccianiga & Severgnini 2011). In principle, the observed trends could be caused by the fact that, on average, $EW_{H\beta}$ 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. 2024b). However, we still observe the same trend between $EW_{H\beta}$ and radius even after splitting our sample into two redshift bins, which suggests that the trend is real and not due to the observational conspiracy between $EW_{H\beta}$, size, and redshift. Alternatively, if starburst galaxies were to have both higher $EW_{H\beta}$ and smaller size than non-starburst galaxies, this could also produce the decreasing radial gradients of $EW_{H\beta}$ in the stack. This is because, under this hypothesis, starburst galaxies (with their high $EW_{H\beta}$) would dominate the central pixels of the stack, while more extended, less bursty galaxies (with lower $EW_{H\beta}$) 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_{H\beta}$) requires further division of 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- versus wide-band imaging (to measure EW gradients).

Alternatively, IFU observations with NIRSpec can be used to directly measure spectral gradients in galaxies (Rodríguez Del Pino et al. 2024; Arribas et al. 2024; Venturi et al. 2024). However, NIRSpec/IFS is best suited for large, bright galaxies, primarily because of its low multiplexity. Using this instrument configuration to target a large sample of low-mass, high-z galaxies with a depth comparable to that of 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. 2024). For all these reasons, we do not believe NIRSpec/IFS to be the most suitable tool for building large samples – at least in the redshift and mass range explored in this article.
**5. Summary and conclusions**

We studied the radial profiles of the [O III]λ5007/Hβ and
[Ne III]λ3869/[O II]λ3726, 3729 line ratios and of EW<sub>Hβ</sub> in a
sample of 63 galaxies at 4 ≤ z < 10 using publicly available
JWST/NIRSpec data from JADES. We focused on deep, spa-
tially resolved spectroscopy in the GOODS-S extragalactic field,
as our aim is to study spatially resolved emission line gradients.
Our findings can be summarised 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 in the line ratios. Specifically, for the two galaxies 4270 at z = 4.023 and 16745 at z = 5.5666, these trends are detected or are marginally detected up to ~2 kpc from the centre of the galaxy. These could 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 of the lines of interest and to marginalise over the complex galaxy morphologies, we stacked the whole sample of 63 galaxies and also two subsamples, which were created by dividing the larger sample into two redshift bins (4 ≤ z < 5.5, and 5.5 ≤ z < 10). We find a slightly increasing or flat radial trend for [O III]λ5007/Hβ, and a decreasing radial trend for [Ne III]λ3869/[O II]λ3726, 3729 in all the three stacked samples.

- By comparing our findings with photoionisation models, we find that the two redshift bins have different average ionisation parameters (Δ ~ 0.5 dex, increasing with redshift), with tentative evidence of an inverted metallicity gradient. The metallicity is on average ~0.3 dex higher in the z = 4–5.5 subsample, implying a stronger evolution than measured using strong-line methods, for example. This is still valid if adopting the calibrations presented in Curti et al. (2020) to convert line ratios to metallicity. It remains difficult to decipher 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<sub>Hβ</sub>, 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 gradi-
ents and of a negative gradient of EW<sub>Hβ</sub> in a stacked sample
of 4 ≤ z < 10 galaxies. There are several possible avenues
to explore in order to confirm and refine this study. For galax-
ies 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.
2023b). 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 tar-
gets, 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 sub-
traction (to avoid self subtraction) or to the observations (e.g.
the NIRSpec/IFS; Böker et al. 2022; Rodríguez Del Pino et al.
2024; Arribas et al. 2024).

**Acknowledgements.** 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 2022), Matplotlib (Hunter 2007), SciPy (Gommers et al. 2023), emcee (Foreman-Mackey et al. 2013).

**References**

Abazajian, K., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2004, AJ, 128, 502

Alves de Oliveira, C., Birkmann, S. M., Böker, T., et al. 2018, SPIE Conf. Ser., 10704, 107040Q

Appenzeller, I., & Oestreicher, R. 1988, AJ, 95, 45

Arribas, S., Perna, M., Rodríguez Del Pino, B., et al. 2024, A&A, 688, A146

Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

Baker, W. M., Maiolino, R., Belfiore, F., et al. 2023, MNRAS, 519, 1149

Baker, W. M., Tacchella, S., Johnson, B. D., et al. 2024, Nat. Astron. [arXiv:2306.02472]

Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5

Belfiore, F., Maiolino, R., Tremonti, C., et al. 2017, MNRAS, 469, 151

Benjamin, R. A., Skillman, E. D., & Smits, D. P. 1999, ApJ, 514, 307

Berg, D. A., Skillman, E. D., Croxall, K. V., et al. 2015, ApJ, 806, 16

Bischetti, M., Maiolino, R., Carniani, S., et al. 2019, A&A, 630, A59

Bischetti, M., Feruglio, C., Piconcelli, E., et al. 2021, A&A, 645, A33

Böker, T., Arribas, S., Lützgendorf, N., et al. 2022, A&A, 661, A82

Bonaventura, N., Jakobsen, P., Ferruit, P., Arribas, S., & Giardino, G. 2023, A&A, 672, A40

Boettke, K., Bunker, A. J., Curtis-Lake, E., et al. 2024, MNRAS, 535, 1796

Bradač, M., Strait, V., Mowla, L., et al. 2024, ApJ, 961, L21

Bresolin, F., Kennicutt, R. C., & Ryan-Weber, E. 2012, ApJ, 750, 122

Bresolin, F., Kudritzki, R.-P., Urbaneja, M. A., et al. 2016, ApJ, 830, 64

Brinchmann, J. 2023, MNRAS, 525, 2087

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2024, A&A, 690, A288

Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44

Caccianiga, A., & Severgnini, P. 2011, MNRAS, 415, 1928

Calabrò, A., Pentericci, L., Santini, P., et al. 2024, A&A, 690, A290

Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023, A&A, 677, A115

Carniani, S., Venturi, G., Parlanti, E., et al. 2024, A&A, 685, A99

Chabrier, G. 2003, PASP, 115, 763

Choi, E., Brennan, R., Somerville, R. S., et al. 2020, ApJ, 904, 8

Ciesla, L., Elbaz, D., Ilbert, O., et al. 2024, A&A, 686, A128

Cresci, G., Mannucci, F., & Curti, M. 2019, A&A, 627, A42

Curti, M., Maiolino, R., Cirasuolo, M., et al. 2020, MNRAS, 492, 821

Curti, M., Hayden-Pawson, C., Maiolino, R., et al. 2022, MNRAS, 512, 4136

Curti, M., D'Eugenio, F., Carniani, S., et al. 2023, MNRAS, 518, 425
Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2024, A&A, 684, A75

Davé, R., Finlator, K., & Oppenheimer, B. D. 2011, MNRAS, 416, 1354

Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1

de Graaff, A., Rix, H.-W., Carniani, S., et al. 2024, A&A, 684, A87

Dekel, A., Sari, R., & Ceverino, D. 2009, ApJ, 703, 785

D'Eugenio, F., Pérez-González, P. G., Maiolino, R., et al. 2024, Nat. Astron., 8, 1443

Dressler, A., Rieke, M., Eisenstein, D., et al. 2024, ApJ, 964, 150

Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023a, ArXiv e-prints [arXiv:2306.02465]

Eisenstein, D. J., Johnson, B. D., Robertson, B., et al. 2023b, ArXiv e-prints [arXiv:2310.12340]

Endsley, R., Stark, D. P., Whitler, L., et al. 2024, MNRAS, 533, 1111

Faucher-Giguère, C.-A. 2018, MNRAS, 473, 3717

Ferruit, P., Jakobsen, P., Giardino, G., et al. 2022, A&A, 661, A81

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Frankel, N., Andrae, R., Rix, H. W., Povick, J., & Chandra, V. 2024, ArXiv e-prints [arXiv:2403.08516]

Fumagalli, M., Patel, S. G., Franx, M., et al. 2012, ApJ, 757, L22

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091

Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93

Gibson, B. K., Courty, S., Cunnama, D., & Mollá, M. 2013, Asociacion Argentina de Astronomia La Plata Argentina Book Series, 4, 57

Gommers, R., Virtanen, P., Burovski, E., et al. 2023, https://doi.org/10.5281/zenodo.7655153

Gutkin, J., Charlot, S., & Bruzual, G. 2016, MNRAS, 462, 1757

Hashimoto, T., Laporte, N., Mawatari, K., et al. 2018, Nature, 557, 392

Henry, R. B. C., & Worthey, G. 1999, PASP, 111, 919

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80

Ji, Z., & Giavalisco, M. 2023, ApJ, 943, 54

Ji, Z., Williams, C. C., Tacchella, S., et al. 2024a, ApJ, 974, 135

Ji, Z., Williams, C. C., Suess, K. A., et al. 2024b, ArXiv e-prints [arXiv:2401.00934]

Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22

Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 33

Kewley, L. J., & Dopita, M. A. 2002, ApJS, 142, 35

Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121

Kewley, L. J., Rupke, D., Zahid, H. J., Geller, M. J., & Barton, E. J. 2010, ApJ, 721, L48

Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15

Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716

Lara-López, M. A., Bongiovanni, A., Cepa, J., et al. 2010, A&A, 519, A31

Levesque, E. M., & Richardson, M. L. A. 2014, ApJ, 780, 100

Li, H., Mao, S., Cappellari, M., et al. 2018, MNRAS, 476, 1765

Li, Z., Wang, X., Cai, Z., et al. 2022, ApJ, 929, L8

Looser, T. J., D'Eugenio, F., Maiolino, R., et al. 2023, ArXiv e-prints [arXiv:2306.02470]

Magrini, L., Stanghellini, L., Corbelli, E., Galli, D., & Villaver, E. 2010, A&A, 512, A63

Maiolino, R., Scholtz, J., Witstok, J., et al. 2024, Nature, 627, 59

Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115

Marconcini, C., D'Eugenio, F., Maiolino, R., et al. 2024, MNRAS, 533, 2488

Maseda, M. V., van der Wel, A., Rix, H.-W., et al. 2014, ApJ, 791, 17

McClymont, W., Tacchella, S., D'Eugenio, F., et al. 2024, ArXiv e-prints [arXiv:2405.15859]

McGaugh, S. S. 1991, ApJ, 380, 140

Moustakas, J., Kennicutt, R. C., & Tremonti, C. A. 2006, ApJ, 642, 775

Nagao, T., Maiolino, R., & Marconi, A. 2006, A&A, 459, 85

Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, ApJS, 269, 33

Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141

Nelson, E. J., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2016, ApJ, 828, 27

Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, MNRAS, 525, 2864

Oppenheimer, B. D., & Davé, R. 2008, MNRAS, 387, 577

Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Mill Valley: University Science Books)

Pezzulli, G., Fraternali, F., & Binney, J. 2017, MNRAS, 467, 311

Planck Collaboration VI. 2020, A&A, 641, A6

Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJ, 936, L14

Portinari, L., & Chiosi, C. 1999, A&A, 350, 827

Reddy, N. A., Kriek, M., Shapley, A. E., et al. 2015, ApJ, 806, 259

Rieke, M., & Jades Collaboration 2023, Am. Astron. Soc. Meet. Abstr., 55, 212.01

Robbins, R. R. 1968, ApJ, 151, 511

Rodríguez Del Pino, B., Perna, M., Arribas, S., et al. 2024, A&A, 684, A187

Rupke, D. S. N., Kewley, L. J., & Barnes, J. E. 2010a, ApJ, 710, L156

Rupke, D. S. N., Kewley, L. J., & Chien, L. H. 2010b, ApJ, 723, 1255

Samland, M., Hensler, G., & Theis, C. 1997, ApJ, 476, 544

Sanders, R. L., Shapley, A. E., Topping, M. W., Reddy, N. A., & Brammer, G. B. 2024, ApJ, 962, 24

Schawinski, K., Thomas, D., Sarzi, M., et al. 2007, MNRAS, 382, 1415

Scholtz, J., Maiolino, R., D'Eugenio, F., et al. 2023, ArXiv e-prints [arXiv:2311.18731]

Scoville, N., Faisst, A., Weaver, J., et al. 2023, ApJ, 943, 82

Shao, Y., Wang, R., Carilli, C. L., et al. 2019, ApJ, 876, 99

Shen, L., Papovich, C., Matharu, J., et al. 2024, ApJ, 963, L49

Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58

Sobral, D., Best, P. N., Smail, I., et al. 2014, MNRAS, 437, 3516

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Stiavelli, M., Morishita, T., Chiaberge, M., et al. 2023, ApJ, 957, L18

Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41

Suzuki, T. L., Minowa, Y., Koyama, Y., et al. 2019, PASJ, 71, 69

Sánchez Almeida, J., Caon, N., Muñoz-Tuñón, C., Filho, M., & Cerviño, M. 2018, MNRAS, 476, 4765

Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Science, 348, 314

Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 458, 242

Tacchella, S., Carollo, C. M., Förster Schreiber, N. M., et al. 2018, ApJ, 859, 56

Tacchella, S., Forbes, J. C., & Caplar, N. 2020, MNRAS, 497, 698

Thielemann, F.-K., Kratz, K.-L., Pfeiffer, B., et al. 1994, Nucl. Phys. A, 570, 329

Tripodi, R., Lelli, F., Feruglio, C., et al. 2023, A&A, 671, A44

Tripodi, R., Scholtz, J., Maiolino, R., et al. 2024a, A&A, 682, A54

Tripodi, R., Feruglio, C., Fiore, F., et al. 2024b, A&A, 689, A220

Troncoso, P., Maiolino, R., Sommariva, V., et al. 2014, A&A, 563, A58

Trussler, J. A. A., Adams, N. J., Conselice, C. J., et al. 2023, MNRAS, 523, 3423

Übler, H., Naab, T., Oser, L., et al. 2014, MNRAS, 443, 2092

Übler, H., Maiolino, R., Pérez-González, P. G., et al. 2024, MNRAS, 531, 355

Vallini, L., Witstok, J., Sommovigo, L., et al. 2024, MNRAS, 527, 10

van der Wel, A., Straughn, A. N., Rix, H. W., et al. 2011, ApJ, 742, 111

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28

Venturi, G., Carniani, S., Parlanti, E., et al. 2024, A&A, 691, A19

Wang, X., Jones, T. A., Treu, T., et al. 2019, ApJ, 882, 94

Wang, X., Jones, T. A., Treu, T., et al. 2020, ApJ, 900, 183

Wang, X., Jones, T., Vulcani, B., et al. 2022, ApJ, 938, L16

Williams, C. C., Tacchella, S., Maseda, M. V., et al. 2023, ApJS, 268, 64

Zhang, H.-X., Hunter, D. A., Elmegreen, B. G., Gao, Y., & Schruba, A. 2012, AJ, 143, 47

Zolotov, A., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 450, 2327
# 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).

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. 2024). 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 σ of the source Gaussian, for which we quote the 50th percentile as fiducial value, and half the 16th–84th inter-percentile range as uncertainty. An example galaxy profile is reported in Fig. A.1. The σ of each galaxy in our sample is shown as a dot in Fig A.2.

<figure><img src="image_10.png" /><figcaption>**Fig. A.1.** 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 σ (bottom right).</figcaption></figure>

For the MSA PSF, we use the model from webbpsf, for which we report the Gaussian-equivalent σ, defined by σ ≡

FWHM/√ln 256 (dashed line in Fig. A.2). We tested that following the approach of de Graaff et al. (2024) gives consistent results.

<figure><img src="image_11.png" /><figcaption>**Fig. A.2.** Comparison of the radial extent of our data (using [O III]λ5007, circles) to the PSF model from webbpsf (dashed line) and to the model from the data reduction pipeline (dotted line). The galaxies in our sample extend beyond the PSF σ.</figcaption></figure>

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 σ of webbpsf, but for each realisation, we sum 20 exposures with an added random offset of ±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. A.2), 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 σ = 0.03 arcsec.

# Appendix B: Analysis of individual sources

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 \le z < 10$, we initially select all the publicly available galaxies belonging to the JADES survey with $S/N>\approx 5$ on the flux of each spectral line[^7] at $4 \le z < 10$: we find that 19 objects satisfy these criteria. Figure 1b 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 ~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). In Fig. B.1, we show all the 1D spectra extracted with their best-fitting models, and we report the fluxes in Table B.1.

<footnote>[^7]: These fluxes are also publicly available, and are computed from the line fitting of the 1D spectrum of each source (Bunker et al. 2024).</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. B.1.** 1D spectra extracted from 7 pixels (from pixel -3 to +3, see Fig. 1b) of the 2D spectrum of 4270 at z = 4.023. The cyan (violet) line is the best-fit function for the continuum and [O III]λ5007, Hβ ([O II]λ3727, [Ne III]λ3869) spectral lines.</figcaption></figure>

<table><caption>**Table B.1.** Fluxes derived from the spectra extracted from 7 pixels of the 2D spectrum of 4270 (see Fig. B.1).</caption><thead><tr><th>Line</th><th>Pix#-3</th><th>Pix#-2</th><th>Pix#-1</th><th>Pix#0</th><th>Pix#1</th><th>Pix#2</th><th>Pix#3</th></tr><tr><td></td><td colspan="7">1e-20 erg s<sup>-1</sup> cm<sup>-2</sup></td></tr></thead><tbody><tr><td>[O III]</td><td>3.2<sup>+0.8</sup><sub>-0.7</sub></td><td>18.6 ± 0.9</td><td>66.4 ± 0.96</td><td>104.6 ± 1.0</td><td>69.6 ± 0.9</td><td>41.1<sup>+1.5</sup><sub>-1.0</sub></td><td>40.6 ± 1.2</td></tr><tr><td>Hβ</td><td>1.6<sup>+1.0</sup><sub>-0.9</sub></td><td>5.5<sup>+1.5</sup><sub>-1.4</sub></td><td>12.4 ± 1.2</td><td>18.8 ± 1.2</td><td>11.5 ± 0.9</td><td>4.7<sup>+2.2</sup><sub>-0.9</sub></td><td>4.5<sup>+1.6</sup><sub>-1.5</sub></td></tr><tr><td>[Ne III]λ3869</td><td>—</td><td>&lt; 1.8</td><td>4.76<sup>+1.8</sup><sub>-1.6</sub></td><td>7.7 ± 1.8</td><td>5.8 ± 1.5</td><td>2.8<sup>+1.0</sup><sub>-0.9</sub></td><td>3.3<sup>+1.9</sup><sub>-1.7</sub></td></tr><tr><td>[O II]λ3727</td><td>—</td><td>6.1<sup>+2.0</sup><sub>-1.8</sub></td><td>10.4<sup>+2.0</sup><sub>-1.9</sub></td><td>14.8 ± 2.1</td><td>13.2 ± 1.85</td><td>6.11<sup>+1.4</sup><sub>-1.3</sub></td><td>5.1<sup>+2.5</sup><sub>-2.1</sub></td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Fig. B.2.** Five best cases of radial profiles of line ratios. Cyan (violet) dots mark the logarithm of [O III]λ5007/Hβ ([Ne III]λ3869/[O II]λ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. 2b). The ID of each object and the redshift are specified in the top left corner of each plot.</figcaption></figure>
<figure><img src="image_14.png" /></figure>

**Fig. B.3.** [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 ~ 0 SDSS AGNs (galaxies) as blue (pink) colormap with contours, z ~ 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. 2024 for details). Spatially resolved ratios for our stacked samples of galaxies at 4 ≤ z < 10, at 4 ≤ z < 5.5, and at 5.5 ≤ 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 photoionisation models of Gutkin et al. (2016) at hydrogen densities log n[cm⁻³] = 2.0. The dotted and solid coloured lines show the variation of the ionisation parameter at fixed metallicity and the variation of metallicity at fixed ionisation 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, Hη 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ζ+Hη) of ≈0.3 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.

Another interesting object is 16745, which also has a very bright [O III]λ5007 emission, extended over 7 pixels that corresponds to ~ 4 kpc at the galaxy redshift (see Figs. 2b and 2c for 2D and 1D spectra). de Graaff et al. (2024) found disc-like kinematics in 16745, with v/σ ~ 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.

<figure><img src="image_15.png" /></figure>

**Fig. B.4.** 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. 2b).

Top left and top right panels of Fig. B.2 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. 1a and 2a), 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. B.2, 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-centre pixels, especially for the [Ne III]λ3869/[O II]λ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. ~ 2 kpc from the galaxy centre). Figure B.3 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 photoionisation 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 ~ 2 and z ~ 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 ionisation 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.

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. 1a and 2a) 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. B.3, 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 photoionisation models, at fixed hydrogen density (log n[cm⁻³] = 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 ionisation parameter, with higher ionisation towards larger radii, although the measurement uncertainties are large. Indeed, since the major photoionisation sources (e.g. stars, AGN) are typically concentrated in the centre, this scenario of higher ionisation in the outskirts appears unlikely, even if considering the presence of a companion. Therefore, we conservatively conclude that the ionisation parameter radial profile is consistent with flat, given the large errorbars along the x-axis. For this target, we also find that EW<sub>Hβ</sub> decreases with radius (Fig. B.4), which would be in agreement with both the AGN scenario as well as a central starburst, fuelled by inflow of relatively metal-poor gas.

# Appendix C: Spectra of the stacked sample

In Fig. C.1 we show the stacked spectra extracted from pixels -2,-1,1,3 with their best-fitting models.

<figure><img src="image_16.png" /></figure>

**Fig. C.1.** Stacked spectra using the NIRSpec/MSA prism spectra of galaxies at $4 \le z < 10$ belonging to the JADES survey. Left panel: spectrum extracted from two pixel off central downwards (i.e. ~ 1.2 kpc radius) of the 2D stacked spectrum. The cyan (violet) line is the best-fit function for the continuum and [O III]λ5007, Hβ ([O II]λ3727, [Ne III]λ3869) spectral lines. Central panel: spectrum extracted from one pixel off central downwards of the 2D stacked spectrum. Colour coding and labels are the same as left panel. Right panel: spectrum extracted from one pixel off central upwards of the 2D stacked spectrum. Colour coding and labels are the same as left panel.

# Appendix D: Alternative data reduction

Our background subtraction strategy based on nodding ensures the highest accuracy, but could introduce spurious trends due to self subtraction of faint, extended features. To test the possible impact of our default background subtraction strategy, in this section we repeat our analysis starting from an alternative data reduction. Unlike the default procedure, this method uses only pairs of nod observations separated by two microshutters. This data reduction reduces the overall signal-to-noise of the data. Overall, the recovered radial trends (see Fig. D.1) are fully consistent with the findings presented in Fig. 4. The only exception is the [O III]5007/Hβ ratio at pixel #3 in the $5.5 \le z < 10$ sample: in this new reduction it is in agreement with a flat trend. Therefore, we conclude that the drop of pix #3 in the right panel of in Fig. 4 is not due to physical reasons but just to our original data reduction, which increases overall the S/N but part of the extended emission can be missed. In any case, this does not affect our main conclusions.
<figure><img src="image_17.png" /></figure>

**Fig. D.1.** 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β ([Ne III]λ3869/[O II]λ3727) ratio as a function of radius expressed in pixels (pixel 0 is the centre of the galaxy). The left panel shows the results stacking all galaxies at $4 \le z < 10$; the central panel stacking all galaxies at $4 \le z < 5.5$; the right panel stacking all galaxies at $5.5 \le z < 10$.

## Appendix E: [Ne III]λ3869/[O II] and (He I λ3889+Hζ)/[O II] maps of MACS-1149-JD1

Marconcini et al. (2024) present JWST/NIRSpec integral-field spectroscopy observations of the lensed galaxy MACS1149-JD1, as part of the GA-NIFS programme. This galaxy, at $z = 9.113$ and with magnification $\mu \approx 10$ (Hashimoto et al. 2018; Stiavelli et al. 2023; Bradač et al. 2024) is significant in being one of very few well resolved sources at $z > 9$, and with available integral-field spectroscopy data (IFS). This system has been characterised as a low-metallicity starburst (Stiavelli et al. 2023), with a rotating-disc component (Bradač et al. 2024; although see Marconcini et al. 2024 for a different interpretation). NIRSpec/IFS enables to measure spatially resolved emission-line ratios and star-formation rates, and Marconcini et al. (2024) find $3-\sigma$ evidence of a spatially resolved anticorrelation of the gas-phase metallicity and the SFR density in MACS1149-JD1.

In Fig. E.1, we present new emission-line ratio maps of [Ne III]λ3869/[O II] and (He I λ3889+Hζ)/[O II] for MACS-1149-JD1, which were not shown in Marconcini et al. (2024). These ratios were measured using Gaussian fitting (Marconcini et al. 2024), and the maps are clipped below S/N=3. Given the high resolution of the JWST/NIRSpec G395H spectrum, we resolve spectrally [Ne III]λλ3869, 3967 from the blend of He I λ3889+Hζ. We find an integrated (He I λ3889+Hζ)/[O II] = 0.5, with too low S/N for measuring a clear radial trend. In contrast, there is tentative evidence of a variation in [Ne III]λ3869/[O II], with the centre enhanced with respect to the outskirts, similar to what we find in our stacked analysis. The finding of a decreasing radial trend in [Ne III]λ3869/[O II] lends credence to our interpretation of the spectrally unresolved gradient in ([Ne III]λ3869+He I λ3889+Hζ)/[O II] as implying a gradient in [Ne III]λ3869/[O II].

<figure><img src="image_18.png" /></figure>

**Fig. E.1.** [Ne III]λ3869/[O II] (left) and (He I λ3889+Hζ)/[O II] (right) maps of MACS-1149-JD1. Grey solid lines are arbitrary [O III] levels. Spaxels with S/N≤ 3 are masked. Y-axis is the same for both panels.