Chemical abundance patterns at z=1-7: Carbon, Nitrogen, Neon, Sulphur, and Argon abundances from JWST spectroscopy

Introduction

The build-up of heavy elements provides one of the most direct observational windows into the processes that regulate galaxy evolution. Metals trace the cumulative effects of star formation (SF), gas accretion, and stellar feedback, while simultaneously shaping the physical properties of the interstellar medium (ISM) by regulating cooling, SF efficiency, dust content and opacity (e.g., [159, 104, 21, 50, 108, 130]). Following the evolution of chemical abundances across cosmic time is therefore key to reconstructing the baryon cycle in galaxies.

Given that oxygen is the most abundant metal in mass, its overall abundance is often used as a tracer of the integrated enrichment history. Relative abundance ratios of different elements, instead, offer a differential probe of the underlying nucleosynthetic channels. Different elements are, in fact, produced by stars of different masses and on different timescales, making their ratios sensitive to the star formation history (SFH), the initial mass function (IMF), and the occurrence of inflows and outflows (e.g., [108, 38]). In massive stars ($\gtrsim8$ $M_{\odot}$) oxygen is efficiently released into the ISM through stellar winds and core-collapse supernovae (CCSNe) on short timescales of $\approx3$–$40$ Myr (e.g., [178, 35, 105, 100]). Under typical H ii region conditions, the total oxygen abundance is well approximated by the sum of O$^{+}$ and O$^{2+}$, with higher ionisation stages usually providing negligible corrections (e.g., [155]). The emission lines associated with these ionic species lie in the rest-frame optical, providing a rich set of diagnostics for oxygen abundances; O/H is therefore commonly adopted as a proxy for the gas-phase metallicity of the ISM ([95]).

The nitrogen-to-oxygen ratio (N/O) is one of the most extensively studied abundance ratios in the ISM of star-forming galaxies, as the emission lines of nitrogen are also readily accessible in the same optical spectral range as the prominent oxygen lines. Nitrogen can be produced through two main channels, referred to as primary and secondary production. In the primary channel, nitrogen originates from carbon and oxygen freshly synthesized within the star itself: these newly formed nuclei are mixed into the hydrogen-burning regions and converted into nitrogen through the CNO cycle. As a result, the nitrogen yield in this case is largely independent of the star’s initial metallicity. Primary nitrogen can be produced both by massive stars on short timescales—tens of Myr—where fast rotational and convective mixing processes are highly efficient (e.g., [74, 118, 34]), and by low- and intermediate-mass stars (LIMS hereafter; $M \lesssim8\,M_{\odot}$), during their asymptotic giant branch (AGB) phase (e.g., [174]). The latter happens on much longer enrichment timescales—hundreds of Myr to several Gyr.

In contrast, the secondary channel for N production relies on C and O already present in the progenitor gas cloud. In this case, the CNO cycle transforms pre-existing heavy elements into nitrogen, leading to a metallicity-dependent yield (e.g., Henry et al. [74]; Vincenzo et al. [174]; Maiolino & Mannucci [108]). This contribution arises mainly from LIMS in their AGB phase (Kobayashi et al. [101]; Nomoto et al. [125]; Ventura et al. [172]). Both channels thus operate over a wide range of stellar masses and timescales (e.g., Vincenzo et al. [174]; Maiolino & Mannucci [108]). The N/O ratio is therefore a sensitive clock of chemical enrichment and has attracted considerable attention in recent JWST studies of high-$z$ galaxies (e.g., Bunker et al. [24]; Ji et al. 2024; Rogers et al. [138]; Arellano-Córdova et al. [5]; Curti et al. [41]; Hayes et al. [73]; Napolitano et al. [122]; Navarro-Carrera et al. [123]; Stiavelli et al. [157]; Topping et al. [167]; Berg et al. [13]; Cataldi et al. [32]; Cameron et al. [26]; Morel et al. [119]; Ji et al. 2026; Schaerer et al. [145]; Zhang et al. [180]).

Carbon plays a crucial role in the chemical evolution of galaxies, along with oxygen and nitrogen. It is produced in the interior of stars through the triple-$\alpha$ process during helium burning, a process that is independent of the star’s initial metallicity. Carbon can be released into the ISM through two main pathways, with by different delay times. LIMS produce carbon during He burning through the triple-$\alpha$ process and during thermal pulses in the AGB phase. The newly produced C is transported to the stellar envelope through third dredge-up episodes and subsequently released into the ISM through stellar winds. The carbon yield generally peaks for stellar masses around $2$–$4\,M_\odot$ and decreases towards higher masses because hot-bottom burning (HBB) converts C into N (Kobayashi et al. [101]; Karakas & Lattanzio [89]). Massive stars ($M > 8\,M_\odot$), on the other hand, release carbon on shorter timescales (a few Myr), primarily through CCSNe explosion or possibly through the strong stellar winds of the Wolf-Rayet (WR) phase (e.g., Woosley et al. [178]; Meynet & Maeder [118]).

Finally, Ne, S, and Ar are all $\alpha$-elements: their most abundant isotopes have mass numbers that are multiples of four, and they are synthesised through $\alpha$-capture reactions in massive stars on the same timescales (up to few tens Myr) as oxygen (e.g., Izotov et al. 2006; Croxall et al. [37]; Esteban et al. [52]; Rogers et al. [136]; Arellano-Córdova et al. [6]). If this was the only nucleosynthetic channel, their abundance ratios would remain approximately constant with metallicity. However, nucleosynthesis models also predict that Type Ia supernovae (SNe Ia) contribute non-negligibly to the S and Ar budget (e.g., Iwamoto et al. [78]; Johnson [86]; Kobayashi et al. [100]). This contribution may manifest itself as an observable trend of S/O and Ar/O with metallicity. The nucleosynthetic time delay of $\sim100$ Myr – Gyr between CCSNe and SNe Ia opens up the possibility of using O/Ar (and, possibly, O/S) as analogues of the classical [$\alpha$/Fe] diagnostic – a well-established tracer of chemical abundances in resolved stellar populations (e.g., Tinsley [164]; Matteucci & Brocato [114]; Tolstoy et al. [165]; Matteucci [113]). In this framework, the ISM abundance ratio between an $\alpha$-element produced predominantly in CCSNe (such as oxygen) and one with a non-negligible SNe Ia contribution (such as argon or sulfur) carries an imprint of the star formation history of the host system.

In the local Universe, most observational studies find that Ne/O, S/O, and Ar/O are broadly constant with metallicity, within the scatter and measurement uncertainties (e.g., Kennicutt et al. [93]; Izotov et al. 2006; Croxall et al. [37]; Arellano-Córdova et al. [3]; Berg et al. [14]; Berg et al. 2022; Rogers et al. [136]; Arellano-Córdova et al. [4]; Esteban et al. [54]), while C/O rises with oxygen abundance (e.g., Berg et al. [16]; Berg et al. [12]). However, several works report mild trends with metallicity for S/O and Ar/O for specific samples (e.g., Dors et al. [49]; Díaz & Zamora [48]; Arnaboldi et al. [8]; Brazzini et al. [23]; Scholte et al. 2026), suggesting that the constancy of these abundance ratios with metallicity is not firmly established.

At high-$z$, the situation is even less constrained: the nebular emission lines used for tracing some of these elements (e.g. [Ar iii] $\lambda$7136, [Ne iii] $\lambda$3869, [S ii] $\lambda\lambda$6716,6730, and [S iii] $\lambda$9069,9532), have been detected in distant galaxies only thanks to the sensitivity of the Near InfraRed Spectrograph (NIRSpec) on board JWST. Early results, based on samples of only a dozen galaxies, suggest that the Ne/O abundance ratio remains broadly consistent with the local Universe, showing little or no evidence for redshift evolution (Arellano-Córdova et al. [3, 6]; Isobe et al. [76]; Stanton et al. [154]). The situation is less clear for Ar and S. Bhattacharya et al. [18] first used the O/Ar–Ar/H plane to investigate the relative contributions of CCSNe and SNe Ia to the chemical enrichment of high-$z$ star-forming galaxies. Later studies reported predominantly subsolar Ar/O ratios, consistent with a limited contribution from SNe Ia relative to CCSNe (Stanton et al. [154]; Rogers et al. [139]; Foley et al. 2026). Rogers et al. [139] also found subsolar S/O, which they interpreted within the same framework of predominantly CCSN-driven enrichment. Overall, existing studies are based on small samples, highlighting the need for larger statistics to robustly establish the evolution of abundance ratios across cosmic time.

In this work, we exploit the capabilities of JWST to assemble a spectroscopic sample of 244 star-forming galaxies at $z \sim1$–$7$ with direct metallicity measurements based on the temperature-sensitive [O iii] $\lambda$4363 line. Within this parent sample, we identify subsamples with detections of emission lines tracing C, Ne, N, S, or Ar. Building on the observational framework established in Cataldi et al. [32] (C26 hereafter) for nitrogen, we extend the analysis to the determination of direct-method abundances for Ar, C, S and Ne. We further present calibrations for N/O, S/O, Ar/O and Ne/O based on strong line ratios. With this large sample of abundance measurements at high redshift, we aim to shed new light on the behaviour of abundance trends as a function of metallicity.

The paper is structured as follows. § 2 describes the observational dataset. § 3 presents the abundance measurement methods. § 4 derives and discusses the strong-line calibrations, including an investigation of their secondary dependencies. § 5 presents the observed abundance patterns in the different X/O–O/H planes and compares them with galactic chemical evolution (GCE) models and local relations. § 6 explores possible physical explanations for the observed trends and assesses the impact of observational selection effects. Finally, § 7 summarises our main findings.

Observational data

We assembled a sample of star-forming galaxies at $z \sim1$–$7$ with JWST/NIRSpec spectroscopy. This work builds upon the work of C26, which focused on nitrogen, and substantially expands to include Ar, C, S and Ne. We selected objects with a signal-to-noise ratio S/N $> 5$ in both H$\alpha$ and H$\beta$ to ensure a robust H$\alpha$/H$\beta$ ratio and avoid introducing large uncertainties in the dust attenuation correction. Depending on the diagnostic of interest, we subsequently imposed a threshold of S/N $> 3$ on the other emission lines of interest: C iii] $\lambda$1907,1909, [O ii] $\lambda\lambda$3727,3729, [Ne iii] $\lambda$3869, [O iii] $\lambda$4363, [S ii] $\lambda\lambda$6716,6730, and [Ar iii] $\lambda$7135. When using the T$_e$ method, we additionally imposed S/N $> 3$ on the auroral [O iii] $\lambda$4363 line. Not all objects in the sample have simultaneous detections of all emission lines of interest. Therefore, the number of galaxies contributing to each abundance ratio measurement varies accordingly. As a

