Chemical Abundance Methodology

Comprehensive documentation of the emission-line abundance derivation pipeline implemented in jwspecabund. The methodology follows the six-step procedure of Berg et al. (2025) for UV nitrogen abundances, extended with strong-line calibrations from Sanders et al. (2025) and a Bayesian forward model following Cullen et al. (2025).

1. Overview

Three independent methods are available, orchestrated by compute_abundances() in jwspecabund._core:

Method	Key reference	When to use
Direct T_e	Berg et al. (2025)	Auroral line detected (e.g. [OIII] 4363)
Strong-line	Sanders et al. (2025)	No auroral detection; uses O3, O2, R23, O32
Bayesian forward model	Cullen et al. (2025)	MCMC/nested sampling over (T_e, n_e, ionic abundances)

All three methods begin with the same dust correction step. The direct T_e method is documented in the most detail below, as it is the primary approach for galaxies with UV spectroscopy.

2. Dust Correction

Source: jwspecabund.dust

2.1 Attenuation from the Balmer decrement

The dust attenuation A_V is derived from a Balmer line ratio (typically Hgamma/Hbeta, as Halpha may fall outside the observed wavelength range at high redshift). The observed ratio R_obs is compared to the intrinsic Case B recombination ratio R_int:

A_V = 2.5 × log10(R_obs / R_int) / (f_den − f_num)

where f_num and f_den are the normalised attenuation curve values A(lambda)/A_V at the numerator and denominator wavelengths. Negative A_V values are clamped to zero.

Intrinsic Case B ratios at T_e = 10^4 K, n_e = 100 cm^-3 (Osterbrock & Ferland 2006):

Ratio	Value
Halpha / Hbeta	2.86
Hgamma / Hbeta	0.468
Hdelta / Hbeta	0.259

2.2 Extinction / attenuation laws

Two laws are implemented:

Cardelli, Clayton & Mathis (1989): Milky Way extinction curve with R_V = 3.1. Polynomial fits in 1/lambda across IR, optical/NIR, and UV regimes. Used for individual galaxy spectra following Berg et al. (2025).
Salim et al. (2018) / Noll et al. (2009): Modified Calzetti et al. (2000) attenuation curve with a UV bump and power-law slope deviation. Default parameters: R_V = 3.15, delta = -0.35, bump strength B = 2.27.

2.3 Multi-Balmer A_V

When multiple Balmer lines are detected (Hgamma through H10), individual A_V values are derived from each line’s ratio to Hbeta and combined as an SNR-weighted average via compute_Av_multi_balmer(). The weights are the inverse-variance of each A_V estimate.

Excluded lines (blended at typical NIRSpec resolution):

Hepsilon (3971 Å): blended with [NeIII] 3968 Å (Δλ = 3 Å)
H8 (3890 Å): blended with HeI 3889 Å (Δλ = 0.4 Å)

Included lines:

Line	λ_rest (Å)	Intrinsic ratio to Hβ
Hgamma	4342.90	0.468
Hdelta	4102.89	0.259
H9	3836.48	0.0731
H10	3799.00	0.0530

Single-pair option. Passing balmer_pair=(numerator, denominator) to compute_abundances (or calling compute_Av_balmer_pair() directly) restricts the derivation to one chosen decrement — e.g. ("Ha", "HBETA") for Hα/Hβ — bypassing the SNR-weighted multi-line fit. This is useful when only one Balmer ratio is high-SNR and the fainter lines would otherwise bias or inflate the A_V error. The intrinsic Case-B ratio is taken from the ladder above (extended with Hα = 2.86 and Hβ = 1.00). A_V is still applied as a single point value to all draws; its formal error is reported on the result but not marginalised into the abundances unless A_V resampling is explicitly enabled (§2.4).

2.4 Dust correction of line fluxes

Each emission-line flux is corrected using:

F_corr = F_obs × 10^(0.4 × A_lambda)

