`jwspecfit` — least-squares fitting

jwspecfit is the least-squares engine: it fits resolution-aware, bin-averaged Gaussian profiles to all observable emission lines simultaneously, and returns a FitResult with fluxes, SNRs, centroids, widths, and equivalent widths.

Note

For science-quality uncertainties the authors recommend the Bayesian MCMC fitter, jwspecmcmc. Use this least-squares engine for quick looks, initial guesses, and fast BIC model selection.

Loading spectra

import jwspecfit

spec = jwspecfit.read_fits("spectrum.fits", z=6.0)

read_fits accepts both the JWST SPEC1D BinTable convention (default) and arbitrary 1-D FITS files: it auto-detects column names (wave/wavelength/lam/lambda/loglam, flux, err/sigma/ivar), reads units from the TUNIT keywords (μm/Å/nm, μJy/mJy/Jy, erg/s/cm²/Å) and falls through to image HDUs with WCS keywords (CRVAL1, CDELT1, CRPIX1, CTYPE1, CUNIT1, BUNIT). Override the auto-detection with hdu=, wave_col=, flux_col=, err_col= if needed.

See the quickstart guide for the full list of readers (read_fits, read_dict, read_npz).

Quick interactive view

To pop a spectrum open in plotly without first running a fit, call plot_spectrum_interactive directly on a path or a Spectrum:

fig = jwspecfit.plot_spectrum_interactive("spectrum.fits", z=6.0)
fig = jwspecfit.plot_spectrum_interactive("stack.npz", rest_frame=True)
fig = jwspecfit.plot_spectrum_interactive(spec, flux_unit="flam")
fig.show()

The plot has zoom/pan, hover read-outs, and a ±1σ band whenever an error array is available.

Fitting

result = jwspecfit.fit_lines(spec, z=6.0)

Common options:

result = jwspecfit.fit_lines(
    spec, z=6.0,
    lines=["OIII_5007", "OIII_4959", "HBETA"],   # restrict line list
    wave_range_A=(35000, 50000),                  # wavelength window
    wave_windows_A=[(13000, 17000), (35000, 52000)],  # multi-window
    n_boot=500,                                    # bootstrap iterations
    mode="auto",                                   # broad-Balmer mode
    niv_doublet_ratio=None,                        # fix NIV ratio
    ciii_doublet_ratio=None,                       # fix CIII] ratio
)

Inspecting results

Each LineResult contains:

Field	Description
`flux`	Integrated line flux (erg / s / cm²)
`flux_err`	1σ uncertainty from bootstrap
`snr`	Signal-to-noise ratio
`ew_A`	Rest-frame equivalent width (Å)
`amplitude`	Peak amplitude
`centroid_A`	Observed centroid wavelength (Å)
`sigma_A`	Gaussian sigma (Å, incl. instrumental broadening)

Upper limits for undetected lines:

ul = result.flux_upper_limit("OIII_4363", n_sigma=3.0)
uls = result.flux_upper_limits(n_sigma=3.0)

Broad Balmer components

Four nested models are tested and selected via BIC when mode="auto":

Narrow only.
Narrow + intermediate broad (FWHM ~ 500–2000 km s⁻¹).
Narrow + very broad (FWHM ~ 2000–5000 km s⁻¹).
Narrow + both broad components.

result = jwspecfit.fit_with_broad(spec, z=2.0, n_boot_bic=100)
print(result.selected_model)     # "narrow" | "broad1" | "broad2" | "both"

UV and absorption lines

Lines prefixed with abs_ are fit as negative Gaussians (low-ionisation ISM absorption, DLA flanks):

uv_lines = [
    "NV_1", "NV_2", "CIV_1", "CIV_2", "HEII_1640",
    "CIII]_1907", "CIII]",
    "abs_SiII1260", "abs_CII1334", "abs_SiIV1394", "abs_SiIV1403",
]
result = jwspecfit.fit_lines(spec, z=6.0, lines=uv_lines)

DLA fitting

jwspecfit.fit_NHI performs a joint Bayesian fit of the local DLA column density and the IGM damping wing, following the methodology of Pollock et al. (2026) Eq. 4:

\[F(\lambda) = \big[F_0 (\lambda/\lambda_\text{piv})^{\beta_\text{UV}} + \mathrm{Ly}\alpha(\lambda)\big] \;\exp(-\tau_\text{DLA})\;\exp(-\tau_\text{IGM})\]

The posterior over (log_F0, β_UV, log_NHI) (and optionally x_HI, plus a Lyα emission profile) is sampled with dynesty nested sampling. The Voigt profile is exact (Faddeeva); the IGM term uses the Miralda-Escudé (1998) closed form integrated from igm_z_min (default 5.3, Bosman+22) up to the source redshift.