craven? result, the subsamples used for determining direct abundances of N/O, Ar/O, C/O, S/O and Ne/O analyses are not identical, and their sizes are reported in Table 1. The redshift distribution of each subsample is instead illustrated in Figure 1.

	Tracer	$N$	$\langle z\rangle$	$\alpha$	$\beta$	$\gamma$	$\sigma_{\mathrm{int}}$
N/O	N2O2	114	2.8	0.03 ±0.06	1.06 ±0.04	–	0.07
Ar/O	Ar3O3	90	2.6	−0.24 ±0.07	1.33 ±0.04	0.57 ±0.03	0.05
S/O	S23/R23	126	2.6	−0.46 ±0.09	1.53 ±0.12	0.59 ±0.07	0.08
Ne/O	Ne3O3	180	3.3	0.53 ±0.04	1.06 ±0.04	–	0.04

AGN contamination was assessed and excluded using the standard emission line diagnostics, namely [O iii]/Hβ versus [N ii]/Hα adopting the redshift-dependent calibration of Kewley et al. [94]. We also tested the diagnostic of Mazzolari et al. [115] based on [O iii] λ4363. This diagnostic yields a substantially different classification for our sample: approximately 35% of the galaxies classified as star-forming according to Kewley et al. [94] fall within the AGN region of the [O iii] λ4363-based diagnostic (although most lie in the proximity of the classification boundary). Nevertheless, we verified that this difference in sample selection does not affect the observational results of this work: although the alternative selection results in a smaller sample, all the observational trends of chemical abundance patterns remain clearly present. We therefore adopt the widely established Kewley et al. [94] diagnostic as our fiducial classification, while noting that our main conclusions are robust to the choice of AGN diagnostic.

The MARTA survey

The Measuring Abundances at high Redshift with the $T_{\mathrm{e}}$ Approach (MARTA; PID 1879, PI Curti) programme is a deep JWST/NIRSpec Micro Shutter Array (MSA) campaign targeting star-forming galaxies at Cosmic Noon ($z \sim2$–3). A key goal of the programme is the detection of the temperature-sensitive auroral line [O iii] λ4363 with the G140M grating ($R \sim1000$), enabling robust direct-method metallicity determinations. The total integration time is 31.9 hours with G140M and 7.4 hours with G235M. The MARTA survey is described in detail in Cataldi et al. [31], hereafter C25. We obtained approximately 110 spectra of galaxies with confirmed spectroscopic redshift.

JWST direct sample

We supplemented the MARTA data with a compilation of galaxies at $z > 1$ drawn from the literature, selecting additional galaxies which satisfy the same criteria in terms of emission line selection and S/N thresholds, obtained from the publicly available JADES DR4 catalogue (Curtis-Lake et al. [42]; Scholtz et al. [147]), and from the DAWN JWST Archive v4.4 (DJA; Valentino et al. [170])¹, a public repository of reduced JWST NIRSpec MSA spectroscopy processed with MSAExp (Brammer [22]), from which we adopted the standard emission-line flux catalogues directly². For all galaxies drawn from JADES DR4 and the DJA, the emission-line fluxes used here were measured from JWST/NIRSpec medium-resolution grating spectra (G140M, G235M, and G395M). The main programs used in this work are summarised in Appendix § 8.

We also compiled published oxygen, sulfur, argon and neon abundances of galaxies at cosmic noon from recent works of Rogers et al. [139]; Stanton et al. [154] and Foley et al. (2026), which we use for comparison with our analysis.

JWST strong-line sample

In addition to the direct-method sample described above, we assembled a larger sample of 514 galaxies satisfying the same selection criteria, except for the auroral-line detection requirement. This sample includes galaxies from MARTA, JADES, and the DJA with S/N above the adopted thresholds for all relevant emission lines. For these galaxies, we derived chemical abundances using the strong-line calibrations described in Sect. 3.3.

Local calibration sample

A local calibration sample is used to give a comparative view of all the main results presented in § 4 and § 5. It comprises low-redshift galaxies (Vale Asari et al. [169]; Curti et al. [39]; Khoram & Belfiore [96]) and H ii regions (Izotov et al. [79]; Guseva et al. [72]; Guseva et al. [71]; Berg et al. [15, 14]; Rogers et al. [137]; Zurita et al. [182]; Rogers et al. [136]) with direct-method oxygen, nitrogen, sulphur, argon abundances. A detailed description of the various subsamples is presented in C25 and C26. For all of these objects, ionic abundances were derived self-consistently by recovering the published line fluxes and recomputing electron temperatures following the methodology of § 3.2. This local sample spans a wide range in metallicity (12+log(O/H) ∈ [7, 9]) and constitutes the empirical anchor for strong-line calibrations.

For carbon, we instead compiled published abundance measurements directly from the literature. Specifically, we included a sample of 10 nearby metal-poor star-forming galaxies observed with HST (Senchyna et al. [149]), 11 additional LyC-leaking galaxies observed with HST (Izotov et al. [80]), and 32 nearby low-metallicity, high-ionization dwarf galaxies observed with HST (Berg et al. [16]; Berg et al. [12]), as well as carbon and oxygen abundances derived from recombination-line measurements in giant extragalactic H ii regions (Esteban et al. [55]; Esteban et al. [53]). Finally, we included measurements for Galactic and extragalactic H ii regions from García-Rojas & Esteban [63], and recombination-line abundances for H ii regions in the blue compact dwarf galaxy NGC 5253 obtained with VLT/UVES (López-Sánchez et al. [106]).

Chemical abundance measurements

Dust correction

Observed emission-line fluxes of all the samples described in § 2 are corrected for dust attenuation using the Balmer decrement. For each object, we fit all Balmer lines detected at S/N $>5$ (i.e. a subset of Hα, Hβ, Hγ, and Hδ depending on the specific target) to their theoretical Case B recombination values, assuming $T_{\mathrm{e}}=10^{4}$ K and $n_{\mathrm{e}}=100\ \mathrm{cm}^{-3}$, representative of the physical conditions in typical H ii regions and commonly adopted (Osterbrock & Ferland [127]). We adopt the Cardelli et al. [27] extinction law with $R_{\mathrm{V}}=3.1$ to derive the color excess $E(B-V)$ via a global $\chi^{2}$ minimisation.

As a consistency check, we also explored a simultaneous Markov Chain Monte Carlo (MCMC) fit of the dust attenuation and gas physical properties. In this approach, $E(B-V)$ is constrained by the observed Balmer-line ratios, while $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ are constrained independently by the available diagnostic lines: auroral lines for $T_{\mathrm{e}}$ and the [S ii] $\lambda6716/\lambda6731$ doublet for $n_{\mathrm{e}}$. When these diagnostics are not available, $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$ are fixed to their fiducial values. We found that the resulting $E(B-V)$ values change negligibly with varying physical conditions, with variations of less than 2%. Since this approach can only be applied to the galaxies with all the required diagnostic lines, we do not adopt it as our fiducial attenuation-correction method.

Employing alternative extinction laws, such as Calzetti et al. [25] or the Large Magellanic Cloud law (Gordon et al. [70]), produces negligible differences, on average of $\sim12\%$ on the computed $E(B-V)$ which translates in changes of less than 1% both on inferred temperatures and metallicities.

We also tested the impact of adopting a high-redshift nebular attenuation curve, as recently proposed by Reddy et al. [134], using their median attenuation law (their Eq. 22). This prescription and our fiducial Cardelli et al. [27] law yield very similar colour excesses, with a median relative offset of $\sim-10\%$. This difference has a negligible impact on the derived physical conditions and abundances, again resulting in variations of less than 1% in temperatures and abundances. It is also worth noting that $\approx45\%$ of our targets have $E(B-V)=0$, limiting the impact of the adopted attenuation curve on our results.

Direct method

For our sample of galaxies with auroral line detections (§ 2.1, § 2.2, § 2.4), ionic abundances were derived using the $T_{\mathrm{e}}$ method, which exploits the temperature sensitivity of collisionally excited emission lines: electron temperatures are inferred from auroral-to-nebular line ratios and then used to determine the ionic abundances.

The final products are 12 + log(O/H) from oxygen lines, and the abundance ratios log(N/O), log(Ar/O), log(C/O), log(S/O), and log(Ne/O) from the respective ionic ratios:

• N2O2 = log([N ii] $\lambda6584$/[O ii] $\lambda\lambda3727$, 29),

• Ar3O3 = log([Ar iii] $\lambda7135$/[O iii] $\lambda5007$),

• C3O3 = log([C iii] $\lambda1907$, 1909/[O iii] $\lambda5007$),

• S23/R23 = log $\frac{[\mathrm{S\ II}]\ \lambda\lambda6717,6731+[\mathrm{S\ III}]\ \lambda\lambda9067,9533}{[\mathrm{O\ II}]\ \lambda\lambda3727,3729+[\mathrm{O\ III}]\ \lambda5007}$,

• Ne3O3 = log([Ne iii] $\lambda3869$/[O iii] $\lambda5007$).

each corrected via the appropriate ICF.

Electron density and temperatures

The electron density $n_{\mathrm{e}}$ is derived from the [S ii] $\lambda\lambda6716,6730$ doublet ratio when both lines are available at S/N$>3$, which is the case for $\sim50\%$ of our high-$z$ targets. For the remaining objects, a fiducial value of $n_{\mathrm{e}}=400\ \mathrm{cm}^{-3}$ is adopted, corresponding to the median density measured for targets in our sample with a detected [S ii] doublet and consistent with the typical values reported for high-$z$ galaxies (e.g., Sanders et al. [142]; Stanton et al. [154]; Topping et al. [166]). As a consistency check, for targets in which the [S ii] doublet is not detected at 3σ, we also computed redshift-dependent densities using the parametrisation from Topping et al. [166]. Adopting this redshift-dependent parametrisation instead of the fiducial value changes the resulting electron temperatures by less than 1% and the inferred metallicities by less than 0.1%, so we retain the fiducial value throughout.