where A_lambda = A_V × f(lambda) from the chosen attenuation law.

2.5 Error propagation

The uncertainty on A_V is propagated from the fractional error on the observed ratio:

sigma_R = R_obs × sqrt[(sigma_num / F_num)^2 + (sigma_den / F_den)^2]
sigma_AV = |dA_V/dR| × sigma_R
         = 2.5 / [R_obs × ln(10) × (f_den − f_num)] × sigma_R

3. Electron Temperature

Source: jwspecabund.direct

3.1 Primary diagnostic: [OIII] auroral-to-nebular ratio

The electron temperature of the high-ionisation zone, T_e(O2+), is derived from the collisionally-excited line ratio:

R_OIII = [OIII] 4363 / ([OIII] 5007 + [OIII] 4959)

PyNEB’s Atom("O", 3).getTemDen() inverts the ratio to solve for T_e. We use to_eval="(L(4363))/(L(5007)+L(4959))" with log-space root-finding (log=True, start_x=3.0, end_x=5.0) to ensure convergence at T_e > 25,000 K, which is typical of metal-poor high-z galaxies.

3.2 UV fallback diagnostic: O III] 1666 / ([OIII] 5007 + 4959)

When [OIII] 4363 is unavailable or has insufficient SNR, the O III] 1666 Å intercombination line provides an alternative T_e diagnostic:

R_1666 = O III] 1666 / ([OIII] 5007 + [OIII] 4959)

The 1666 Å line arises from the 5→2 transition with a 7.5 eV energy gap (compared to 2.8 eV for 4363), making this ratio more temperature-sensitive. PyNEB’s Atom("O", 3).getTemDen() is used with to_eval="(L(1666))/(L(5007)+L(4959))" and solver range extended to 250,000 K (start_x=3.0, end_x=5.4). The ratio is monotonically increasing with T_e, so the solution is unique.

This diagnostic is automatically used as a fallback in compute_abundances() when [OIII] 4363 is not detected but O III] 1666 and the [OIII] 5007+4959 nebular lines are available.

3.3 Secondary diagnostic: [NII] auroral-to-nebular ratio

When [NII] 5756 is detected, T_e(N+) is computed from:

R_NII = [NII] 5756 / [NII] 6585

using PyNEB’s Atom("N", 2).getTemDen().

3.4 T_e–T_e relations for zone temperatures

A single auroral line measures T_e in one ionisation zone. Empirical relations map T_e(O2+) (the high zone) to the intermediate and low zones:

Intermediate zone (O+, C2+, Si2+) — Garnett (1992):

T_int = [0.243 + t × (1.031 − 0.184 × t)] × 10^4 K

where t = T_high / 10^4.

Low zone (N+, S+) — two options:

DESI DR2 (arXiv:2601.02463, recommended):
```
T_low = 0.648 × T_high + 3270 K
```
Classical (Campbell, Terlevich & Melnick 1986; used in Garnett 1992):
```
T_low = 0.7 × T_high + 3000 K
```

4. Electron Density

Source: jwspecabund.direct, notebook 08

4.1 Density-sensitive doublets

Five diagnostic doublets span the full ionisation structure, each measured via PyNEB’s getTemDen() at the appropriate zone temperature:

Doublet	Ion	Zone	T_e used	PyNEB wavelengths
[S II] 6718/6732	S+	Low	T_low	6718, 6732
[O II] 3726/3729	O+	Low	T_low	3726, 3729
C III] 1907/1909	C2+	Intermediate	T_int	1907, 1909
Si III] 1883/1892	Si2+	Intermediate	T_int	1883, 1892
N IV] 1483/1486	N3+	High	T_high	1483, 1487
[Ar IV] 4711/4740	Ar3+	High	T_high	4711, 4740

4.2 Multi-phase density assignments (Berg et al. 2025, Section 4.1)

Each ionisation zone is assigned a representative density:

n_e,low: from [S II] (preferred) or [O II]
n_e,int: from C III] (preferred) or Si III]
n_e,high: from N IV] (preferred) or [Ar IV]
Fallback: n_e = 100 cm^-3 if no diagnostic is available

5. Oxygen Abundance — 12 + log(O/H)

Source: jwspecabund.direct.compute_ionic_abundances, compute_total_abundances

5.1 Ionic abundance formula

For any ion X^i+, the ionic abundance relative to hydrogen is:

X^i+/H+ = (F_line / F_Hbeta) × (epsilon_Hbeta / epsilon_line)

where epsilon denotes the volume emissivity computed by PyNEB at the appropriate (T_e, n_e).

For doublets (e.g. [O II] 3726+3729), the emissivities of both components are summed before dividing.

5.2 Hbeta emissivity

T_e <= 30,000 K: PyNEB RecAtom("H", 1).getEmissivity(Te, ne, wave=4861), based on the Storey & Hummer (1995) recombination tables.
T_e > 30,000 K: The Storey & Hummer tables do not extend beyond 30,000 K. We use the Aller (1984) Case B formula:
```
alpha_Hbeta_eff = 3.03 × 10^-14 × (T_e / 10^4)^(-0.874)   [cm^3 s^-1]
epsilon_Hbeta   = alpha_Hbeta_eff × h × nu_Hbeta
```
where h × nu_Hbeta = hc / lambda_Hbeta with lambda_Hbeta = 4862.68 A.

5.3 O2+/H+ from [OIII] 5007

Parameter	Value
Line	[OIII] 5007
PyNEB ion	`Atom("O", 3)`, wave=5007
Zone	T_high, n_e,high

5.4 O+/H+ from [OII] 3726+3729

Parameter	Value
Lines	[OII] 3726 + [OII] 3729 (summed)
PyNEB ion	`Atom("O", 2)`, waves=[3726, 3729]
Zone	T_low, n_e,low

If only the unresolved doublet is available (e.g. from prism data), the blended flux is used directly.

5.5 Total oxygen abundance

No ionisation correction factor is required for oxygen:

O/H = O+/H+ + O++/H+
12 + log(O/H) = 12 + log10(O/H)

The fraction O+/O quantifies the contribution of the low-ionisation zone, which feeds into the ICF calculations for other elements.

6. Nitrogen-to-Oxygen Ratio — log(N/O)

Source: jwspecabund.direct.compute_ionic_abundances, jwspecabund.martinez25_icf

Nitrogen abundance determination requires up to four ionic species, each measured in a different ionisation zone. An ionisation correction factor (ICF) then maps the observed ionic ratio to the total N/O.

6.1 N+ from [NII] 6585

Parameter	Value
Line	[NII] 6585
PyNEB ion	`Atom("N", 2)`, wave=6584
Zone	T_low, n_e,low

6.2 N2+ from NIII] 1749+1752

Parameter	Value
Lines	NIII] 1749 + NIII] 1752 (summed)
PyNEB ion	`Atom("N", 3)`, waves=[1749, 1752]
Zone	T_high, n_e,high

6.3 N3+ from NIV] 1483+1486

Parameter	Value
Lines	NIV] 1483 + NIV] 1486 (summed)
PyNEB ion	`Atom("N", 4)`, waves=[1483, 1487]
Zone	T_high, n_e,high

6.4 N4+ from NV 1239+1243

PyNEB has no atomic data for N V. We implement a manual three-level Li-like atom model:

Level 1: 2s ^2S_{1/2} (g = 2)
Level 2: 2p ^2P_{1/2} (g = 2), lambda = 1242.804 A, A_{21} = 3.37 × 10^8 s^-1
Level 3: 2p ^2P_{3/2} (g = 4), lambda = 1238.821 A, A_{31} = 3.40 × 10^8 s^-1

