jwspecabund

jwspecabund — Chemical abundance calculations for JWST spectra.

Computes oxygen abundances (O/H), nitrogen-to-oxygen ratios (N/O), and other element ratios from emission-line fluxes produced by jwspecfit or jwspecmcmc.

Three pathways are supported:

1. Direct T_e method — when auroral lines (e.g. [OIII] 4363) are detected, uses PyNEB for electron temperature / density diagnostics and ionic abundance calculations.

2. Bayesian forward model — samples ionic abundances, T_e, and n_e via MCMC/nested sampling, predicting line ratios with PyNEB CEL emissivities and the Aller (1984) Hbeta formula (Cullen+25 approach).

3. Strong-line calibrations — Sanders et al. (2025) simultaneous polynomial fit across O3, O2, R23, O32 diagnostics.

Example

>>> import jwspecfit, jwspecabund
>>> spec = jwspecfit.read_fits("spectrum.fits")
>>> result = jwspecfit.fit_lines(spec, z=6.0)
>>> abund = jwspecabund.compute_abundances(result, z=6.0)
>>> print(abund.OH)          # 12 + log(O/H)
>>> print(abund.method)      # "direct" or "strong_line"

class jwspecabund.AbundanceResult(method, OH, OH_err, NO=None, NO_err=None, CO=None, CO_err=None, Te_high=None, Te_high_err=None, Te_low=None, Te_low_err=None, ne=None, Av=None, Av_err=None, Av_posterior=None, ionic=None, OH_posterior=None, NO_posterior=None, CO_posterior=None, ratios_used=None, chi2=None, SO=None, SO_err=None, NeO=None, NeO_err=None, ArO=None, ArO_err=None, logU=None, logU_err=None, ne_low=None, ne_mid=None, ne_high=None, icf_method=None, NO_icf_name=None, lya_f_esc=None, lya_f_esc_err=None, lya_f_esc_posterior=None, lya_f_esc_details=None, excluded_lines=None, ionic_upper_limits=None, ionic_ul_details=None, NO_tiers=None, icf_values=None, failures=None, diagnostics=None, alt_results=None, _forward_result=None)

Bases: object

Container for a chemical abundance measurement.

Parameters:

• method (str) – "direct", "forward", or "strong_line".

• OH (float) – 12 + log(O/H).

• OH_err (float or tuple of float) – Symmetric error or (lo, hi) 68 % CI half-widths.

• NO (float or None) – log(N/O), if nitrogen lines available.

• NO_err (float or tuple of float or None) – Error on log(N/O).

• CO (float or None) – log(C/O), if UV lines present.

• CO_err (float or tuple of float or None) – Error on log(C/O).

• Te_high (float or None) – T_e(O++) in K (direct method only).

• Te_low (float or None) – T_e(O+/N+) in K (direct method only).

• ne (float or None) – Electron density in cm^-3 (direct method only).

• Av (float or None) – Dust attenuation A_V.

• ionic (dict or None) – Ionic abundance dict, e.g. {"O+/H+": val, "O++/H+": val, ...}.

• OH_posterior (np.ndarray or None) – Full posterior samples of 12+log(O/H) (MCMC input).

• NO_posterior (np.ndarray or None) – Full posterior samples of log(N/O).

• CO_posterior (np.ndarray or None) – Full posterior samples of log(C/O).

• ratios_used (list of str or None) – Diagnostic ratios used (strong-line method).

• chi2 (float or None) – Goodness-of-fit chi-squared (strong-line method).

• SO (float or None) – log(S/O) if [SII] and [SIII] available.

• NeO (float or None) – log(Ne/O) if [NeIII] available.

• ArO (float or None) – log(Ar/O) if [ArIII] available.

• excluded_lines (list of str or None) – Line names excluded by the per-line SNR filter.

• failures (dict or None) – Reasons why specific abundance ratios could not be computed, e.g. {"N/O": "no nitrogen ions detected"}.

• Te_high_err (float | tuple[float, float] | None)

• Te_low_err (float | tuple[float, float] | None)

• Av_err (float | None)

• Av_posterior (ndarray | None)

• SO_err (float | tuple[float, float] | None)

• NeO_err (float | tuple[float, float] | None)

• ArO_err (float | tuple[float, float] | None)

• logU (float | None)

• logU_err (float | tuple[float, float] | None)

• ne_low (float | None)

• ne_mid (float | None)

• ne_high (float | None)

• icf_method (str | None)

• NO_icf_name (str | None)

• lya_f_esc (float | None)

• lya_f_esc_err (float | tuple[float, float] | None)

• lya_f_esc_posterior (ndarray | None)

• lya_f_esc_details (dict | None)

• ionic_upper_limits (dict[str, float] | None)

• ionic_ul_details (dict[str, dict] | None)

• NO_tiers (dict[str, float] | None)

• icf_values (dict[str, dict] | None)

• diagnostics (dict[str, str] | None)

• alt_results (dict[str, AbundanceResult] | None)

• _forward_result (dict[str, Any] | None)

ArO: float | None = None

ArO_err: float | tuple[float, float] | None = None

Av: float | None = None

Av_err: float | None = None

Av_posterior: ndarray | None = None

CO: float | None = None