Ionic abundances are computed from dust-extinction-corrected fluxes, $T_{\mathrm{e}}$, and $n_{\mathrm{e}}$ using the PyNeb getIonAbundance function (Luridiana et al. [107]). Specifically, the dust-corrected line fluxes are converted into ionic abundances by dividing by the corresponding emissivities, which depend on the adopted $T_{\mathrm{e}}$ and $n_{\mathrm{e}}$. The electron temperature assigned to each ionic species is guided by the structure of the ionisation zones, although the mapping is not always one-to-one. A summary of the lines used and the atomic data adopted is provided in Appendix § 9.

In particular, we derive the electron temperature of the high-ionisation zone, $T_{\mathrm{e}}(\mathrm{O}^{2+})$, from the [O iii] $\lambda4363$/[O iii] $\lambda5007$ auroral-to-nebular line ratio. When both [S iii] lines are detected at S/N $>3$, we simultaneously constrain $T_{\mathrm{e}}(\mathrm{O}^{2+})$ and $n_{\mathrm{e}}$ using the getCrossTemDen routine in PyNeb, which takes the [O iii] $\lambda4363$/[O iii] $\lambda5007$ and [S ii] $\lambda6731$/[S ii] $\lambda6716$ ratios as input and return both the electron temperature and density; otherwise, we derive $T_{\mathrm{e}}(\mathrm{O}^{2+})$ from the [O iii] aur/neb ratio alone using getTemDen, adopting the density constraint described above.

For the low-ionisation zone, we assume that N$^{+}$, S$^{+}$, and O$^{+}$ arise mostly from a common region, again motivated by the proximity of their ionisation potentials. Accordingly adopt a single electron temperature, $T_{\mathrm{e}}(\mathrm{O}^{+})$, to derive N$^{+}$/H$^{+}$ and S$^{+}$/H$^{+}$ as well. When the [O ii] $\lambda\lambda7319,7329$ auroral lines are detected at S/N $>3$ ($\sim5\%$ of our targets), $T_{\mathrm{e}}(\mathrm{O}^{+})$ is derived directly from their auroral-to-nebular ratio; for the remaining targets, we instead infer $T_{\mathrm{e}}(\mathrm{O}^{+})$ from the fiducial $T_{\mathrm{e}}(\mathrm{O}^{+})$–$T_{\mathrm{e}}(\mathrm{O}^{2+})$ relation of C25, calibrated on high-redshift systems.

The case of intermediate-ionisation species such as S$^{2+}$, Ar$^{2+}$ and C$^{2+}$ is less straightforward, as their ionisation zones partially overlap with both O$^{+}$ and O$^{2+}$. We therefore tested two different approaches to estimate their ionic abundances: one adopting the directly measured $T_{\mathrm{e}}(\mathrm{O}^{2+})$, and the other using $T_{\mathrm{e}}(\mathrm{S}^{2+})$ as a proxy for the electron temperature of the intermediate-ionisation zone. Since $T_{\mathrm{e}}(\mathrm{S}^{2+})$ is not directly constrained by auroral-line measurements in our sample, we infer it from empirical calibrations from $T(\mathrm{O}^{2+})$ using the relation presented in Garnett [64], which has recently been shown by Rogers et al. [139] to provide an accurate description of cosmic-noon galaxies based on direct temperature measurements from the CEERS, Strom et al. [158] and AURORA samples. Applying this calibration, we find that S, Ar and C abundances derived using $T_{\mathrm{e}}(\mathrm{S}^{2+})$ are systematically higher than those obtained with $T_{\mathrm{e}}(\mathrm{O}^{2+})$, with a median offset of 0.09 dex. Nevertheless, although the absolute abundance scale is shifted, the observed trends remain unchanged and all the main conclusions of this work are unaffected. For this reason, we adopt $T_{\mathrm{e}}(\mathrm{O}^{2+})$ throughout the paper. While this choice may introduce a systematic offset in the absolute Ar$^{2+}$, S$^{2+}$, and C$^{2+}$ abundances, it relies on directly measured temperatures rather than on empirically inferred ones, thereby avoiding the additional uncertainties associated with calibration-based estimates.

For neon, Ne$^{2+}$ and O$^{2+}$ predominantly coexist in the highly ionised inner regions of the nebula, although their ionisation zones are not strictly coincident (see Table in Appendix § 9). Ne$^{2+}$ preferentially traces somewhat more highly ionised gas than O$^{2+}$, placing it between O$^{2+}$ and higher-ionisation species such as C$^{3+}$ and N$^{3+}$. The choice of the appropriate electron temperature is therefore not straightforward, since neither direct auroral-line constraints nor widely established empirical calibrations of $T_{\mathrm{e}}$ for C$^{3+}$ and N$^{3+}$ are currently available. Given the substantial overlap between the Ne$^{2+}$ and O$^{2+}$ ionisation zones, we therefore adopt $T_{\mathrm{e}}(\mathrm{O}^{2+})$ for Ne$^{2+}$.

Ionisation correction factors

Ionisation correction factors (ICFs) are applied to convert ionic abundances into total elemental abundances by accounting for unobserved ionisation stages. We adopt the prescriptions of Amayo et al. [2], derived from a large grid of photoionisation models selected to reproduce the properties of 985 giant H ii regions in spiral galaxies and blue compact galaxies. To assess the sensitivity of our results to this choice, we recalculated the abundances using the widely adopted ICFs of Izotov et al. (2006). The two prescriptions give broadly consistent results: the median offsets are $\lesssim0.06$ dex for all abundance ratios except Ne/O, for which the offset is approximately 0.1 dex. A similar difference for Ne/O was reported by Esteban et al. [54], who found mean values differing by 0.12 dex between the Amayo and Izotov prescriptions in a sample of local star-forming regions. In all cases, however, the overall abundance trends remain unchanged.

Strong-line method

For the high-$z$ sample with no auroral line detections (§ 2.3) we inferred chemical abundances using empirical strong-line calibrations. We determined oxygen abundances by inverting the high-redshift metallicity calibrations of C25 through $\chi^2$ minimisation, following the methodology described in C26. We then inferred the other abundance ratios using the strong-line relations calibrated in § 4.

Strong-line calibrations for N, Ar, S, and Ne

Strong-line diagnostics are essential tools for deriving abundance ratios in galaxies where auroral lines are undetected, which constitutes the majority of high-redshift objects. While most previous work has focused on calibrations for oxygen abundance or N/O, in this section we also present calibrations for Ar/O, S/O and Ne/O using the strong-line ratios N2O2, S23/R23, Ar3O3, and Ne3O3. These ratios are chosen to minimise the impact of ICFs on the derived abundance ratios. The total abundance of an element relative to oxygen is most robustly constrained when the observed ionic ratio traces ionisation zones that overlap as closely as possible with those of [O ii] and [O iii]; this condition is best satisfied when the IPs of the relevant ionic species are similar.

Our inferred strong-line calibrations are shown in Figure 2. Each panel displays the logarithm of an abundance ratio (N/O, Ar/O, C/O, S/O and Ne/O) as a function of the corresponding emission-line diagnostic. The figure also shows the possible secondary dependence of these diagnostics on O32 = $\log([\mathrm{O\,III}]5007/[\mathrm{O\,II}]\lambda\lambda3727,3729)$, used as a tracer of the ionisation parameter, as colour coded by the colour bar on the right. We first visually inspect the distribution of each abundance diagnostic as a function of O32 to assess whether a secondary dependence is present.

Columns report the number of galaxies with direct abundance measurements used in each fit ($N$), the best-fit coefficients of the relation $\log(X/\mathrm{O}) = \alpha+ \beta Y + \gamma\mathrm{O32}$, where $\mathrm{O32} = \log([\mathrm{O\,III}]\lambda5007/[\mathrm{O\,II}]\lambda\lambda3727,3729)$; $Y$ is the adopted strong-line ratio, while $\sigma_{\mathrm{int}}$ is the intrinsic scatter of the relation, i.e. the scatter after accounting for the observational uncertainties in dex. The $\gamma$ parameter is not fitted for the N/O and Ne/O relations, as we find no evidence for a dependence on O32 in either case; including this additional dimension does not provide any improvement in the intrinsic scatter. The C/O relation is not fitted owing to the limited sample size of only 22 objects.

Across the panels, the objects with the largest deviations from the fitted relations are predominantly local galaxies with very low O32, typically O32 $\lesssim-0.2$. These values lie well below those sampled by most of the high-$z$ galaxies, illustrating that the local and high-$z$ datasets probe different regions of this parameter space. Where their O32 ranges overlap, however, the two samples broadly follow the same abundance–line relations. For Ar/O and S/O, including O32 accounts for much of the apparent spread, although the most weakly ionised local objects remain outliers. Whether the relations also apply to high-$z$ galaxies at such low O32 cannot be tested with the present sample: such objects are absent from it and may be particularly difficult to detect with the emission-line requirements adopted here. No fit was performed for the relation between C/O and C3O3, given the small number of direct measurements available for this diagnostic (22 objects), and the substantial observed scatter.

The fits are performed with LINMIX (Kelly [91]), a Bayesian linear regression method that accounts for measurement uncertainties in both variables and provides an estimate of the intrinsic scatter. The details of the individual diagnostics are discussed in the following sections. The best-fit parameters performed on the high-$z$ objects are reported in Table 1.

N/O: the N2O2 diagnostic

The N2O2 ratio has long been recognised as a particularly tight tracer of N/O because N$^{+}$ and O$^{+}$ have very similar IPs at both the lower and higher ionisation thresholds (see Table B.1), so they likely coexist within the same ionisation zone of the nebula and only have a limited dependence on its ionisation structure (e.g., Florido et al. [57]). The relation between N/O and N2O2 (Figure 2, panel a) is remarkably tight ($\sigma_{\mathrm{int}} \simeq0.08$ dex) and displays a slope close to unity, indicating an almost linear mapping between the emission-line ratio and the abundance ratio.