Effective collision strengths Upsilon(T) are interpolated from the tabulation of Aggarwal & Keenan (2004) / CHIANTI v10. Statistical equilibrium is solved for the three-level system, and the emissivity follows from:

epsilon = (n_upper / n_total) × A × h × nu / n_e

6.5 Direct-sum N/O (Berg et al. 2025, Eq. 5)

The raw ionic sum is:

(N/O)_raw = (N+ + N2+ + N3+ + N4+) / (O+ + O2+)

This is biased because the ionisation structure is non-uniform — the N3+/N ratio differs from O2+/O. An ICF is required to correct this (see Section 8).

6.6 ICF-corrected N/O

The final N/O is obtained by applying the best available ICF from Martinez et al. (2025):

log(N/O) = log10(ionic_ratio × ICF)

The ICF selection procedure is described in Section 8.2.

7. Carbon-to-Oxygen Ratio — log(C/O)

Source: jwspecabund.direct.compute_ionic_abundances, compute_total_abundances

7.1 C2+ from CIII] 1907+1909

Parameter	Value
Lines	CIII] 1907 + CIII] 1909 (summed)
PyNEB ion	`Atom("C", 3)`, waves=[1907, 1909]
Zone	T_high, n_e,high

7.2 C3+ from CIV 1548+1551

Parameter	Value
Lines	CIV 1548 + CIV 1551 (summed)
PyNEB ion	`Atom("C", 4)`, waves=[1548, 1551]
Zone	T_high, n_e,high

Note that CIV can have both stellar wind and nebular components. Only the nebular component should be used; broad stellar P Cygni profiles must be excluded.

7.3 C+ from CII] 2326

Parameter	Value
Lines	CII] 2324 + CII] 2326 (multiplet sum)
PyNEB ion	`Atom("C", 2)`, waves=[2323, 2325, 2326, 2327, 2328]
Zone	T_low, n_e,low

The CII] λ2326 multiplet (five components from 2323 to 2328 Å) traces C+ in the low-ionisation zone. When detected, the C+ ionic abundance is included in the carbon budget.

7.4 Total C/O

Two paths are available depending on whether C+ is detected:

With CII] 2326 (Garnett+1997 ICF):

When C+ is detected, the total C/O is computed using the Garnett (1997) ionisation correction factor:

C/O = ICF_C × (C2+/H+ + C3+/H+) / O2+/H+
ICF_C = (O+ + O2+) / O2+

This corrects for any remaining unobserved carbon ions by assuming the carbon and oxygen ionisation structures are similar.

Without CII] 2326 (Jones et al. 2023):

When C+ is not detected, C/O is computed as a direct ionic sum:

C/O = (C2+/H+ + C3+/H+) / (O2+/H+)

This assumes that C2+ and C3+ account for essentially all carbon in the O2+ zone, following Jones et al. (2023).

8. Ionisation Correction Factors

8.1 Izotov et al. (2006) — classical optical ICFs

Source: jwspecabund.icf

These are polynomial fits in x = O+/O_total, calibrated against photoionisation models. Used as the fallback when UV lines or log(U) are unavailable.

Element	Equation	Formula	Applied as
Nitrogen	Eq. 18	ICF_N = -0.825x^2 + 0.718x + 0.853	N/O = ICF_N × (N+/O+)
Neon	Eq. 19	ICF_Ne = -0.385x^2 + 0.405x + 0.335	Ne/O = ICF_Ne × (Ne2+/O2+)
Sulfur	Eq. 20	ICF_S = 0.013 + x(5.986 + x(-21.085 + x(26.462 - 11.282x)))	S/O = ICF_S × (S+ + S2+)/O
Argon	Eqs. 22-23	ICF_Ar = 0.157 + x(3.119 + x(-6.185 + 4.517x))	Ar/O = ICF_Ar × (Ar2+/O2+)

8.2 Martinez et al. (2025) — density-dependent ICFs

Source: jwspecabund.martinez25_icf