CO_err: float | tuple[float, float] | None = None

CO_posterior: ndarray | None = None

NO: float | None = None

NO_err: float | tuple[float, float] | None = None

NO_icf_name: str | None = None

NO_posterior: ndarray | None = None

NO_tiers: dict[str, float] | None = None

NeO: float | None = None

NeO_err: float | tuple[float, float] | None = None

OH: float

OH_err: float | tuple[float, float]

OH_posterior: ndarray | None = None

SO: float | None = None

SO_err: float | tuple[float, float] | None = None

Te_high: float | None = None

Te_high_err: float | tuple[float, float] | None = None

Te_low: float | None = None

Te_low_err: float | tuple[float, float] | None = None

alt_results: dict[str, AbundanceResult] | None = None

chi2: float | None = None

diagnostics: dict[str, str] | None = None

excluded_lines: list[str] | None = None

failures: dict[str, str] | None = None

icf_method: str | None = None

icf_values: dict[str, dict] | None = None

ionic: dict[str, float] | None = None

ionic_ul_details: dict[str, dict] | None = None

ionic_upper_limits: dict[str, float] | None = None

logU: float | None = None

logU_err: float | tuple[float, float] | None = None

lya_f_esc: float | None = None

lya_f_esc_details: dict | None = None

lya_f_esc_err: float | tuple[float, float] | None = None

lya_f_esc_posterior: ndarray | None = None

method: str

ne: float | None = None

ne_high: float | None = None

ne_low: float | None = None

ne_mid: float | None = None

ratios_used: list[str] | None = None

summary()

Return a human-readable summary string.

Returns:

Multi-line summary of the abundance measurement, including a diagnostics section explaining how each physical quantity was derived.

Return type:

str

jwspecabund.Te_low_from_high(Te_high, relation='desi')

Derive T_e(low) from T_e(high) using an empirical T_e-T_e relation.

Parameters:

• Te_high (float) – T_e(O++) in K.

• relation (str) – "desi" (default) — DESI DR2 (arXiv:2601.02463): T_low = 0.648 * T_high + 3270 "classical" — Garnett (1992): T_low = 0.7 * T_high + 3000

Returns:

T_e(low) in K.

Return type:

float

jwspecabund.cardelli_extinction(wave_A, Av, Rv=3.1)

Cardelli, Clayton & Mathis (1989) MW extinction curve.

Parameters:

• wave_A (np.ndarray) – Rest-frame wavelength in Angstroms.

• Av (float) – V-band extinction A_V.

• Rv (float) – Total-to-selective ratio (default 3.1).

Returns:

A(lambda) in magnitudes at each wavelength.

Return type:

np.ndarray

jwspecabund.compute_Av_balmer_pair(fluxes, errors, pair, law='salim', **kwargs)

Derive A_V from a single, explicitly-chosen Balmer decrement.

Like compute_Av_multi_balmer() but restricted to exactly one (numerator, denominator) pair — useful when you want to force the dust correction onto a high-SNR pair (e.g. ("Ha", "HBETA")) and ignore lower-SNR Balmer lines.

Parameters:

• fluxes (dict) – Observed (not dust-corrected) line fluxes and 1σ errors.

• errors (dict) – Observed (not dust-corrected) line fluxes and 1σ errors.

• pair (tuple of str) – (numerator, denominator) line names from the Balmer ladder ("Ha", "HBETA", "HGAMMA", "HDELTA", "H9", "H10"), e.g. ("Ha", "HBETA") for Hα/Hβ.

• law (str) – "salim" (default) or "cardelli".

• **kwargs – Extra arguments to the attenuation law (Rv, delta, B).

Returns:

Same shape as compute_Av_multi_balmer(): "Av", "Av_err", "individual", "n_lines" (1 if usable, else 0), "anchor" (the denominator line).

Return type:

dict

Raises:

ValueError – If either line name is not in the Balmer ladder.

jwspecabund.compute_Av_from_balmer(flux_num, flux_den, flux_num_err, flux_den_err, law='salim', intrinsic_ratio=0.468, wave_num_A=4341.68, wave_den_A=4862.68, **kwargs)

Derive A_V from a Balmer decrement.

By default uses Hgamma/Hbeta, but any pair of Balmer lines can be used by specifying the wavelengths and intrinsic ratio.

Parameters:

• flux_num (float) – Observed flux of the numerator line (e.g. Hgamma).

• flux_den (float) – Observed flux of the denominator line (e.g. Hbeta).

• flux_num_err (float) – Error on numerator flux.

• flux_den_err (float) – Error on denominator flux.

• law (str) – "salim" (default) or "cardelli".

• intrinsic_ratio (float) – Intrinsic flux ratio (default 0.468 for Hgamma/Hbeta).

• wave_num_A (float) – Rest wavelength of numerator line in Angstroms (default Hgamma).

• wave_den_A (float) – Rest wavelength of denominator line in Angstroms (default Hbeta).

• **kwargs – Extra arguments to the attenuation law.

Returns:

(Av, Av_err).

Return type:

tuple of float

jwspecabund.compute_Av_multi_balmer(fluxes, errors, law='salim', snr_min=3.0, anchor='HBETA', **kwargs)