Our best fit relation is consistent within the errors with the calibration derived using the slightly smaller dataset (92 objects) presented by C26 (slope = 1.04 ± 0.09 and intercept of −0.03 ± 0.06). In particular, the updated relation exhibits a moderately steeper slope, further strengthening the hint of redshift evolution that was already noted in C26. Indeed, this relation shows a clear dependence on the electron temperature: galaxies with higher $T_{\mathrm{e}}$ systematically lie toward the upper end of the relation, following a steeper trend. This difference is likely driven by the enlarged sample used in the present analysis, which includes a larger number of high-redshift, high-temperature galaxies from the DJA. Their inclusion shifts the best-fitting relation toward a steeper slope. Cameron et al. [26] also favours a steeper slope by studying the same relation focused on the JADES sample.

Ar/O: the Ar3O3 diagnostic

The Ar/O ratio, traced by Ar3O3, is shown in panel b of Figure 2. For argon, we choose Ar3O3 as tracer, because Ar$^{2+}$ (IP 27.6 eV) pairs with O$^{2+}$ (35.1 eV), as discussed in the previous section. The Ar3O3 diagnostic shows a clear dependence on O32, and our best-fit is in agreement with the one from Foley et al. (2026) who reported similar coefficients. After accounting for uncertainties, the intrinsic scatter is 0.05 dex, significantly lower than the 0.14 dex obtained from a linear fit based solely on Ar3O3.

S/O: the S23/R23 diagnostic

For sulphur, deriving the total abundance from low-ionisation lines alone can be problematic. In particular, the S2O2 ≡ $\log([\mathrm{S\,II}]/[\mathrm{O\,II}])$ diagnostic has been shown to provide a poor constraint on the S/O abundance ratio (e.g., [2, 54]). This limitation is related to the ionisation structure of sulphur. The ionisation stages relevant in H ii regions and star-forming galaxies span a broad range of IPs, with S$^{+}$ and S$^{2+}$ producing strong emission lines in the optical, while higher ionisation stages can also contribute significantly under typical nebular conditions. In particular, S$^{2+}$ can be further ionised to S$^{3+}$, whose ionisation potential of 34.8 eV lies within the range of energies readily provided by the ionising radiation field of H ii regions.

Amayo et al. [2] developed an ICF for sulphur based solely on [S ii] emission, but found substantial discrepancies when compared with ones obtained using both S$^{+}$ and S$^{2+}$. For this reason, we adopt instead the more widely used S23/R23 diagnostic, defined in § 3.2, which incorporates both low- and intermediate-ionisation species. This choice comes at the cost of reducing the available sulphur sample, as it requires the detection of both [S ii] and [S iii] lines. In our sample the number of objects decreases from 175 to 126. A discussion on the different sulphur diagnostics (S2O2, S3O3, and S23/R23), and their impact on the inferred S/O ratios, is presented in Appendix § 10.

Panel d of Figure 2 shows the relation between S/O and S23/R23. The intrinsic scatter of this relation is 0.08 dex, based on a sample of 126 objects. The best-fit relation significantly benefits from the inclusion of the dependence on O32, which reduces the intrinsic scatter by ∼40%. However, a substantial scatter is still observed, likely reflecting the complex ionisation structure of sulfur and the difficulty of accurately tracing its different ionic phases. In Appendix § 10, we present additional S/O calibrations based on S2O2 and S3O3 (≡ $\log([\mathrm{S\,III}]\lambda9071,9533/[\mathrm{O\,III}]\lambda5007)$). These diagnostics require fewer emission lines and may therefore be applicable to a broader range of observations, although they exhibit moderately larger intrinsic scatter than the relation presented here. We focus here on S23/R23 for consistency with the other panels, as it is the diagnostic adopted to infer S/O in our analysis.

Ne/O: the Ne3O3 diagnostic

Since Ne$^{2+}$ and O$^{2+}$ predominantly coexist in the highly ionised regions of the nebula (see discussion in § 3.2.1), Ne3O3 is a good tracer for the Ne/O ratio, as shown in panel e of Figure 2. A dependence of Ne/O on O32 is also visible, particularly in the full sample. However, when restricting the analysis to the high-$z$ sample, the Ne/O relation is already relatively tight in terms of intrinsic scatter. Including O32 as an additional parameter yields only a modest improvement in the intrinsic scatter (∼0.01 dex), which does not justify the introduction of a third parameter. We therefore retain the two-parameter relation for Ne/O.

Observed abundance patterns

In this section we examine the abundance patterns as a function of oxygen abundance. We compare the high-$z$ measurements with our local benchmarks and chemical evolution models for the Milky Way from Kobayashi & Taylor [102] and M31 from Arnaboldi et al. [8]; ); Kobayashi et al. [98]. FigFigure 3 provides an overview of the observed patterns. Their implications for chemical enrichment and galaxy evolution are discussed jointly in § 6.

Nitrogen

Figure 3a confirms the enhancement in N/O at low metallicities in high-$z$ galaxies with respect to their lower-redshift counterparts previously reported by C26 and Cameron et al. [26]. Our updated compilation contains 133 direct measurements (19 of which are from Rogers et al. [139]), compared with 92 in C26, owing to the continued expansion of the DJA archive, plus 477 objects from the strong-line compilation. The sample spans $7.0 \lesssim12 + \log(\mathrm{O}/\mathrm{H}) \lesssim8.8$, extending towards lower metallicities than the minimum of approximately 7.3 reached in our previous work. The expanded dataset confirms that high-redshift galaxies tend to have higher N/O than local galaxies at low oxygen abundance, with the excess becoming more pronounced at lower metallicity. The UV nitrogen emitters shown in Figure 3 occupy a more extreme region, with $\log(\mathrm{N}/\mathrm{O}) \gtrsim-0.8$. This distinction is also present in the analysis of Cameron et al. [26], who measured N/O in 588 JADES galaxies at $z \sim1.5$–7, including 40 with direct $T_{\mathrm{e}}$-based abundances. They also find a moderate N enhancement at low metallicity. Other recent observations instead yield less pronounced differences: using deep, medium-resolution rest-optical spectra of 25 typical star-forming galaxies at $z \sim3$ from the LyC22 survey, Schaerer et al. [145] find no significant N/O evolution within the scatter of local H II-region measurements.

We discuss the possible connection between these populations of N-enhanced galaxies and the additional constraints provided by carbon in § 6.

Argon

We collect 138 galaxies with direct-method Ar/O measurements (90 from our sample and additional 48 objects from Stanton et al. [154]; Foley et al. (2026); Rogers et al. [139]), complemented by 92 galaxies with strong-line abundances. Panel b of Figure 3 shows a decreasing trend of Ar/O with increasing oxygen abundance, particularly evident in the direct-method measurements. The strong-line measurements exhibit larger scatter and suggest a somewhat flatter behaviour at higher metallicities. A similar decrease in Ar/O with metallicity was reported locally by Esteban et al. [54], both for star-forming galaxies and H II regions, using abundances derived from either [Ar III] alone or [Ar III] together with [Ar IV]. Recently, Scholte et al. (2026) identified a negative correlation and derived a linear relation from approximately 50,000 star-forming galaxies at $z \sim0$ in the DESI survey (DESI Collaboration et al. [45]), with which our overall median trend is consistent within the uncertainties.

Recent high-$z$ measurements of Ar/O broadly agree with the sequence traced by our sample. The CECILIA galaxies of Rogers et al. [139], also included in our archival dataset shown in Figure 3, have subsolar Ar/O on average and show no clear dependence on O/H within their sample, but they align well with our trend and extend the available constraints towards higher metallicities. The EXCELS measurements of Stanton et al. [154] and the JADES stacks of Isobe et al. [75] are likewise consistent with this sequence.

In Figure 4, we examine the distribution of galaxies in the O/Ar–Ar/H plane and its dependence on stellar mass. Following the work of Bhattacharya et al. [19], we divide the sample into two stellar-mass bins, $\log(M_\star/M_\odot) \leq9$ and $\log(M_\star/M_\odot) > 9$, and compute the median O/Ar in adaptive bins of Ar/H, each containing at least five galaxies. The shaded regions indicate the 16th–84th percentile range of the observed O/Ar distribution within each bin. Despite the limited sample size, a clear mass-dependent separation emerges. At fixed Ar/H, the median sequence of higher-mass galaxies lies systematically $\sim0.1$–$0.2$ dex above that of lower-mass systems. The higher-mass sequence broadly follows the GCE tracks shown in the figure, although its last two bins at the highest Ar/H exhibit a downturn towards lower O/Ar. Lower-mass galaxies, spanning approximately $7.5 \lesssim\log(M_\star/M_\odot) \lesssim9$, occupy systematically lower O/Ar values over the overlapping Ar/H range.

Possible interpretations of the Ar/O behaviour and its mass dependence in the O/Ar–Ar/H plane are discussed in more detail in § 6.

Carbon

We measure C/O in 22 galaxies, restricted to systems at $z \gtrsim3$ for which the relevant rest-frame UV carbon lines fall within the observed spectral range. Because O III] $\lambda\lambda1666$ is generally unavailable, we infer C/O using C III] $\lambda\lambda1907,1909$ and [O III] $\lambda5007$. The large wavelength separation makes the C III]/[O III] ratio sensitive to differential dust attenuation. However, most galaxies in the sample have low $E(B-V)$, which limits the impact of the assumed attenuation curve. The largest value, $E(B-V) = 0.45$, corresponds to a correction of approximately a factor of seven in C III]/[O III] relative to the uncorrected ratio. This is a substantial correction, but only $\sim$5% of the galaxies have $E(B-V) > 0.3$, so corrections of this size affect a small fraction of the sample.

Figure 3c suggests an increase of C/O with oxygen abundance. This behaviour resembles that observed in local star-forming galaxies and H II regions ([66, 103, 55, 53, 16, 124, 12]). The trend is also broadly consistent with the chemical evolution models for the Milky Way and M31 shown in the figure, and the high-$z$ values lie within the uncertainties of the local C/O–O/H relation of Nicholls et al. [124]. However, the limited sample size prevents firm conclusions. The implications of this result, particularly in combination with the elevated N/O, are discussed in Sect. 6.

Sulfur