These ICFs are calibrated on Cloudy photoionisation models and account for the dependence of the ionisation correction on electron density, ionisation parameter, and metallicity. They are the recommended approach for UV N/O determinations (Berg et al. 2025, Section 4.2).

Ionisation parameter diagnostics

Before applying an ICF, the ionisation parameter log(U) must be estimated from one of two diagnostics:

N43 diagnostic (recommended):
```
log(N43) = log10(NIV] 1486 / NIII] 1750)
log(U_high) = f(log(N43), Z/Zsun; n_e)
```
Density-insensitive. Preferred when NIV] and NIII] are both detected. Valid range: -3.0 <= log(N43) <= 0.5.

O32 diagnostic (fallback):

log(O32) = log10([OIII] 5007 / [OII] 3727)
log(U_int) = f(log(O32), Z/Zsun; n_e)

Density-sensitive. Valid range: -1.0 <= log(O32) <= 2.5.

Both diagnostics use bicubic surface fits (Martinez et al. 2025, Table 3) interpolated between electron density grid points at log10(n_e) = [2, 3, 4, 5, 6]:

f(x, y) = A + Bx + Cy + Dxy + Ex^2 + Fy^2 + Gxy^2 + Hyx^2 + Ix^3 + Jy^3

The metallicity input is Z/Z_sun = 10^(12+log(O/H) - 8.69), where 12 + log(O/H)_sun = 8.69 (Asplund et al. 2009).

Overall validity ranges: 0.05 <= Z/Z_sun <= 0.40, -3.5 <= log(U) <= -1.0.

Five N/O ICFs (Table 4)

Each ICF corrects a specific ionic ratio to total N/O:

ICF	Ionic ratio	Notes
ICF 5	(N2+ + N3+) / O2+	Pure UV. Most robust at log(U) > -2.75 with corrections < 5%. Most preferred.
ICF 4	(N+ + N2+) / (O+ + O2+)	Mixed optical + UV. Robust across a wide range of conditions.
ICF 2	N2+ / O2+	UV single-ion. Overestimates at low log(U).
ICF 1	N+ / O+	Classical optical. Simplest but least accurate at high ionisation.
ICF 3	N3+ / O2+	Large corrections (100–15,000%); returned in log-space. Least preferred.

Selection priority

compute_NO_martinez25() selects the best ICF based on which ionic abundances are available, in order of decreasing preference:

ICF 5 — requires N2+, N3+, O2+
ICF 4 — requires N+, N2+, O+, O2+
ICF 2 — requires N2+, O2+
ICF 1 — requires N+, O+
ICF 3 — requires N3+, O2+ (last resort)

The final N/O is:

N/O = ionic_ratio × ICF(log(U), Z/Z_sun, n_e)

Inputs outside the calibration bounds

The Martinez et al. (2025) coefficients are polynomial fits over a bounded calibration domain (the validity ranges above). The bicubic surface diverges if it is evaluated by extrapolation, so an input that falls outside its range — log(O32), log(N43), Z/Z_sun, or the resulting log(U) — is rejected, not clipped: the surface returns NaN rather than a boundary value or an extrapolated number. This keeps unphysical, uncalibrated values out of the reported N/O entirely.

Because only the N/O ICF depends on log(U) and Z/Z_sun, this rejection affects N/O alone. O/H (direct T_e), T_e, and the other X/O ratios (C/O, S/O, Ne/O, Ar/O) do not use the Martinez ICF and are unaffected — see Section 11.3 for how this is handled in the Monte Carlo / posterior loops.

9. Strong-Line Calibrations

Source: jwspecabund.strong_line

When no auroral line is detected, 12 + log(O/H) is estimated from strong-line ratios using the simultaneous polynomial calibrations of Sanders et al. (2025).

9.1 Diagnostic ratios