Derive A_V from every available Balmer decrement.

Ratios each detected Balmer line bluer than the chosen anchor against it, and returns individual and weighted-average A_V values. Only lines at shorter wavelength than the anchor are used: anchor Hβ therefore uses Hγ/Hβ, Hδ/Hβ, H9/Hβ, H10/Hβ (Hα is excluded), while anchor Hα uses Hβ/Hα, Hγ/Hα, Hδ/Hα, H9/Hα, H10/Hα (all of them).

Parameters:

• fluxes (dict) – Emission-line fluxes (observed, not dust-corrected).

• errors (dict) – Corresponding 1σ errors.

• law (str) – "salim" (default) or "cardelli".

• snr_min (float) – Minimum SNR for a Balmer line to be included (default 3.0).

• anchor (str) – Reference (denominator) line for the decrement. "HBETA" (default) uses only bluer lines (Hγ/Hβ, Hδ/Hβ, H9/Hβ, H10/Hβ), excluding Hα. "Ha" uses every other Balmer line (Hβ/Hα, Hγ/Hα, Hδ/Hα, H9/Hα, H10/Hα).

• **kwargs – Extra arguments to the attenuation law (Rv, delta, B).

Returns:

Keys: "Av" (weighted mean), "Av_err" (error on mean), "individual" (list of dicts with per-line results), "n_lines" (number of lines used), "anchor" (anchor name).

Return type:

dict

jwspecabund.compute_NO_martinez25(ionic, logU, Z_Zsun, ne=100.0)

Compute total N/O using the best available Martinez+2025 ICF.

Selects the ICF based on which ionic abundances are present, with the following priority (most to least preferred):

1. ICF 5: (N2+ + N3+)/O2+ — pure UV, best at high logU

2. ICF 4: (N+ + N2+)/(O+ + O2+) — mixed UV+optical

3. ICF 2: N2+/O2+ — UV only, single ion

4. ICF 1: N+/O+ — classical optical

5. ICF 3: N3+/O2+ — large corrections, least preferred

Parameters:

• ionic (dict) – Ionic abundance dict, e.g. {"O+/H+": val, "O++/H+": val, "N+/H+": val, "N++/H+": val, "N+++/H+": val}.

• logU (float) – Ionisation parameter log(U).

• Z_Zsun (float) – Gas-phase metallicity in solar units.

• ne (float) – Electron density in cm^-3 (default 100).

Returns:

(N/O, icf_name) where icf_name identifies the ICF used (e.g. "NppNppp_Opp"). Returns (None, None) if no suitable ionic abundances are available.

Return type:

tuple

jwspecabund.compute_Te_NII(flux_5756, flux_6585, ne)

Compute T_e(N+) from the [NII] auroral/nebular ratio.

Parameters:

• flux_5756 (float) – [NII] 5756 flux.

• flux_6585 (float) – [NII] 6585 flux.

• ne (float) – Electron density in cm^-3.

Returns:

T_e(N+) in K.

Return type:

float

jwspecabund.compute_Te_OIII(flux_4363, flux_5007, flux_4959, ne)

Compute T_e(O++) from the [OIII] auroral/nebular ratio.

Uses PyNEB getTemDen() on the O++ atom with the standard diagnostic ratio [OIII] 4363 / ([OIII] 5007 + [OIII] 4959).

Parameters:

• flux_4363 (float) – [OIII] 4363 flux.

• flux_5007 (float) – [OIII] 5007 flux.

• flux_4959 (float) – [OIII] 4959 flux.

• ne (float) – Electron density in cm^-3.

Returns:

T_e(O++) in K.

Return type:

float

jwspecabund.compute_Te_OIII_1666(flux_1666, flux_5007, flux_4959, ne)

Compute T_e(O++) from the [OIII] UV/optical ratio 1666/(5007+4959).

Uses the O III] 1666 Å intercombination line (5→2 transition) as a UV auroral diagnostic when [OIII] 4363 is unavailable or low-SNR. The emissivity ratio 1666/(5007+4959) is monotonically increasing with T_e and more temperature-sensitive than 4363/(5007+4959) due to the larger energy gap (7.5 eV vs 2.8 eV).

Parameters:

• flux_1666 (float) – O III] 1666 flux (dust-corrected).

• flux_5007 (float) – [OIII] 5007 flux (dust-corrected).

• flux_4959 (float) – [OIII] 4959 flux (dust-corrected).

• ne (float) – Electron density in cm^-3.

Returns:

T_e(O++) in K.

Return type:

float

jwspecabund.compute_abundances(result, z, *, dust_correct=True, dust_law='salim', Av=None, Av_err=None, Av_prior='gaussian', method='auto', snr_auroral=3.0, snr_line=2.0, ne_high_max=500000.0, snr_ne=3.0, snr_logU=1.5, n_mc=1000, Te_relation='desi', Rv=3.15, delta=-0.35, B_bump=2.27, icf_method='auto', snr_NO=1.5, icf_tier=None, ne_low_override=None, ne_mid_override=None, ne_high_override=None, snr_balmer=3.0, balmer_anchor='HBETA', balmer_pair=None, forward_sampler='emcee', forward_n_walkers=32, forward_n_steps=5000, forward_n_burn=1000, forward_n_live=500, forward_seed=42, forward_progress=True, progress=True, n_posterior=1000)