We collect 143 direct S/O measurements (126 from our sample and 17 from Rogers et al. [139]) in the high-$z$ sample, complemented by 227 galaxies with strong-line abundances. Panel d of Figure 3 shows a decreasing S/O trend with increasing oxygen abundance, broadly similar to that observed for Ar/O. The strong-line measurements follow a similar overall trend, but suggest a flattening at higher metallicities. These measurements dominate the sample above $12 + \log(\mathrm{O}/\mathrm{H}) \approx8.1$, as expected given that they are not subject to the detectability limitations imposed by weak auroral lines. The median S/O values are consistent within the uncertainties with the decreasing linear relations derived for local galaxies by Dors et al. [49] and Brazzini et al. [23], except for the highest-metallicity bin, which is dominated by strong-line measurements and lies above both relations. In this regime, the median S/O values tend to flatten, approaching the relations predicted by local GCE models and the empirical relation derived by Scholte et al. (2026). This flattening could be related to dust depletion, which becomes increasingly important in metal-rich environments. If oxygen is preferentially depleted onto dust grains relative to sulphur, the measured gas-phase S/O ratio would increase, potentially counteracting the declining trend observed at lower metallicities. The local literature presents a mixed picture. Decreasing relations have been reported in several studies ([173, 49, 48, 23, 129]), whereas most find approximately constant or increasing S/O ([65, 93, 14, 54] (Izotov et al. 2006; Scholte et al. 2026)). These differences highlight the importance of the adopted samples and abundance prescriptions. At high redshift, our measurements overlap well with those from CE-CILIA presented by Rogers et al. [139]. Their sample, of approximately 20 galaxies, does not establish a metallicity dependence on its own, but extends the available direct measurements to $12 + \log(\mathrm{O/H}) \sim8.8$.

Neon

Ne/O is the best-sampled abundance ratio in our analysis, with 188 direct measurements (180 from our sample and eight from Stanton et al. [154]), enabled by the widespread detection of [Ne iii] $\lambda3869$ in JWST spectra, complemented by 395 strong-line measurements. As shown in panel e of Figure 3, the two samples occupy broadly similar regions of the Ne/O–O/H plane, although the strong-line measurements exhibit larger scatter. The combined median Ne/O remains approximately constant up to $12+\log(\mathrm{O}/\mathrm{H}) \sim8.2$, with an intrinsic scatter of $\sim0.16$ dex, and shows no significant redshift dependence over $z\sim1$–7. At higher metallicities, where the strong-line sample extends the coverage to $12+\log(\mathrm{O}/\mathrm{H}) \sim8.6$, the median values suggest a mild increase of $\lesssim0.1$ dex.

The nearly constant Ne/O at lower metallicities is consistent with the EXCELS measurements of Stanton et al. [154] and the JADES stacks of Isobe et al. [75]. A weak increase towards higher O/H has also been reported in local studies, including Izotov et al. (2006), who find a variation of approximately 0.1 dex, as well as Esteban et al. [54] and Scholte et al. (2026). Similarly, the local relation of Scholte et al. (2026) and the chemical evolution models shown in Figure 3 predict only a modest variation across the metallicity range probed by our sample. The slight increase suggested by our strong-line measurements above $12+\log(\mathrm{O}/\mathrm{H}) \sim8.2$ is therefore qualitatively consistent with these expectations, although both the local relation and the models lie approximately 0.2 dex below our median measurements.

Part of the offset relative to the local relation may arise from differences in the adopted ionisation correction factors (ICFs). While Scholte et al. (2026) employ the prescription of Izotov et al. (2006), we adopt that of Amayo et al. [2]. As discussed in § 3.2.2, applying the former prescription to our sample would lower the inferred Ne/O by approximately 0.1 dex, reducing the discrepancy. The observed distributions nevertheless show partial overlap. The increase in Ne/O at higher metallicities could also be related to dust depletion. Since neon is not depleted onto dust grains, whereas oxygen can be, an increasing oxygen depletion at higher metallicities would lead to an apparent enhancement of the gas-phase Ne/O ratio. Interestingly, this increase occurs over approximately the same metallicity range where the S/O relation begins to flatten, suggesting that dust depletion could contribute to both features.

Discussion

The abundance patterns presented above provide complementary constraints on the enrichment histories of high-redshift galaxies. N/O is moderately enhanced at low oxygen abundance, while C/O remains broadly consistent with the local metallicity-dependent relation. S/O and Ar/O decrease with increasing O/H, and Ne/O shows no significant dependence on metallicity or redshift. Interpreting these observations jointly can help distinguish changes in stellar enrichment from the effects of gas flows and abundance methodology. A limitation is that the different elements are not measured in identical samples. In particular, the C/O sample is smaller and predominantly probes higher redshifts than the rest-optical measurements. The following discussion distinguishes the implications of the overall distributions from those of simultaneous abundance measurements in individual galaxies.

Constraints from nitrogen and carbon

JWST observations have identified galaxies with strongly enhanced nitrogen abundances, often selected through N IV] $\lambda1486$ and N III] $\lambda1750$ emission (e.g., Morel et al. [119]; Ji et al. (2026)). These extreme UV N-emitters must be distinguished from the moderate N/O enhancement found through rest-optical diagnostics in our sample and by Cameron et al. [26]. The distinction does not necessarily imply unrelated enrichment mechanisms, but differences in the strength, spatial extent, and duration of enrichment may be important.

Stacking analyses suggest that nitrogen enrichment may also extend beyond the rare, most extreme emitters. Hayes et al. [73], using approximately 1000 low-resolution PRISM spectra, report an increase in UV-based N/O at $z>4$. They argue that enrichment mechanisms invoked for individual N-emitters could operate more weakly in the wider star-forming population.

Several mechanisms have been proposed for nitrogen enhancement, including WR, very massive, and supermassive stars (Charbonnel et al. [33]; Nandal et al. [121]; Kobayashi & Ferrara [99]; Marques-Chaves et al. [109]; Watanabe et al. [175]), enrichment in dense young clusters associated with proto-globular cluster formation (D’Antona et al. [43]; Gieles et al. [68]; Berg et al. [13]; Ji et al. (2026)), bursty star-formation histories (Bhattacharya & Kobayashi [20]; McClymont et al. [116]), preferential loss of oxygen-rich ejecta (Rizzuti et al. [135]), and dilution by metal-poor gas (Stiavelli et al. [157]). We refer to C26 for a detailed discussion of these possibilities. Here, we focus on how carbon can provide additional constraints.

Subsolar C/O is common among metal-poor galaxies. High-$z$ measurements generally remain compatible with the local C/O–O/H relation (Arellano-Córdova et al. [7]; Isobe et al. [76]; Jones et al. [87]; Castellano et al. [30]; Topping et al. [168]; Morel et al. [119]). For example, Arellano-Córdova et al. [4] measured C/O at $z\sim8.5$ using the C III] doublet and found a value consistent with the scatter of lower-redshift metal-poor dwarfs. Similarly, Topping et al. [168] measured subsolar C/O in RXCJ2248-ID, compatible with enrichment dominated by low-metallicity massive stars. The increase of C/O towards higher O/H can reflect both metallicity-dependent carbon yields from massive stars and delayed enrichment by lower-mass stars. Stronger stellar winds at higher metallicity can enhance the carbon contribution from massive stars, while LIMS return carbon on longer timescales ([1, 28]). These contributions help explain the local increase in C/O, particularly above 12 + log(O/H) $\sim$ 8.0 ([124]).

Against this background, elevated N/O combined with approximately normal C/O at fixed O/H would indicate nitrogen enrichment without a comparable carbon enhancement. A possible explanation involves the different yields and lifetimes of AGB stars: lower-mass AGB stars, with initial masses of $\sim$ 2–4 $M_\odot$, can release the carbon previously brought to their envelopes through third dredge-up ([88, 90]); more massive AGB stars can instead undergo HBB, which converts carbon into nitrogen through the CNO cycle ([88, 171]). Their shorter lifetimes allow this nitrogen-enriched material to be returned to the ISM before the carbon contribution from lower-mass AGB stars becomes important. The relevant mass ranges and yields depend on metallicity and stellar-evolution prescriptions. This scenario does not require an episodic SFH, although the duration and shape of star formation determine the relative contributions of the different stellar populations. The EXCELS observations of Arellano-Córdova et al. [6] illustrate this possibility. Their two galaxies at $z \simeq5$ have elevated N/O and relatively low C/O, with carbon abundances consistent with CCSN yields. Comparison with the models of Kobayashi et al. [101] favours a bottom-heavy IMF in their analysis, but the uncertainties also permit a standard IMF.

Rapid enrichment by H-burning products provides another explanation, particularly for extreme, young N-emitters. Several systems have C/N and N/O approaching the abundance patterns expected for CNO-processed material ([76, 30, 148]). Potential sources include WN stars, very massive or supermassive stars, and tidal disruption events. By contrast, a substantial contribution from carbon-rich He-burning ejecta would also affect C/O. The observations of Berg et al. [13] provide direct evidence for WR stars in RXCJ2248-ID3, where enhanced N/O and low C/O are interpreted as signatures of WR-dominated enrichment. This establishes the relevance of WR enrichment in at least some systems, but does not imply that it accounts for the moderate enhancement in the wider star-forming population. The duration of the nitrogen-rich phase and the dilution of stellar ejecta into the surrounding gas are essential parts of this interpretation. The compilation of 45 N-emitters by Morel et al. [119] likewise shows elevated N/C, while C/O is generally comparable to local values at similar O/H. This favours enrichment by H-burning products without a corresponding systematic enhancement from carbon-rich He-burning material. The increasing incidence of N-emitters towards higher redshift further motivates enrichment associated with dense young stellar clusters, potentially containing WR, VMSs, SMSs.

Gas flows can modify these patterns, but different processes have distinct effects. In the limiting case of dilution by gas containing negligible metals, O/H decreases while N/O and C/O remain unchanged. A galaxy can therefore appear displaced from the local abundance-ratio relations without newly producing nitrogen. Conversely, preferential oxygen loss, with carbon and nitrogen retained, raises both N/O and C/O while leaving N/C unchanged. Elevated N/C without a comparable C/O enhancement would therefore be difficult to explain through oxygen loss alone. These simplified predictions do not exclude combinations of gas flows and stellar enrichment, but clarify why simultaneous C, N, and O measurements are valuable.

Sulfur and argon in the context of alpha-element enrichment