Ratio	Definition	Lines required
O3	log10([OIII] 5007 / Hbeta)	OIII_5007, HBETA
O2	log10([OII] 3726+3729 / Hbeta)	OII_doublet, HBETA
R23	log10(([OIII] 5007 + [OII] 3726+3729) / Hbeta)	OIII_5007, OII_doublet, HBETA
O32	log10([OIII] 5007 / [OII] 3726+3729)	OIII_5007, OII_doublet

Each ratio R is modelled as a polynomial in Z = 12 + log(O/H):

R(Z) = c0 + c1 × (Z - 8.0) + c2 × (Z - 8.0)^2 + ...

9.2 Simultaneous fit

All available ratios are fit simultaneously by minimising the combined chi-squared:

chi^2(Z) = sum_i [(R_obs,i - R_model,i(Z))^2 / (sigma_obs,i^2 + sigma_cal,i^2)]

where sigma_cal = sqrt(sigma_R,fit^2 + sigma_R,int^2) accounts for both the fit residual and intrinsic scatter in the calibration. The best-fit Z is found via bounded scalar minimisation over 6.5 <= Z <= 9.0.

9.3 Monte Carlo uncertainty

The uncertainty on Z is estimated by perturbing each ratio by its observational error (1000 MC iterations) and re-fitting. The reported uncertainty is the standard deviation of the posterior.

10. Bayesian Forward Model

Source: jwspecabund.forward

An alternative approach following Cullen et al. (2025): free parameters (T_e, n_e, ionic abundances) are sampled via MCMC or nested sampling. For each parameter vector, collisionally-excited line emissivities from PyNEB and the Aller (1984) Hbeta recombination formula predict line flux ratios that are compared to observations via a Gaussian likelihood.

This method self-consistently accounts for the T_e–n_e coupling and non-Gaussian posterior distributions, but is more computationally expensive than the direct method.

11. Uncertainty Propagation

11.1 Monte Carlo approach (direct T_e method)

The full abundance pipeline is wrapped in a Monte Carlo loop (default N_mc = 1000 iterations):

Draw perturbed fluxes: F_perturbed ~ N(F_observed, sigma_F) for each emission line.
Re-derive all quantities: dust correction, T_e, n_e, ionic abundances, log(U), ICFs, total abundances.
Accumulate posterior arrays: {O/H, N/O, C/O, T_e, …}.
Report: median and 16th/84th percentile intervals.

11.2 MCMC posterior propagation

When MCMC chains are available from jwspecmcmc, posterior flux samples are drawn directly from the chains rather than assuming Gaussian errors. The dust correction is applied per-sample:

F_corr,i = F_obs,i × 10^(0.4 × k(lambda) × A_V,i)

where A_V,i is re-derived from the Balmer decrement in each sample. This propagates the full posterior covariance between line fluxes through to the final abundances.

11.3 Out-of-bounds resampling and the O/H–N/O decoupling

The Martinez et al. (2025) log(U)/ICF surfaces are only valid inside their calibration domain (Section 8.2). Inputs outside that domain are rejected (set to NaN), not clipped — see “Inputs outside the calibration bounds” in Section 8.2. In the uncertainty loops this is handled per draw so the reported sample counts stay meaningful:

Only N/O is gated. If a draw’s log(O32)/log(N43), Z/Z_sun, or resulting log(U) falls outside the bounds, log(U) is set to NaN for that draw. The Martinez N/O ICF then returns NaN, so the draw contributes nothing to N/O.
O/H and the other ratios are kept. O/H (direct T_e), T_e, and the X/O ratios that do not use the Martinez ICF (C/O, S/O, Ne/O, Ar/O) are computed and recorded for every drawn sample, in-bounds or not. They are never discarded on account of an out-of-bounds N/O.
Resampling to the requested count. The loop keeps drawing until n_mc (Gaussian MC) or n_posterior (MCMC pool) in-bounds N/O draws have been collected, so N/O is sampled to the requested depth rather than silently under-sampled. Solver failures and missing-line NaNs are not re-drawn (re-drawing cannot fix them) and count toward the target.
Different lengths are expected. Because out-of-bounds draws are re-drawn for N/O but retained for O/H, the O/H posterior generally holds more finite samples than the N/O posterior. The two arrays are independent; do not assume they are the same length.
Safety cap. To guarantee termination for an object centred outside the bounds (e.g. Z/Z_sun < 0.05, where acceptance ≈ 0), the total number of attempts is capped at 20× the requested count. If the target is not reached, a WARNING is logged ("... in-bounds N/O draws ...") and N/O is reported from whatever in-bounds draws were collected; O/H and the other ratios remain fully sampled.