Compute chemical abundances from a fitting result.

Dust handling (default)

By default A_V is computed once from the weighted-mean Balmer decrement of the median posterior fluxes and applied as a deterministic correction to every posterior draw — i.e. A_V is held constant across draws. The Balmer-propagated uncertainty on A_V is reported on the returned AbundanceResult (Av, Av_err) but is not marginalised into the abundance posteriors.

To marginalise over A_V instead (each draw dust-corrected with its own A_V sample), opt in by passing Av_err > 0 explicitly — see Av_err below.

type result:

Any

param result:

Emission-line fitting result from jwspecfit or jwspecmcmc.

type result:

FitResult | BroadFitResult | MCMCResult | MCMCBroadFitResult

type z:

float

param z:

Source redshift.

type z:

float

type dust_correct:

bool

param dust_correct:

Whether to apply dust correction (default True).

type dust_correct:

bool

type dust_law:

str

param dust_law:

"salim" (default) or "cardelli".

type dust_law:

str

type Av:

float | None

param Av:

V-band attenuation. If None, derived from Balmer decrement (see balmer_anchor).

type Av:

float or None

type balmer_anchor:

str

param balmer_anchor:

Denominator line for the multi-Balmer A_V derivation. Only Balmer lines bluer than the anchor (and above snr_balmer) are ratioed against it. "HBETA" (default) therefore uses Hγ/Hβ, Hδ/Hβ, H9/Hβ, H10/Hβ — Hα is excluded. "Ha" uses every other Balmer line (Hβ/Hα, Hγ/Hα, Hδ/Hα, H9/Hα, H10/Hα). So anchor on Hα when you want Hα included, on Hβ when you don’t. Ignored when Av is supplied directly, or when balmer_pair is set.

type balmer_anchor:

str

type balmer_pair:

tuple[str, str] | None

param balmer_pair:

Force the A_V derivation onto a single Balmer decrement instead of the multi-line fit, given as (numerator, denominator) line names from the ladder ("Ha", "HBETA", "HGAMMA", "HDELTA", "H9", "H10"). For example balmer_pair=("Ha", "HBETA") uses only Hα/Hβ — useful to avoid low-SNR Balmer lines. Overrides balmer_anchor / snr_balmer for the A_V step. Ignored when Av is supplied directly. None (default) keeps the multi-Balmer fit.

type balmer_pair:

tuple of str or None

type Av_err:

float | None

param Av_err:

Controls whether A_V is marginalised over. None (default) keeps A_V fixed at its central value — derived once from the median Balmer decrement, or taken from Av if supplied — and applies a single deterministic dust correction to all draws. Passing a positive value here switches on per-draw resampling: each posterior draw is dust-corrected with an A_V drawn from the chosen Av_prior (centred on Av, or on the Balmer-derived value when Av is None). The Balmer-derived error is always reported on the result; it is used for resampling only when you pass it back in explicitly here.

type Av_err:

float or None

type Av_prior:

str

param Av_prior:

Prior shape for A_V sampling: "gaussian" (default) draws from N(Av, Av_err) clipped at 0; "uniform" draws from U(max(Av − Av_err, 0), Av + Av_err). Only matters when Av_err is set.

type Av_prior:

str

type method:

str

param method:

"auto" (default), "direct", "forward", or "strong_line". "auto" uses direct if [OIII] 4363 SNR >= snr_auroral. "forward" runs the Bayesian forward model (Cullen+25) — see forward_model().

type method:

str

type snr_auroral:

float

param snr_auroral:

Minimum SNR for [OIII] 4363 to use the direct method (default 3.0).

type snr_auroral:

float

type snr_line:

float

param snr_line:

Minimum per-line SNR for inclusion in the abundance calculation (default 2.0). Lines below this threshold are removed from the flux dict after dust correction. Does not affect the auroral-line SNR check (controlled by snr_auroral) or lines essential for T_e computation (OIII 4363/5007/4959, Hbeta). Set to 0 to disable filtering.

type snr_line:

float

type ne_high_max:

float

param ne_high_max:

Maximum allowed high-ionisation electron density in cm^-3 (default 5e5). If n_e(high) from a UV doublet exceeds this, the code falls back to n_e(low). Prevents unphysical density estimates from noisy doublet ratios.

type ne_high_max:

float

type snr_ne:

float

param snr_ne:

Minimum SNR for both members of a density-sensitive doublet (default 3.0). Doublets where either member has flux / error < snr_ne are skipped, and the default density (300 cm^-3) is used. Set to 0 to disable.

type snr_ne:

float

type snr_logU:

float

param snr_logU:

Minimum total-doublet SNR for NIV] and NIII] when computing log(U) from N43 (default 1.5). The summed doublet flux is divided by the quadrature-summed error; if this is below the threshold, N43 is skipped and O32 is used instead.

type snr_logU:

float

type n_mc:

int

param n_mc:

Monte Carlo iterations for error propagation (default 1000).

type n_mc:

int

type Te_relation:

str

param Te_relation:

T_e-T_e relation: "desi" (default) or "classical".

type Te_relation:

str

type Rv:

float

param Rv:

Total-to-selective ratio for Salim law (default 3.15).

type Rv:

float

type delta:

float

param delta:

Slope deviation for Salim law (default -0.35).

type delta:

float

type B_bump:

float

param B_bump:

UV bump strength for Salim law (default 2.27).

type B_bump:

float

type icf_method:

str

param icf_method:

ICF scheme for the direct method. "auto" (default): use Martinez+25 when logU is available, fall back to Izotov+06 otherwise. "izotov06": always use Izotov+06 ICFs (N/O = ICF × N⁺/O⁺, independent of logU). "martinez25": force Martinez+25 ICFs (requires logU). "direct_sum": sum all detected nitrogen ions directly (Topping+2024, Yanagisawa+2025, Cameron+2023). No ICF or logU needed; falls back to Izotov+06 if only N⁺ is available.

type icf_method:

str

type snr_NO:

float

param snr_NO:

Minimum total-line SNR for each nitrogen ion when using icf_method="direct_sum" (default 2.0). For doublets (NIII], NIV], NV), the summed flux is divided by the quadrature-summed error. Ions below this threshold are excluded from the direct sum, causing the code to fall through to a lower tier (or Izotov+06).

type snr_NO:

float

type ne_low_override:

float | None

param ne_low_override:

If set, use this value (cm^-3) for the low-ionisation zone density instead of deriving it from [SII] or [OII].

type ne_low_override:

float or None

type ne_mid_override:

float | None

param ne_mid_override:

If set, use this value (cm^-3) for the mid-ionisation zone density instead of deriving it from CIII].

type ne_mid_override:

float or None

type ne_high_override:

float | None

param ne_high_override:

If set, use this value (cm^-3) for the high-ionisation zone density instead of deriving it from NIV] (or the fallback chain). Useful when the CIII]-derived fallback is suspect.

type ne_high_override:

float or None

type forward_sampler:

str

param forward_sampler:

Sampler for forward model: "emcee" or "dynesty" (default "emcee").

type forward_sampler:

str

type forward_n_walkers:

int

param forward_n_walkers:

Number of walkers for emcee forward model (default 32).

type forward_n_walkers:

int

type forward_n_steps:

int

param forward_n_steps:

MCMC steps for emcee forward model (default 5000).

type forward_n_steps:

int

type forward_n_burn:

int

param forward_n_burn:

Burn-in steps for emcee forward model (default 1000).

type forward_n_burn:

int

type forward_n_live:

int

param forward_n_live:

Live points for dynesty forward model (default 500).

type forward_n_live:

int

type forward_seed:

int

param forward_seed:

Random seed for the forward model (default 42).

type forward_seed:

int

type forward_progress:

bool

param forward_progress:

Show progress bar for the forward model (default True). Deprecated — use progress instead.

type forward_progress:

bool

type progress:

bool

param progress:

Show progress bars for MC / posterior loops (default True).

type progress:

bool

type n_posterior:

int

param n_posterior:

Maximum number of posterior samples to propagate through PyNEB / strong-line calculations (default 1000). If the MCMC posterior has more samples, a random subsample is drawn.

type n_posterior:

int

returns:

Chemical abundance measurement.

rtype:

AbundanceResult

Parameters:

• result (Any)

• z (float)

• dust_correct (bool)

• dust_law (str)

• Av (float | None)

• Av_err (float | None)

• Av_prior (str)

• method (str)

• snr_auroral (float)

• snr_line (float)

• ne_high_max (float)

• snr_ne (float)

• snr_logU (float)

• n_mc (int)

• Te_relation (str)

• Rv (float)

• delta (float)

• B_bump (float)

• icf_method (str)

• snr_NO (float)

• icf_tier (str | None)

• ne_low_override (float | None)

• ne_mid_override (float | None)

• ne_high_override (float | None)

• snr_balmer (float)

• balmer_anchor (str)

• balmer_pair (tuple[str, str] | None)

• forward_sampler (str)

• forward_n_walkers (int)

• forward_n_steps (int)

• forward_n_burn (int)

• forward_n_live (int)

• forward_seed (int)

• forward_progress (bool)

• progress (bool)

• n_posterior (int)

Return type:

AbundanceResult

jwspecabund.compute_ionic_abundances(fluxes, Te_high, Te_low, ne, ne_mid=None, ne_high=None)

Compute all available ionic abundances.

Parameters:

• fluxes (dict) – Dust-corrected emission-line fluxes keyed by line name. Must include "HBETA" for normalisation.

• Te_high (float) – T_e(O++) in K.

• Te_low (float) – T_e(O+/N+) in K.

• ne (float) – Electron density in cm^-3 (low-ionisation zone).

• ne_mid (float, optional) – Electron density for the intermediate-ionisation zone (cm^-3). Traced by CIII] 1907/1909 (~24 eV). If None, defaults to ne. Used for O²⁺, Ne²⁺, C²⁺, N²⁺, S²⁺, Ar²⁺. O²⁺ and Ne²⁺ use this intermediate-zone density (not ne_high): the [OIII] 5007/Hβ and [NeIII] 3869/Hβ abundances are density- insensitive below ~10⁴–10⁵ cm⁻³, and CIII] (24–48 eV) overlaps the O²⁺ zone (35–55 eV), whereas NIV] (47–77 eV) traces more highly-ionised gas. This decouples O²⁺/Ne²⁺ from the noisy high-ionisation N IV] density.