The similar behaviour of S/O and Ar/O motivates a joint discussion of these elements, given their related nucleosynthetic origins ([100]). Our sample exhibits decreasing Ar/O with increasing O/H, broadly consistent with the trends reported by Esteban et al. [54]; Scholte et al. (2026). Although some high-redshift studies do not identify a significant metallicity dependence within their individual samples, their measurements broadly overlap with the sequence traced here ([154, 139] (Foley et al. 2026)). The relatively high Ar/O values at low oxygen abundance raise the question of how a substantial argon contribution can develop while the gas remains metal-poor.

Argon is produced by both CCSNe and SNe Ia, whereas oxygen is predominantly supplied by massive stars ([78, 100]). The balance between these enrichment channels can therefore modify Ar/O. Bhattacharya et al. [18] explored an illustrative model involving intermittent star formation and substantial pristine gas accretion to explain several high-redshift galaxies with low O/Ar. Following an initial burst, continued enrichment by SNe Ia during a period of reduced SF can increase Ar/O (and hence lower O/Ar). Subsequent pristine gas accretion can lower O/H while initially preserving the abundance ratio Ar/O, allowing relatively Ar-rich gas to occur at low metallicity. Renewed SF from a second burst then supplies fresh oxygen, producing a loop in the abundance plane and potentially bringing the system back towards the standard enrichment sequence. Such histories can explain the positions of individual galaxies.

Additional information is however provided by the mass dependence observed in Figure 4: at a given Ar/H, lower-mass galaxies have systematically lower O/Ar. The statistical analysis presented above supports this association, suggesting that the distribution is not dominated by displacements of individual systems. A similar mass-dependent trend is identified in local SDSS galaxies by Bhattacharya et al. [19], who interpret it in terms of star-formation efficiency. In their framework, more massive galaxies have more efficient SF, so enrich their gas in oxygen rapidly through CCSNe, reaching higher metallicities before the delayed contribution from SNe Ia becomes substantial. Less efficient, more prolonged SF in lower-mass galaxies allows SNe Ia to contribute while Ar/H is still relatively low, resulting in lower O/Ar at a given Ar/H. Our observed mass dependence is qualitatively consistent with this interpretation. Importantly, Bhattacharya et al. [19] find that the mass-dependence becomes unclear when including galaxies that have undergone perturbations, such as metal-poor gas accretion or outflows. The persistence of a mass-dependent sequence in our sample – and its statistical significance – therefore suggests that large accretion-driven excursions do not dominate the population, although they could of course affect individual objects ([102, 19]).

Slower chemical enrichment, potentially sustained over an extended period, thus provides a plausible interpretation of the relatively high Ar/O values in our lower-mass galaxies. However, the correspondence between these mass-dependent sequences and the overall Ar/O–O/H trend remains to be established quantitatively. The Milky Way and M31 tracks provide illustrative comparisons, while models tailored to the masses, redshifts, and physical conditions of our sample are required to test whether this framework can simultaneously reproduce all the observed distributions.

Dust depletion provides an additional complication for sulfur. Oxygen is likely more strongly depleted than sulfur, so that increasing depletion at higher metallicity would raise gas-phase S/O rather than produce the decrease observed here [83]. More efficient dust formation at higher metallicity could strengthen this effect [62]. Under this picture, oxygen depletion is unlikely to explain the observed S/O trend. However, the flattening suggested by the strong-line measurements at $12 + \log(\mathrm{O/H}) \gtrsim8.0$ could be partly driven by increasing oxygen depletion, which would counteract the declining S/O trend. This could also help explain why the median values at higher metallicities approach the approximately constant relations predicted by local chemical evolution models. Nevertheless, uncertainties in sulfur depletion, ionisation corrections, and temperature estimates prevent us from identifying nucleosynthesis as the sole cause.

Finally, Ne/O offers a useful comparison. Oxygen and neon are predominantly produced by massive stars, with comparatively small contributions from AGB stars and SNe Ia, so their ratio is expected to be less sensitive to delayed enrichment [100] (Izotov et al. 2006). Its approximately constant behaviour up to $12 + \log(\mathrm{O/H}) \sim8.2$ is consistent with this expectation. At higher metallicities, however, the strong-line measurements suggest a mild increase in Ne/O, occurring over a similar metallicity range to the flattening observed in S/O. Both features could be related to increasing oxygen depletion onto dust grains, which would raise the measured gas-phase abundance ratios. This interpretation has also been proposed to explain the mild increase in Ne/O reported in some local samples, since neon is not incorporated into dust grains (e.g., Izotov et al. 2006). However, uncertainties in the ionisation corrections may also contribute to the observed behaviour. In particular, Esteban et al. [54] highlight limitations of Ne ICFs at low ionisation, as predicting the relative ionisation structures of Ne and O depends on the adopted ionising spectra and treatment of stellar opacities and winds [153] (Pérez-Montero et al. 2007).

In summary, the weak overall metallicity dependence of Ne/O, with a possible mild increase at higher O/H, and the absence of significant redshift evolution support the use of neon-based strong-line diagnostics. The bright [Ne III] λ3869 line remains accessible over a wide redshift range, making such diagnostics very useful for JWST observations [144, 146] (Pérez-Montero et al. 2007).

Selection effects in the abundance planes

In this section, we aim to answer the following question: if high-redshift galaxies obey all the currently established high-redshift scaling relations, but their relative abundances follow the local abundance–metallicity relations, would they be detectable with typical integrations of our archival JWST/NIRSpec data? We simulate a population of ‘regular’ high-$z$ galaxies, as described in detail in Appendix § 11, following a forward-modelling approach. We considered a galaxy as ‘detectable’ for a specific abundance measurement if both the auroral line [O III] λ4363 and the emission line associated with the abundance diagnostic under consideration (e.g. [Ar III] for Ar/O) are detected above the adopted signal-to-noise threshold of 3. The other nebular lines entering the abundance diagnostics (e.g. [O III] λ5007) are always brighter, so they do not drive the detectability. We focus this comparison on galaxies with direct-method abundances, since these drive the observed trends at low metallicity. In particular, we aim to assess whether the detectability requirement on [O III] λ4363 could introduce selection effects capable of producing or enhancing these trends.

The results are presented in Figure 5, which reproduces four of the abundance–metallicity planes shown in Figure 3, together with the simulated galaxy population. The contours enclose the 97% highest-density regions of the simulated galaxies satisfying the adopted detectability criteria for a 5-hour JWST/NIRSpec observation, representative of the median exposure time of the observing programmes currently included in the DJA. Different colours correspond to different ΔMS bins.

Several conclusions can be drawn from this figure. Most observed galaxies fall within the regions occupied by simulated galaxies that satisfy the adopted detection criteria; only a few lie outside the plotted contours. These cases are recovered when the contours are extended beyond the adopted 97% enclosed fraction towards the full detectable distribution; moreover, some of these galaxies were observed for longer than the representative 5-hour exposure assumed here, and would be detectable at correspondingly lower flux limits. Detectability depends strongly on a galaxy’s position relative to the star-forming main sequence. Simulated galaxies with ΔMS < 0.3 are detectable down to $12 + \log(\mathrm{O/H}) \sim7.8$–8.0, but are largely absent from the detectable population at lower metallicities. Below this value, the observed galaxies can be reproduced primarily by simulations in the ΔMS $\geq$ 0.3 bin, which extend down to $12 + \log(\mathrm{O/H}) \sim7.2$. Thus, with ∼5-hour JWST/NIRSpec exposures, direct-abundance samples at the lowest metallicities preferentially select galaxies with SFRs elevated above the main sequence.

Within the assumptions of the forward model, this selection does not appear to produce the abundance–metallicity trends by itself in all cases. Although the detection criteria restrict the metallicity range accessible in each ΔMS bin, the detectable simulated galaxies broadly follow the abundance relations assigned to the parent population. In the Ne/O panel, the local relation of Scholte et al. (2026), used to generate the simulated population, lies in median below the observed high-$z$ measurements. This offset possibly reflects differences in the adopted ICFs, since Scholte et al. (2026) adopt the prescriptions of Izotov et al. (2006), whereas we use Amayo et al. [2], which yield median offsets of ∼ 0.1 dex from the former, as previously pointed out in § 3.2.2. Nevertheless, the observed Ne/O values remain consistent with the simulated distribution within its scatter, and neither the observed nor the detectable simulated sample shows a clear Ne/O trend with metallicity. The detection cut mainly determines which part of that relation is accessible, with the lowest-metallicity galaxies appearing almost exclusively in the highest-ΔMS bin. A similar result is found for Ar/O: the simulated and observed populations show little discrepancy, and nearly all galaxies fall within the detectable region of the parameter space.

The N/O plane presents a different picture. At low metallicity, the detectable simulated galaxies preferentially have enhanced N/O, producing a distribution consistent with the trend observed in the direct-method sample. This suggests that selection effects could account for part of the N/O enhancement at low O/H, even if the underlying population follows the local abundance–metallicity relation. A weaker similar preference towards higher abundance ratios at low O/H is also visible in the S/O plane, hinting that selection may contribute to its observed trend.

Overall, these results suggest that detectability shapes the composition of the observed sample, particularly at low metallicity, but does not fully account for the observed abundance patterns of Ne/O and Ar/O, while potentially contributing to the trends observed for N/O and, to a lesser extent, S/O. The comparison remains qualitative because it adopts a single representative 5-hour exposure, while the observations span a range of integration times.

Another noteworthy feature is that essentially no galaxies with super-solar metallicities are predicted to be detectable, for each considered abundance ratio and $\Delta\mathrm{MS}$ bin. This limitation is not driven by the individual diagnostic lines, but rather by the requirement of detecting the temperature-sensitive [O iii] $\lambda4363$ auroral line. As illustrated in Appendix § 11, the [O iii] $\lambda4363$/[O iii] $\lambda5007$ ratio decreases steeply with increasing metallicity. More efficient cooling in metal-rich nebulae lowers the electron temperature, weakening the [O iii] $\lambda4363$ auroral line until it becomes difficult to detect with JWST/NIRSpec. Its measurement can also be complicated by the nearby [Fe ii] $\lambda4360$ emission line, whose contribution becomes particularly important when [O iii] $\lambda4363$ is faint [39, 151]. These limitations make direct-method abundances based on [O iii] $\lambda4363$ particularly challenging at high metallicity. Other auroral lines, however, do not decline as steeply: the [O ii] $\lambda\lambda7320,7330$ doublet may remain measurable in more metal-rich systems, offering an alternative route to constraining the electron temperature.