Note

If you set the jwspecabund logger above WARNING, the under-sampling message is suppressed. To detect objects where N/O could not be sampled to the requested depth, inspect the posterior directly, e.g. np.isfinite(result.NO_posterior).sum().

11.4 Pre-computed emissivity grids

To avoid calling PyNEB ~10^6 times (N_mc × N_lines), emissivity grids are pre-computed as a function of T_e(high) on a fine grid (e.g. 5,000–80,000 K, 500 points). Each MC iteration interpolates into these grids, giving a factor ~100× speedup while preserving the parametric dependence on T_e.

12. Upper Limits for Non-Detected Ions

Source: jwspecabund._core._compute_continuum_rms_limits, _compute_ionic_upper_limits

When an emission line is not detected (flux consistent with zero or below the SNR threshold), we compute a 3σ upper limit on its flux and convert this to an ionic abundance upper limit. This is essential for constraining N/O and C/O when UV lines (e.g. N III], N IV], C IV) fall below the detection threshold.

12.1 Continuum-RMS flux upper limits

The flux upper limit is derived from the noise in the continuum around the expected line position, independent of the emission-line fit:

F_ul = n_sigma × RMS_cont × sigma_inst × sqrt(2π)

where:

RMS_cont: root-mean-square of the fit residuals (data minus model) in a window of ±5σ around the expected line centre, excluding the central ±2σ where the line would sit. If fewer than 3 pixels fall in this annulus, the window is expanded to ±10σ with ±3σ exclusion.
sigma_inst: the instrumental line width (in Angstrom), computed from the grating resolving power R as σ = λ_obs / (2.355 × R). When grating metadata is unavailable (e.g. spectra loaded from npz files), R is estimated from the pixel sampling via R_from_pixels(wave_um).
n_sigma: detection threshold, default 3.

This method is preferred over using the fit error because it measures actual noise at the line position rather than relying on bootstrap or MCMC error estimates, which can be inflated by unconstrained parameters (particularly for faint UV secondary doublet members).

12.2 Dust correction of upper limits

Flux upper limits are corrected for dust attenuation using the same A_V and attenuation law applied to detected lines:

F_ul,corr = F_ul × 10^(0.4 × A_lambda)

This ensures upper limits are on the same (intrinsic) flux scale as the detected line fluxes used for ionic abundances.

12.3 Doublet upper limits

For doublets where both members are undetected (e.g. N III] 1749+1752), the individual member upper limits are combined in quadrature:

F_ul,doublet = sqrt(F_ul,1^2 + F_ul,2^2)

This is the correct treatment for independent noise in each member, rather than linear summation which would overestimate the upper limit.

12.4 Ionic abundance upper limits

The flux upper limit is converted to an ionic abundance upper limit using the same emissivity-ratio formula as for detected lines (Section 5.1):

(X^i+/H+)_ul = (F_ul / F_Hbeta) × (epsilon_Hbeta / epsilon_line)

12.5 Abundance ratio upper limits

Upper limits on N/O and C/O are computed by summing detected ionic abundances with upper limits for non-detected ions, divided by total oxygen:

(N/O)_ul = (N_detected + N_upper_limits) / (O+ + O2+)
log(N/O)_ul = log10((N/O)_ul)

These are reported in the AbundanceResult.summary() output alongside the ionic upper limit details.