• ne_high (float, optional) – Electron density for the high-ionisation zone (cm^-3). Traced by NIV] 1483/1486 (~47 eV). If None, defaults to ne_mid. Used for N³⁺, N⁴⁺, C³⁺.

Returns:

Ionic abundances, e.g. {"O+/H+": val, "O++/H+": val, ...}.

Return type:

dict

jwspecabund.compute_line_ratios(line_fluxes, line_errors, snr_thresh=1.5)

Compute available strong-line diagnostic ratios.

Parameters:

• line_fluxes (dict) – {line_name: flux}.

• line_errors (dict) – {line_name: flux_err}.

• snr_thresh (float) – Minimum SNR for a line to be used (default 1.5).

Returns:

{ratio_name: {"val": log10_ratio, "err": error}}.

Return type:

dict

jwspecabund.compute_lya_escape_fraction(lya_flux, lya_flux_err, fluxes_corr, errors_corr, snr_min=3.0)

Compute Lyα escape fraction from dust-corrected Balmer lines.

Each available Balmer line predicts an intrinsic Lyα flux via Case B recombination (T=10⁴ K, n_e=100 cm⁻³). The escape fraction is the ratio of observed Lyα to intrinsic Lyα.

When multiple Balmer lines are available, f_esc is computed as the inverse-variance weighted mean of the individual estimates.

The observed Lyα flux should NOT be dust-corrected — the escape fraction encapsulates both dust attenuation and resonant scattering.

Parameters:

• lya_flux (float) – Observed (not dust-corrected) Lyα flux.

• lya_flux_err (float) – 1σ error on the observed Lyα flux.

• fluxes_corr (dict) – Dust-corrected emission-line fluxes keyed by line name.

• errors_corr (dict) – Dust-corrected 1σ errors keyed by line name.

• snr_min (float) – Minimum SNR for a Balmer line to be included (default 3.0).

Returns:

Keys: "f_esc" (weighted mean), "f_esc_err", "individual" (list of per-line dicts with line, f_esc, f_esc_err, lya_intrinsic, lya_intrinsic_err), "n_lines" (number used).

Return type:

dict

jwspecabund.compute_lya_escape_fraction_mc(lya_flux, lya_flux_err, fluxes_corr, errors_corr, Av, Av_err, Av_prior='gaussian', dust_law='salim', snr_min=3.0, n_mc=1000, rng=None, **dust_kwargs)

Compute Lyα escape fraction with MC error propagation.

Draws n_mc samples from the observed Lyα flux, dust-corrected Balmer fluxes, and A_V, recomputing the dust correction at each draw. Returns the median and 16th/84th percentile uncertainties.

Parameters:

• lya_flux (float) – Observed (not dust-corrected) Lyα flux.

• lya_flux_err (float) – 1σ error on Lyα flux.

• fluxes_corr (dict) – Dust-corrected Balmer fluxes (used as central values; the MC re-applies dust correction from observed fluxes internally).

• errors_corr (dict) – Dust-corrected 1σ errors.

• Av (float) – Central A_V value.

• Av_err (float) – 1σ error on A_V (Gaussian) or half-width (uniform).

• Av_prior (str) – Prior shape for A_V: "gaussian" (default, draws from N(Av, Av_err) clipped at 0) or "uniform" (draws from U(max(Av − Av_err, 0), Av + Av_err)).

• dust_law (str) – "salim" or "cardelli".

• snr_min (float) – Minimum SNR for a Balmer line to be included (default 3.0).

• n_mc (int) – Number of Monte Carlo draws (default 1000).

• rng (np.random.Generator or None) – Random number generator (for reproducibility).

• **dust_kwargs – Extra arguments to the dust law (Rv, delta, B).

Returns:

Keys: "f_esc" (median), "f_esc_err" (lo, hi) as 16th/84th percentile half-widths, "f_esc_posterior" (1D array of MC samples), "individual" (per-line results), "n_lines".

Return type:

dict

jwspecabund.compute_ne(flux_line1, flux_line2, doublet='SII', Te_guess=10000.0)

Compute electron density from a density-sensitive doublet.

Parameters:

• flux_line1 (float) – Flux of the blue member (e.g. [SII] 6718 or [OII] 3726).

• flux_line2 (float) – Flux of the red member (e.g. [SII] 6732 or [OII] 3729).

• doublet (str) – "SII" (default) or "OII".

• Te_guess (float) – Electron temperature guess in K (default 10^4).

Returns:

Electron density n_e in cm^-3.

Return type:

float

jwspecabund.compute_ne_CIII(flux_1907, flux_1909, Te_guess=10000.0)

Compute electron density from the CIII] 1907/1909 ratio.

Probes the intermediate-ionisation zone.

Parameters:

• flux_1907 (float) – CIII] 1907 flux.

• flux_1909 (float) – CIII] 1909 flux.

• Te_guess (float) – Electron temperature guess in K (default 10^4).

Returns:

Electron density n_e in cm^-3.

Return type:

float