Conclusions

We have presented a systematic study of the abundance patterns of C, N, Ne, S, and Ar in star-forming galaxies at $z \sim1$--$7$, combining ultra-deep MARTA observations with archival JWST/NIRSpec spectroscopy. Our direct-method sample, based on detections of the auroral [O iii] $\lambda4363$ line, spans $7 \lesssim12 + \log(\mathrm{O/H}) \lesssim8.5$, allowing us to investigate multiple abundance ratios tracing different nucleosynthetic channels and enrichment timescales. Our main findings are summarised below:

1. In § 4, we derived new empirical strong-line calibrations for N/O, Ne/O, Ar/O, and S/O, anchored to the high-redshift direct-method sample. N2O2 and Ne3O3 provide particularly tight relations, with intrinsic scatters of $\sim0.07$ and $\sim0.04$ dex, respectively, while Ar/O and S/O show a significant secondary dependence on O3O2, which we include in their calibrations. These relations extend the abundance analysis to larger samples and higher metallicities, where auroral lines become difficult to detect. The limited C/O sample prevents us from deriving a robust calibration.

2. In § 5.1, we confirm the systematic N/O enhancement of high-redshift galaxies relative to local systems, particularly at low metallicities, previously reported in C26. The strong-line measurements support this result and extend the comparison to a larger sample and wider metallicity range.

3. In § 5.2, direct-method measurements reveal decreasing Ar/O with increasing oxygen abundance, broadly supported by strong-line measurements, which suggest a possible flattening at higher metallicities. We also identify a stellar-mass dependence in the O/Ar--Ar/H plane, with lower-mass galaxies exhibiting lower O/Ar at fixed Ar/H. This is qualitatively consistent with less efficient star formation and a larger relative contribution from delayed SNe Ia enrichment in lower-mass systems.

4. In § 5.3, the available C/O measurements tentatively suggest an increase with oxygen abundance, broadly consistent with the local relation and chemical-evolution models. This may reflect increasing contributions from delayed carbon production by intermediate-mass stars and AGB winds, although the limited sample prevents firm conclusions.

5. In § 5.4, S/O shows a decreasing trend with metallicity, similar to Ar/O. Strong-line measurements suggest a flattening above $12 + \log(\mathrm{O/H}) \approx8.1$, where median values approach local chemical-evolution predictions. Although oxygen depletion onto dust cannot explain the overall decline, it may contribute to the flattening at higher metallicities by increasing gas-phase S/O.

6. In § 5.5, Ne/O remains approximately constant up to $12 + \log(\mathrm{O/H}) \sim8.2$, with an intrinsic scatter of $\sim0.16$ dex and no significant redshift evolution over $z \sim1$--$7$. Strong-line measurements extend to $12 + \log(\mathrm{O/H}) \sim8.6$ and suggest a mild increase ($\lesssim0.1$ dex) at higher metallicities, potentially related to oxygen depletion onto dust. The relative stability of Ne/O supports the use of neon-based diagnostics at high redshift.

7. In § 6, we investigate observational selection effects through a forward model (Appendix § 11). The requirement of detecting [O iii] $\lambda4363$ and the relevant abundance-sensitive lines limits the accessible metallicity range, particularly below the star-forming main sequence. While selection primarily restricts the observable range for Ne/O and Ar/O, it preferentially selects galaxies with enhanced N/O at low metallicity, producing a distribution consistent with observations. A weaker effect is found for S/O, suggesting that selection may contribute to these abundance trends alongside intrinsic chemical enrichment.

Taken together these findings reveal distinct abundance--metallicity behaviours among elements with different nucleosynthetic origins. The S/O and Ar/O trends and the stellar-mass dependence of O/Ar provide new constraints for chemical-evolution models; our forward-modelling highlights the need to account for observational selection when interpreting abundance patterns, especially for N/O at low metallicity. Larger samples with homogeneous sensitivity and GCE models tailored to high-$z$ galaxies will be essential to disentangle nucleosynthetic processes, enrichment timescales, and observational biases.

Acknowledgements. FB, FM, AM, MG, IL, CM, AF, EB, GC acknowledge support from the INAF Fundamental Astrophysics programme 2022, 2023 and 2024. IL, FB, and AM acknowledge support from PRIN-MUR project “PROMETEUS” financed by the European Union – Next Generation EU, Mission 4 Component 1 CUP B53D23004750006 and C53D2300800006. GC acknowledges the COSMOWEBB project “A new era for cosmology: exploiting the JWST revolution”, funded by the Italian National Institute for Astrophysics (INAF) through the PRORIS programme (PROgramma di RIcerca Spaziale di base), CUP C83C26000300005. GV acknowledges financial support from the Italian National Institute for Astrophysics (INAF) under the IAF - Astrophysics Fellowships in Italy grant CUP C59J21034720001 - “AD MAJORA”. Some of the data products presented herein were retrieved from the Dawn JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center (DAWN), which is funded by the Danish National Research Foundation under grant DNRF140.

JWST archival data sources

Survey / PI	Program ID	Reference
JADES	1180, 1181, 1210, 1286, 1287, 3215	D’Eugenio et al. [46]
GTO (M. Stiavelli)	1199	Stiavelli et al. [156]
GTO (G. Rieke)	1207	—
GTO WIDE	1212, 1215	Maseda et al. [111]
GTO (C. Willott)	1222	Christensen et al. [36]
GLASS-ERS	1324	Mascia et al. [110]
CEERS-ERS	1345	Finkelstein et al. [56]
(B. Maseda)	1671	Maseda et al. [112]
Bluejay	1810	Belli et al. [10]
LyC-22 (D. Schaerer)	1869	—
AURORA	1914	Shapley et al. [150]
SUSPENSE	2110	Slob et al. [152]
(D. Stark)	2478	Topping et al. [167]
CECILIA	2593	Strom et al. [158]
SMACS-0723 ERO	2736	Pontoppidan et al. [132]
GTO (M. Stiavelli)	2758	—
JADES Ultra-Deep	3215	Eisenstein et al. [51]
ASPIRE (F. Wang)	3325	—
EXCELS	3543	Carnall et al. [29]
DEEPDIVE	3567	Ito et al. [77]
(E. Nelson)	4106	—
RUBIES	4233	de Graaff et al. [44]
(J. Antwi-Danso)	4318	—
(S. Fujimoto)	4762	—
GLIMPSE (S. Fujimoto)	9223	—

Atomic Parameters

Ion	IP inf [eV]	IP sup [eV]
$O^{+}$	13.6	35.1
$S^{+}$	10.4	23.3
$N^{+}$	14.5	29.6
$O^{2+}$	35.1	54.9
$Ne^{2+}$	40.9	63.4
$Ar^{2+}$	27.6	40.7
$C^{2+}$	24.4	47.9
$S^{2+}$	23.3	34.8

Notes. Ionisation potentials are taken from Osterbrock & Ferland [127]. For each ion, we report both the lower IP, corresponding to the energy required to produce that ionisation stage from the preceding one, and the upper IP, corresponding to the energy required to further ionise it to the next stage.

Ion	Atomic data	Collision strengths
$O^{2+}$	Tayal & Zatsarinny [163]	Palay et al. [128]
$S^{2+}$	Froese Fischer et al. [60] cited by Podobedova et al. [131]	Tayal & Gupta [161]
$S^{+}$	Rynkun et al. [140]	Tayal & Zatsarinny [162]
$O^{+}$	Zeippen [179]; Wiese et al. [176]	Kisielius et al. (2009)
$Ne^{2+}$	Galavís, M. E. et al. [61]	McLaughlin & Bell [117]
$Ar^{2+}$	Munoz Burgos et al. [120]	Munoz Burgos et al. [120]
$C^{2+}$	Glass [69]; Nussbaumer & Storey [126]; Wiese et al. [176]	Berrington et al. [17]
$N^{+}$	Froese Fischer & Tachiev [59]	Tayal [160]

Alternative sulfur diagnostics

Here we explore alternative strong-line calibrations for S/O. As discussed in § 4, sulfur abundance determinations can be affected by dust depletion and uncertainties in the ionisation correction. In particular, while the available spectra provide access to S$^{+}$ and S$^{2+}$ through the [S II] and [S III] lines, respectively, they do not directly constrain S$^{3+}$, whose [S IV] emission line at 10.5 µm lies outside their wavelength coverage. This ion may contribute to the total sulfur abundance in highly ionised nebulae ([2] (Izotov et al. 2006)), since the energy required to ionise S$^{2+}$ to S$^{3+}$ is approximately 34.9 eV, nearly identical to the 35.1 eV required to ionise O$^{+}$ to O$^{2+}$ ([11]). It is therefore useful to assess how different combinations of sulfur and oxygen lines trace S/O.

We first compare the S/O abundances inferred using S2O2, which involves only low-ionisation lines, with those obtained using S23/R23, which includes both observed sulfur and oxygen ionisation stages (Figure C.1). The two estimates show a systematic offset: in the high-redshift sample, S23/R23 yields lower S/O values than S2O2, whereas the local comparison sample exhibits the opposite behaviour. For both samples, the departures from the one-to-one relation appear to depend on excitation, as traced by O32. In particular, the highly excited systems, predominantly belonging to the high-redshift sample, show substantially lower S/O when inferred using S23/R23.

These differences may reflect limitations of the adopted ionisation corrections. Over the range of $W = \mathrm{O}^{+}/(\mathrm{O}^{+}+\mathrm{O}^{2+})$ probed by our samples, both corrections based on Amayo et al. [2] vary weakly. However, the correction applied to S2O2 has a higher normalisation, accounting for the contribution of ionisation stages not traced by the low-ionisation lines. By contrast, the correction applied to S23/R23 is smaller, since this diagnostic already includes S$^{+}$ and S$^{2+}$. This difference in normalisation contributes to the offset between the inferred abundances, but does not by itself explain its dependence on excitation or its reversal between the local and high-z samples. A possible explanation is that S$^{3+}$ contributes more substantially in the highly excited high-redshift systems than is accounted for by the adopted correction to S23/R23. A harder ionising spectrum could increase the fraction of sulfur in this unobserved stage, leading to an underestimate of S/O from the observed S$^{+}$ and S$^{2+}$ abundances if the correction is insufficient. The excitation-dependent discrepancy is consistent with this interpretation, although the comparison alone cannot establish the S$^{3+}$ contribution or determine which diagnostic provides the more accurate total abundance.