Install the optional dependency first:

pip install jwspecfit[dla]   # installs dynesty

Standard call

import jwspecfit

spec = jwspecfit.read_npz("stack_zgt6_Muv19_21_prism.npz", z=0.0, R=400.0)

result = jwspecfit.fit_NHI(
    spec.wave_A, spec.flux_ujy, spec.err_ujy,
    z=0.0,                       # 0 for rest-frame stacks
    R=spec.R,                    # LSF convolution
    fit_x_HI=False,              # set True for individual z>6 spectra
    n_live=400, seed=42,
)
print(result.summary())
result.plot()

The defaults match Pollock+26: priors log_NHI ∈ [18, 24], x_HI ∈ [0, 1], β_UV ∈ [-4, 0], log_F0 auto-centred ±3 dex around the data; fit window 1216–3000 Å rest-frame; emission lines masked with mask_width_A=20 Å; Lyα emission optionally added inside the exp(−τ) factor.

Useful options

Argument	Purpose
`fit_x_HI`	Sample IGM neutral fraction jointly (Pollock+26 z > 6 mode)
`igm_z_min`	Lower-z bound of the IGM integration (default 5.3)
`b_kms`	Voigt Doppler parameter (default 30 km s⁻¹; insensitive at log N_HI > 20.3)
`prior_log_NHI`, `prior_x_HI`, `prior_beta_UV`, `prior_log_F0`	Override the uniform prior bounds
`mask_lines`, `mask_width_A`	Default `True`, 20 Å (PRISM-appropriate)
`fit_range_A`	Default (1216, 3000) Å rest-frame
`mask_regions_A`	Extra rest-frame regions to ignore
`Av`, `dust_law`	Pre-correct the spectrum for foreground attenuation
`fit_lya`, `lya_params`	Joint or fixed Lyα emission profile
`n_live`, `dlogz`, `seed`	dynesty controls

The returned DLAResult carries posterior samples for every free parameter (result.samples is a dict of arrays), median ± 16/84 % percentiles, log-evidence, and a 95th-percentile upper limit on log_NHI (result.log_NHI_upper95, result.is_upper_limit=True) for spectra where the wing is unconstrained — Pollock+26 do exactly this for 21 of their 48 sources.

D_Lyα equivalent-width statistic (Heintz+25)

When the N_HI–β degeneracy bites at low spectral resolution, the more robust diagnostic is the rest-frame equivalent width of the Lyα absorption integrated 1180–1350 Å (Heintz+25, Eq. 1). It is β-prior–independent to first order and directly comparable across the high-z DLA literature.

DLya = jwspecfit.compute_D_Lya(spec.wave_A, spec.flux_ujy, spec.err_ujy, z=0.0)
print(f"D_Lyα = {DLya['D_Lya']:.1f} ± {DLya['D_Lya_err']:.1f} Å")
# Interpretation (Heintz+25 §3.1):
#   D_Lyα < 0    → Lyα emitter
#   ~0 to 35 Å   → weak / no absorption
#   35 to 50 Å   → IGM-dominated regime
#   > 50 Å       → DLA-dominated (log N_HI ≳ 21.5)

Caveats at PRISM resolution

The N_HI–β posterior covariance is r ≈ −0.96 at PRISM R~100 (Keating+25), so log N_HI shifts by ~1 dex when the β prior is tightened from [-4, 0] to a star-forming-galaxy range like [-3, -1.5]. For a stack at z > 6 with PRISM data, report D_Lyα as the primary statistic and treat any specific log N_HI as a prior-conditioned posterior. Robust column densities require medium-resolution gratings (G140M / G235M, R ≳ 1000).

Plotting

fig = jwspecfit.plot_fit(result, save_path="fit.pdf")
fig = jwspecfit.plot_fit_interactive(result)   # plotly

See the Plotting & visualisation page for the full set of static and interactive helpers, line markers, and the 2-D + 1-D viewer.

Save / load

jwspecfit.save_result(result, "fit.npz")
loaded = jwspecfit.load_result("fit.npz")
jwspecfit.export_lines_txt(result, "lines.txt")

Why resolution-aware?

NIRSpec prism resolving power R varies from ~30 at 1 μm to ~300 at 5 μm, so narrow emission lines are often sub-pixel in the blue and super-pixel in the red. Evaluating a Gaussian at pixel centres biases the recovered flux systematically downward for under-sampled lines. jwspecfit computes the analytic integral of each Gaussian across each pixel’s bin edges using the error function, eliminating this bias.

jwspecfit — least-squares fitting