jwspecabund.compute_ne_NIV(flux_1483, flux_1486, Te_guess=10000.0)

Compute electron density from the NIV] 1483/1486 ratio.

Probes the high-ionisation zone.

Parameters:

• flux_1483 (float) – NIV] 1483 flux.

• flux_1486 (float) – NIV] 1486 flux.

• Te_guess (float) – Electron temperature guess in K (default 10^4).

Returns:

Electron density n_e in cm^-3.

Return type:

float

jwspecabund.compute_total_abundances(ionic, logU=None, Z_Zsun=None, ne=None, icf_method='auto', ionic_upper_limits=None, _lock_NO_icf=None)

Derive total element abundances from ionic abundances + ICFs.

Parameters:

• ionic (dict) – Ionic abundance dict from compute_ionic_abundances().

• logU (float, optional) – Ionisation parameter log(U). Required for Martinez+25 ICFs.

• Z_Zsun (float, optional) – Gas-phase metallicity in solar units. Required for Martinez+25 ICFs.

• ne (float, optional) – Electron density in cm^-3 for Martinez+25 ICF density interpolation.

• icf_method (str) – "auto" (default): use Martinez+25 for N/O when logU is provided, fall back to Izotov+06 otherwise. "martinez25": force Martinez+25 ICFs (requires logU and Z_Zsun). "izotov06": use Izotov+06 ICFs only. "direct_sum": sum all detected nitrogen ions directly (Topping+2024, Yanagisawa+2025, Cameron+2023). Tiered fallback: Tier 1 (N⁺ + N²⁺ + N³⁺) / (O⁺ + O²⁺), Tier 2/3 (N²⁺ + N³⁺) / O²⁺, Tier 4 Izotov+06 optical fallback.

• ionic_upper_limits (dict[str, float] | None)

• _lock_NO_icf (str | None)

Returns:

Total abundance ratios: "O/H", "N/O", "S/O", "Ne/O", "Ar/O", "C/O", "N/O_UV" as available. When Martinez+25 is used, also includes "NO_icf_name" and "icf_method" keys.

Return type:

dict

jwspecabund.dust_correct_fluxes(line_fluxes, Av, law='salim', **kwargs)

Apply dust correction to a set of emission-line fluxes.

Parameters:

• line_fluxes (dict) – {line_name: (flux, flux_err, rest_wave_A)}.

• Av (float) – V-band attenuation.

• law (str) – "salim" (default) or "cardelli".

• **kwargs – Extra keyword arguments passed to the attenuation law (e.g. Rv, delta, B for Salim).

Returns:

{line_name: (corrected_flux, corrected_flux_err)}.

Return type:

dict

jwspecabund.forward_model(line_fluxes, line_errors, *, sampler='emcee', n_walkers=32, n_steps=5000, n_burn=1000, n_live=500, seed=42, progress=True)

Run a Bayesian forward model on emission-line fluxes.

Free parameters are log(T_e), log(n_e), and one log(ionic abundance) per detected ion. The model predicts line/Hbeta flux ratios using PyNEB CEL emissivities and the Aller (1984) Hbeta formula, and the Gaussian likelihood compares predictions to observations.

Parameters:

• line_fluxes (dict) – {line_name: flux}. Must include "HBETA".

• line_errors (dict) – {line_name: flux_err}.

• sampler (str) – "emcee" (default) or "dynesty".

• n_walkers (int) – Number of walkers for emcee (default 32).

• n_steps (int) – Total MCMC steps for emcee (default 5000).

• n_burn (int) – Burn-in steps to discard for emcee (default 1000).

• n_live (int) – Number of live points for dynesty (default 500).

• seed (int) – Random seed (default 42).

• progress (bool) – Show progress bar (default True).

Returns:

Result dictionary with keys:

• param_names, param_labels — parameter metadata

• samples — raw posterior samples, shape (n, n_dim)

• Te, Te_err, Te_posterior

• ne, ne_err, ne_posterior

• OH, OH_err, OH_posterior

• NO, NO_err, NO_posterior (if nitrogen detected)

• ionic — dict of median log-ionic abundances

• ionic_posteriors — dict of log-ionic posterior arrays

• ratio_names — observable labels

Return type:

dict

Examples

>>> import jwspecabund
>>> fluxes = {
...     "HBETA": 1.0, "OIII_5007": 5.0,
...     "OIII_4363": 0.05, "NeIII_3869": 0.3,
... }
>>> errors = {k: 0.05 * v for k, v in fluxes.items()}
>>> res = jwspecabund.forward_model(fluxes, errors, n_steps=2000)
>>> print(f"12+log(O/H) = {res['OH']:.2f}")

jwspecabund.hbeta_emissivity_aller84(Te)

Hbeta volume emissivity coefficient using the Aller (1984) Case B formula.

\[\varepsilon_{H\beta} = \alpha_{H\beta}^{\rm eff}(T)\;h\nu_{H\beta}\]

where \(\alpha_{H\beta}^{\rm eff} = 3.03\times10^{-14}\,(T/10^4)^{-0.874}\) cm³ s^-1.

Unlike PyNEB’s recombination tables, this formula is valid at all temperatures — no upper-bound limitation at 30,000 K.