12.6 Fallback to fit errors

If the continuum-RMS method cannot be applied (e.g. the spectrum object is unavailable, or the residuals do not cover the line position), the pipeline falls back to the fit-error method:

F_ul = n_sigma × sqrt(sum(sigma_i^2))

where sigma_i are the per-member flux errors from bootstrap or MCMC. The ionic_ul_details dictionary records which method was used for each ion ("method": "continuum_rms" or "method": "fit_error").

13. References

Reference	Context
Aggarwal, K. M. & Keenan, F. P., 2004, A&A, 427, 763	N V effective collision strengths (CHIANTI v10)
Aller, L. H., 1984, Physics of Thermal Gaseous Nebulae (Reidel)	Hbeta Case B recombination formula for T_e > 30,000 K
Asplund, M. et al., 2009, ARA&A, 47, 481	Solar oxygen abundance: 12 + log(O/H)_sun = 8.69
Berg, D. A. et al., 2025, arXiv:2511.13591	Six-step UV nitrogen abundance procedure
Calzetti, D. et al., 2000, ApJ, 533, 682	Starburst attenuation curve (base for Salim+18)
Campbell, A., Terlevich, R. & Melnick, J., 1986, MNRAS, 223, 811	T_e–T_e relation: T_low = 0.7 × T_high + 3000
Cardelli, J. A., Clayton, G. C. & Mathis, J. S., 1989, ApJ, 345, 245	Milky Way extinction law
Cullen, F. et al., 2025, (EXCELS)	Bayesian forward model for emission-line abundances
DESI Collaboration, 2026, arXiv:2601.02463	Updated T_e–T_e relation: T_low = 0.648 × T_high + 3270
Garnett, D. R., 1992, AJ, 103, 1330	Classical T_e–T_e relations for intermediate/low zones
Garnett, D. R. et al., 1997, ApJ, 489, 63	C/O ionisation correction factor: ICF_C = (O+ + O++)/O++
Izotov, Y. I. et al., 2006, A&A, 448, 955	ICF equations 18–23 for N, Ne, S, Ar
Jones, T. et al., 2023	C/O from direct ionic sum: (C2+ + C3+) / O2+
Martinez, M. A., Berg, D. A. et al., 2025, arXiv:2510.21960	Density-dependent ICFs and log(U) diagnostics
Noll, S. et al., 2009, A&A, 499, 69	UV bump + slope modification to Calzetti curve
Osterbrock, D. E. & Ferland, G. J., 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (University Science Books)	Case B recombination ratios
Salim, S. et al., 2018, ApJ, 859, 11	UV slope and bump calibration for dust attenuation
Sanders, R. L. et al., 2025	Strong-line metallicity calibrations (O3, O2, R23, O32)
Storey, P. J. & Hummer, D. G., 1995, MNRAS, 272, 41	H I recombination emissivity tables (valid to 30,000 K)

Appendix: Zone Assignment Summary

Ion	Lines	Zone	T_e	n_e
O2+	[OIII] 5007	High	T_high	n_e,high
O+	[OII] 3726+3729	Low	T_low	n_e,low
N+	[NII] 6585	Low	T_low	n_e,low
N2+	NIII] 1749+1752	High	T_high	n_e,high
N3+	NIV] 1483+1486	High	T_high	n_e,high
N4+	NV 1239+1243	High	T_high	n_e,high
C+	CII] 2324+2326	Low	T_low	n_e,low
C2+	CIII] 1907+1909	High	T_high	n_e,high
C3+	CIV 1548+1551	High	T_high	n_e,high
S+	[SII] 6718+6732	Low	T_low	n_e,low
S2+	[SIII] 9069	Mid	(T_high+T_low)/2	n_e,low
Ne2+	[NeIII] 3869	High	T_high	n_e,high
Ar2+	[ArIII] 7136	Mid	(T_high+T_low)/2	n_e,low