We next compare three strong-line diagnostics, S23/R23, S2O2, and S3O3, which involve either both observed sulfur ionisation stages or each stage individually. Their relations with S/O are shown in Figure C.2. In all three cases, the reference S/O abundances used to calibrate the diagnostics are derived from both the [S II] and [S III] lines, together with the adopted ionisation correction. For this comparison, we adopt linear relations without the additional dependence on O32 included in the fiducial calibration. The points are nevertheless colour-coded by log(O32) to illustrate how excitation contributes to the observed dispersion.

We fit the combined MARTA and JWST samples using LinMix, adopting the functional form

and accounting for measurement uncertainties in both coordinates and intrinsic scatter in log(S/O). The local comparison

Ratio	$\alpha$	$\beta$	$\sigma_{\mathrm{int}}$
S23/R23	$-0.94 \pm0.09$	$0.88 \pm0.02$	$0.124$
S2O2	$-0.89 \pm0.07$	$0.74 \pm0.07$	$0.135$
S3O3	$-1.07 \pm0.06$	$0.52 \pm0.06$	$0.135$

sample is shown for reference but is not included in the fits. The resulting coefficients and intrinsic scatters are listed in Table C.1.

These relations provide alternative estimates of S/O. In particular, S2O2 and S3O3 may be useful when only [S II] or [S III], respectively, is available, together with the corresponding oxygen lines. Their intrinsic scatters of approximately 0.12–0.14 dex should be taken into account when applying these calibrations. The colour coding also reveals a pronounced dependence on O32 in the S3O3 and S23/R23 relations, suggesting that including an explicit excitation term can improve these calibrations, as adopted for our fiducial relation. By contrast, S2O2 shows no similarly clear dependence on O32, although appreciable scatter remains.

Finally, it’s worth noticing that the MARTA galaxies occupy an intermediate locus between the local comparison sample and the highest-z systems, both in excitation and in the relative S/O estimates obtained from different diagnostics. This suggests that the discrepancy between the diagnostics may vary continuously with ionisation conditions, with typical Cosmic Noon galaxies bridging the two regimes. Such a pattern would be consistent with a gradual change in the ionic fractions and, potentially, in the hardness of the ionising spectrum, rather than a sharp distinction between local and high-z systems.

Modelling observational selection effects

The abundance trends presented in § 5 are based on a very heterogeneous sample, from a variety of observing programs (see Appendix § 8). For a galaxy to appear in any of the four diagnostic planes, its emission lines must be detected above the survey sensitivity threshold. Since the S/N cuts on the emission lines effectively impose a selection on the strength of the nebular emission, and hence on the SFR, the resulting sample may be biased with respect to the underlying galaxy population.

Here, we describe how we developed a forward-modelling framework to quantify the impact of selection effects on each of the four abundance diagnostics, and to assess which elements provide the most robust constraints on intrinsic chemical enrichment.

Forward model description

We construct a synthetic population of galaxies at a fixed redshift over a grid of stellar masses spanning $10^{7}$–$10^{11}\ M_{\odot}$ sampled uniformly in logarithmic space with bin width of 0.1 dex. This range is slightly broader than that of our observational sample, where the bulk of our objects reside at $M=7.5$–$10.5\ M_{\odot}$. We draw 300 Monte Carlo realisations per mass bin. For each physical quantity derived step-by-step in the code, we explicitly account for uncertainties by perturbing the value predicted by the various scaling relation discussed below with Gaussian noise, centred on the relation and scaled by its intrinsic scatter.

Star-formation rate. Each galaxy is assigned a star formation rate according to the empirical MS of Di Cesare et al. [47], calibrated on star forming galaxies at $z\sim4$–5, which gives a better representation to our sample than the main sequence presented by Popesso et al. [133]. An intrinsic scatter of $\sigma_{\mathrm{MS}}=0.37$ dex, as reported by Di Cesare et al. [47], is included.

Gas-phase metallicity. Oxygen abundances are drawn from the mass–metallicity relation (MZR). The parametrisation is taken from Curti et al. [40], adopting the $z=3$–6 relation over the $z=6$–10 one, as it better matches the redshift distribution of our sample. The relation has a reported intrinsic scatter of 0.07 dex. As a consistency check, we also implemented the MZR from Sanders et al. [141], adopting the $z=3.3$ relation over the $z=2.3$ one for the same reason. We find consistent results between the two MZR prescriptions.

Emission-line fluxes. The intrinsic H$\beta$ flux is derived from the SFR via the Kennicutt [92] calibration (renormalised to a Chabrier IMF). The [O II] and [O III] emission-line fluxes are predicted from the gas-phase metallicity O/H using the high-redshift strong-line calibrations of C25 for R2 and R3, respectively. We also tested the calibrations of Sanders et al. [143] and found that adopting them does not produce any significant differences in the resulting trends.

Finally, the [O III] $\lambda4363$ auroral-line flux is predicted from the galaxy metallicity by adopting a recalibrated relation between the [O III] $\lambda4363$/[O III] $\lambda5007$ line ratio and O/H. This relation has been derived in this work using the compilation of direct-method measurements described in § 2, including both low- and high-$z$ datasets. The resulting calibration is shown in Figure D.1, fitted with a second-order polynomial of the form:

where $x=12+\log(\mathrm{O/H})-8.69$. The best-fitting coefficients are

$a=-0.56\pm0.05,$

$b=-1.79\pm0.08,$

$c=-2.79\pm0.03.$

The abundance ratios X/O for each element are assigned from local empirical log(X/O)–metallicity relations (see below), and the corresponding emission-line fluxes are obtained by inverting the strong-line calibrations of § 4.

Abundance ratio and metallicity relations. We use the local abundance relations as a baseline for the intrinsic properties of the synthetic galaxies. Specifically, we assume that, before observational selection, galaxies at high redshift follow the same log(X/O)–$12+\log(\mathrm{O/H})$ relations as the local calibration sample described in § 2.4. We then apply the observational selection to this synthetic population and ask whether it can reproduce the trends measured in our high-$z$ sample. This test therefore assesses the effect of selection under the assumption of no intrinsic evolution in the abundance relations.

!Relation between the [O III] auroral-to-nebular line ratio, [O III] λ4363/[O III] λ5007, and the gas-phase metallicity bbox=0.12,0.07,0.48,0.29

For each element, we derive the baseline relation from robust linear fits to the local sample in metallicity bins of width 0.3 dex. For nitrogen, we adopt the relation from C26:

with $Z_0=-1.385$, $(\mathrm{O/H})_0=4 \times10^{-4}$, $\gamma=2.8$, $\beta=2.3$, and an intrinsic scatter of 0.2 dex.

For sulphur, argon, and neon, we adopt the recent calibrations of Scholte et al. (2026), which are derived from a large, statistical sample of $\sim49{,}000$ local galaxies and therefore provide a suitable reference for comparison with our high-$z$ measurements.

Dust attenuation. Observed fluxes are obtained from intrinsic values by applying a Cardelli et al. [27] extinction law ($R_V=3.1$) with $E(B-V)$ drawn from the mass–$E(B-V)$ relation of Woodrum et al. [177]. The scatter on $E(B-V)$ is propagated from both the intrinsic dispersion of the relation ($\sigma_{E(B-V)}=0.13$) and the uncertainty in stellar mass ($\sigma_{\log M_\star}=0.37$). Each line is attenuated at its rest-frame wavelength independently.

Detection thresholds

Detection thresholds for JWST/NIRSpec are computed following the approach of Zhu et al. [181], adapted to the medium-resolution gratings. We use the publicly available NIRSpec sensitivity curves provided by the JWST User Documentation³, which give the continuum flux density required to achieve a signal-to-noise ratio of 10 in a 10 ks exposure (see also Giardino et al. [67]). We assume that the emission lines are spectrally unresolved, such that their observed width is set by the instrumental resolution⁴. The corresponding $1\sigma$ noise level for a 10 ks exposure is therefore $\sigma(10\,\mathrm{ks})=f_\lambda/10$, which is then rescaled to an arbitrary exposure time $t$ assuming background-limited observations,

The minimum detectable integrated line flux at a significance level of $N\sigma$ is then $F_{\mathrm{lim}}=N\,\sigma(t)\,\Delta\lambda$, where we adopt $N=3$. Following the assumption of spectrally unresolved line, the observed line width $\Delta\lambda$ is $\Delta\lambda=\lambda_{\mathrm{obs}}/R(\lambda_{\mathrm{obs}})$, where $R(\lambda_{\mathrm{obs}})$ is the wavelength-dependent spectral resolving power of NIRSpec. For each emission line, we evaluate this width at the wavelength that gives the most conservative flux limit within the relevant grating coverage, using the tabulated resolution curves for G140M and G235M from the JWST User Documentation⁵.

With this information, we compute the expected signal-to-noise values for [N II] λ6583, [S II] λ6716,6730, [Ar III] λ7135, [Ne III] λ3869 and [O III] λ4363 emission lines, assuming representative integration times of 1, 5, 10, and 30 hours, as shown in Figure D.2. The resulting flux limits are also reported in Table D.1. We do not perform this analysis for carbon, as the currently observed sample is too small to draw robust conclusions about any trends that might emerge in the simulated population.

Figure D.2 provides a sanity check of our forward model: the simulated flux distributions broadly reproduce the observed line fluxes, and the adopted detectability thresholds are consistent with the region occupied by the observed sample. This agreement supports the use of the simulated population to assess how line detectability affects the observed abundance distributions.

Line	1h	5h	10h	30h
[N II] λ6583	5.6	2.5	1.8	1.0
[S II] λ6716,6730	5.1	2.3	1.6	0.9
[Ar III] λ7135	6.6	2.9	2.1	1.2
[Ne III] λ3869	9.6	4.3	3.0	1.8
[O III] λ4363	8.1	3.7	2.5	1.6

Notes. All fluxes are in units of $10^{-19}$ erg/s/cm$^2$.