Parameters:

Te (float) – Electron temperature in K.

Returns:

\(\varepsilon_{H\beta}\) in erg cm³ s^-1.

Return type:

float

jwspecabund.icf_argon(O_plus, O_total)

ICF for argon (Izotov+06 eqs. 22/23).

Ar/O = ICF_Ar * Ar++/O++.

Parameters:

• O_plus (float) – Ionic abundance O+/H+.

• O_total (float) – Total oxygen abundance O/H.

Returns:

ICF_Ar.

Return type:

float

jwspecabund.icf_neon(O_plus, O_total)

ICF for neon (Izotov+06 eq. 19).

Ne/O = ICF_Ne * Ne++/O++.

Parameters:

• O_plus (float) – Ionic abundance O+/H+.

• O_total (float) – Total oxygen abundance O/H.

Returns:

ICF_Ne.

Return type:

float

jwspecabund.icf_nitrogen(O_plus, O_total)

ICF for nitrogen (Izotov+06 eq. 18).

N/O = ICF_N * N+/O+, where ICF_N accounts for the contribution of higher ionisation states of nitrogen.

Parameters:

• O_plus (float) – Ionic abundance O+/H+.

• O_total (float) – Total oxygen abundance O/H (= O+/H+ + O++/H+).

Returns:

ICF_N. Returns 1.0 for the standard assumption N/O ~ N+/O+.

Return type:

float

jwspecabund.icf_sulfur(O_plus, O_total)

ICF for sulfur (Izotov+06 eq. 20).

S/O = ICF_S * (S+ + S++)/O.

Parameters:

• O_plus (float) – Ionic abundance O+/H+.

• O_total (float) – Total oxygen abundance O/H.

Returns:

ICF_S.

Return type:

float

jwspecabund.log_U_from_N43(log_N43, Z_Zsun, ne=100.0)

Compute log(U_high) from the N43 diagnostic.

Parameters:

• log_N43 (float) – log10(NIV] 1486 / NIII] 1750).

• Z_Zsun (float) – Gas-phase metallicity in solar units.

• ne (float) – Electron density in cm^-3 (default 100).

Returns:

log(U_high).

Return type:

float

Notes

Valid for 0.05 <= Z/Zsun <= 0.4, -3.0 <= log(N43) <= 0.5. N43 is density-insensitive — recommended over O32 when NIV] and NIII] are available (Martinez+2025 Section 4.2). Martinez+2025 Table 3 (lower).

jwspecabund.log_U_from_O32(log_O32, Z_Zsun, ne=100.0)

Compute log(U_int) from the O32 diagnostic.

Parameters:

• log_O32 (float) – log10([OIII] 5007 / [OII] 3727).

• Z_Zsun (float) – Gas-phase metallicity in solar units.

• ne (float) – Electron density in cm^-3 (default 100).

Returns:

log(U_int).

Return type:

float

Notes

Valid for 0.05 <= Z/Zsun <= 0.4, -1.0 <= log(O32) <= 2.5. O32 is density-sensitive; prefer N43 when available. Martinez+2025 Table 3 (upper).

jwspecabund.salim_attenuation(wave_A, Av, Rv=3.15, delta=-0.35, B=2.27)

Salim+18/Noll+09 modified Calzetti attenuation curve.

A(lambda) = E(B-V) * [k’_Calzetti * (lambda/0.55)^delta + D_bump * B]

Parameters:

• wave_A (np.ndarray) – Rest-frame wavelength in Angstroms.

• Av (float) – V-band attenuation A_V.

• Rv (float) – Total-to-selective ratio (default 3.15).

• delta (float) – Power-law slope deviation (default -0.35).

• B (float) – UV bump strength (default 2.27).

Returns:

A(lambda) in magnitudes at each wavelength.

Return type:

np.ndarray

jwspecabund.sanders25_metallicity(line_fluxes, line_errors, n_mc=1000, snr_thresh=1.5, seed=42, progress=True)

Derive 12+log(O/H) via simultaneous Sanders+25 calibrations.

Parameters:

• line_fluxes (dict) – {line_name: flux}.

• line_errors (dict) – {line_name: flux_err}.

• n_mc (int) – Monte Carlo iterations for error propagation (default 1000).

• snr_thresh (float) – Minimum SNR for a ratio to be used (default 1.5).

• seed (int) – Random seed for reproducibility.

• progress (bool) – Show a tqdm progress bar (default True).

Returns:

(Z_best, Z_lo, Z_hi, chi2, ratios_used, Z_mc_samples) where Z_best is 12+log(O/H), Z_lo/Z_hi are 16th/84th percentiles, chi2 is the best-fit chi-squared, ratios_used lists the diagnostics, and Z_mc_samples is the full MC array.

Return type:

tuple

Raises:

ValueError – If no valid diagnostic ratios are available.

Modules

direct	Direct T_e method for ionic and total abundances.
dust	Dust correction for emission-line fluxes.
forward	Bayesian forward model for emission-line abundance determination.
icf	Ionisation correction factors.
martinez25_icf	Density-dependent ICFs and logU diagnostics from Martinez+2025.
result	Abundance result containers.
strong_line	Strong-line metallicity calibrations from Sanders et al. (2025).