jwspecfit.dla

DLA column density fitter — measure N_HI from Lya damping wings.

Fits a UV power-law continuum attenuated by a damped Lyman-alpha absorber (DLA) to derive the neutral hydrogen column density N_HI in the galaxy’s ISM.

The model (following Pollock et al. 2026, Eq. 4) is:

F(lambda) = F0 * (lambda/lambda_pivot)^beta_UV * exp(-tau_DLA)

where tau_DLA = C * a * N_HI * H(a, x) uses the Voigt-Hjerting function H(a, x) evaluated via the Faddeeva function (exact).

The model is convolved with the instrumental line-spread function (Gaussian with FWHM = lambda / R) before comparison to data, as in Pollock et al. (2026).

Sampling is performed with dynesty nested sampling (matching the paper’s methodology).

References

Pollock, C. L., et al. 2026, A&A, arXiv:2602.11783.: Method and application to z > 9 galaxies.
Tepper-Garcia, T. 2006, MNRAS, 369, 2025.: Voigt-Hjerting function framework.

Functions

`compute_D_Lya`(wave_A, flux, flux_err, z, *)	Heintz+25 D_Lyα equivalent-width statistic (Eq.
`fit_NHI`(wave_A, flux, flux_err[, z, ...])	Joint Bayesian fit of the DLA + IGM damping wing (Pollock+26-style).
`tau_DLA`(wave_A, log_NHI[, z, b_kms])	Compute DLA optical depth at each wavelength.
`tau_IGM_DW`(wave_A, x_HI, z_source[, z_min, ...])	IGM damping-wing optical depth from a uniformly neutral IGM.
`voigt_H`(a, u)	Voigt-Hjerting function H(a, u) via the Faddeeva function.

Classes

DLAResult(log_NHI, log_NHI_err, beta_UV, ...)

Result of a DLA column density fit.

class jwspecfit.dla.DLAResult(log_NHI, log_NHI_err, beta_UV, beta_UV_err, log_F0, log_F0_err, Sigma_HI, samples, wave_fit, flux_fit, flux_err_fit, model_best, z=0.0, Av=0.0, log_evidence=0.0, _lya_subtracted=False, lya_params=None, x_HI=0.0, x_HI_err=(0.0, 0.0), fit_x_HI=False, igm_z_min=5.3, log_NHI_upper95=0.0, is_upper_limit=False, lam_pivot_A=None)[source]

Bases: object

Result of a DLA column density fit.

Parameters:

log_NHI (float)
log_NHI_err (tuple[float, float])
beta_UV (float)
beta_UV_err (tuple[float, float])
log_F0 (float)
log_F0_err (tuple[float, float])
Sigma_HI (float)
samples (dict[str, ndarray])
wave_fit (ndarray)
flux_fit (ndarray)
flux_err_fit (ndarray)
model_best (ndarray)
z (float)
Av (float)
log_evidence (float)
_lya_subtracted (bool)
lya_params (ndarray | None)
x_HI (float)
x_HI_err (tuple[float, float])
fit_x_HI (bool)
igm_z_min (float)
log_NHI_upper95 (float)
is_upper_limit (bool)
lam_pivot_A (float | None)

log_NHI

Median posterior log10(N_HI / cm^-2).

Type:: float

log_NHI_err

(lower, upper) 68% CI half-widths.

Type:: tuple of float

beta_UV

Median posterior UV spectral slope.

Type:: float

beta_UV_err

(lower, upper) 68% CI half-widths.

Type:: tuple of float

log_F0

Median posterior log10(F0) continuum normalisation.

Type:: float

log_F0_err

(lower, upper) 68% CI half-widths.

Type:: tuple of float

Sigma_HI

H I gas surface density in M_sun pc^-2 (from log_NHI).

Type:: float

samples

Full posterior samples: {“log_NHI”: arr, “beta_UV”: arr, “log_F0”: arr}.

Type:: dict

wave_fit

Wavelengths used in fit (after masking), observed frame.

Type:: np.ndarray

flux_fit

Dust-corrected fluxes used in fit.

Type:: np.ndarray

flux_err_fit

Dust-corrected flux errors used in fit.

Type:: np.ndarray

model_best

Best-fit model evaluated on wave_fit.

Type:: np.ndarray

z

Redshift used in fit.

Type:: float

Av

Dust correction applied.

Type:: float

log_evidence

Log-evidence (logZ) from dynesty.

Type:: float

Av: float = 0.0

Sigma_HI: float

beta_UV: float

beta_UV_err: tuple[float, float]

corner(**kwargs)[source]

Plot a corner plot of the posterior samples.

Requires the corner package.

Parameters:: **kwargs (Any) – Passed to corner.corner().
Returns:: The corner plot figure.
Return type:: matplotlib Figure

fit_x_HI: bool = False

flux_err_fit: ndarray

flux_fit: ndarray

igm_z_min: float = 5.3

is_upper_limit: bool = False

lam_pivot_A: float | None = None

log_F0: float

log_F0_err: tuple[float, float]

log_NHI: float

log_NHI_err: tuple[float, float]

log_NHI_upper95: float = 0.0

log_evidence: float = 0.0

lya_params: ndarray | None = None

model_best: ndarray

plot(ax=None, show_residuals=True, flux_unit='fnu', **kwargs)[source]

Plot the DLA fit over the data.

Parameters:

ax (matplotlib Axes, optional) – Axes to plot on. If None, creates a new figure. If show_residuals is True, this is ignored and a new figure with two panels is created.
show_residuals (bool) – If True, show a residual panel below the main plot.
flux_unit (str) – "fnu" for F_nu (default, same units as input) or "flam" for F_lambda (converted via F_lam = F_nu * c / lam^2).
**kwargs (Any) – Passed to the data plot (e.g. color, alpha).

Returns:

The figure object.

Return type:

matplotlib Figure

samples: dict[str, ndarray]

summary()[source]

Return a formatted summary string.

Return type:: str

wave_fit: ndarray

x_HI: float = 0.0

x_HI_err: tuple[float, float] = (0.0, 0.0)

z: float = 0.0

jwspecfit.dla.compute_D_Lya(wave_A, flux, flux_err, z, *, cont_range_A=(1250.0, 2600.0), int_range_A=(1180.0, 1350.0), mask_lines=True, mask_width_A=20.0, n_mc=1000, seed=42)[source]

Heintz+25 D_Lyα equivalent-width statistic (Eq. 1).

Fits a power-law continuum F_cont(λ) = F_1550 (λ/1550)^β over cont_range_A (lines masked) then integrates

D_Lyα = ∫ (1 − F_λ / F_cont(λ)) dλ / (1 + z)

over int_range_A (rest-frame). D_Lyα < 0 indicates a Lyα emitter; ≳50 Å indicates a strong damped absorber.

Parameters:

wave_A (np.ndarray) – Spectrum in observed-frame Å and arbitrary linear flux units.
flux (np.ndarray) – Spectrum in observed-frame Å and arbitrary linear flux units.
flux_err (np.ndarray) – Spectrum in observed-frame Å and arbitrary linear flux units.
z (float) – Source redshift.
cont_range_A (tuple) – Rest-frame range used to fit the power-law continuum (default 1250–2600 Å, matching Heintz+25).
int_range_A (tuple) – Rest-frame range over which the EW integral is taken (default 1180–1350 Å, matching Heintz+25 Eq. 1).
mask_lines (bool) – Whether to mask known UV emission lines from the continuum fit (default True). The integral is always taken over the full int_range_A — emission/absorption inside it is the signal.
mask_width_A (float) – Half-width of line masks for the continuum fit (rest-frame Å).
n_mc (int) – Monte-Carlo iterations for D_Lyα uncertainty propagation.
seed (int) – RNG seed.

Returns:

{"D_Lya": float, "D_Lya_err": float, "beta_UV": float, "beta_UV_err": float, "F_1550": float, "F_1550_err": float}

Return type:

dict

jwspecfit.dla.fit_NHI(wave_A, flux, flux_err, z=0.0, *, fit_x_HI=False, igm_z_min=5.3, b_kms=30.0, continuum=None, Av=0.0, dust_law='cardelli', Rv=3.1, R=None, mask_lines=True, mask_width_A=20.0, mask_lya_emission_width_A=8.0, fit_range_A=(1216.0, 3000.0), mask_regions_A=None, fit_lya=False, lya_params=None, emission_lines=None, lam_pivot_A=None, prior_log_NHI=(18.0, 24.0), prior_x_HI=(0.0, 1.0), prior_beta_UV=(-4.0, 0.0), prior_log_F0=None, prior_log_A_lya=(-22.0, -15.0), prior_sigma_lya=(1.0, 30.0), prior_alpha_lya=(0.0, 5.0), n_live=500, dlogz=0.5, seed=42)[source]

Joint Bayesian fit of the DLA + IGM damping wing (Pollock+26-style).

Implements the Pollock et al. (2026) joint nested-sampling fit over (log_F0, β_UV, log_NHI[, x_HI[, Lyα]]) with an exact Voigt-Hjerting profile and the Miralda-Escudé (1998) IGM damping wing in the Totani+06 closed form.

Parameters:

wave_A (np.ndarray) – Wavelength (observed-frame Å) and flux density arrays.
flux (np.ndarray) – Wavelength (observed-frame Å) and flux density arrays.
flux_err (np.ndarray) – Wavelength (observed-frame Å) and flux density arrays.
z (float) – Source redshift (use 0 for rest-frame stacks).
fit_x_HI (bool) – If True, sample over IGM neutral fraction x_HI in addition to the local DLA. Recommended at z > 6.
igm_z_min (float) – Lower-redshift bound of the IGM integration (default 5.3, Bosman+22 reionisation end).
b_kms (float) – Doppler parameter (km/s) for the Voigt profile. The damping wing is insensitive to b for log N_HI > 20.3.
continuum (np.ndarray, optional) – Pre-smoothed continuum estimate. When provided the model is fit against this smoothed array instead of the raw flux — emission lines are then implicitly removed.
Av (float, str, float) – Optional foreground dust correction applied before fitting.
dust_law (float, str, float) – Optional foreground dust correction applied before fitting.
Rv (float, str, float) – Optional foreground dust correction applied before fitting.
R (float or None) – Spectral resolving power for LSF convolution.
mask_lines (bool) – If True (default), mask known UV emission lines. Lyα, the NV λλ1238.8, 1242.8 doublet, and any abs_* low-ionisation absorption features are excluded from this mask: Lyα and NV sit inside the damping wing (1216–1263 Å rest) where the N_HI constraint lives, and any DLA strong enough to fit absorbs NV at τ ≫ 1 (no emission left to contaminate); abs_* features are foreground-gas signatures aligned with the DLA.
mask_width_A (float) – Half-width of the general emission-line mask (rest-frame Å). Default 20 Å is appropriate for NIRSpec PRISM resolution.
mask_lya_emission_width_A (float) – Half-width of the narrow Lyα-emission core mask in rest-frame Å (default 8 Å). Set this small enough not to eat the damping wing. Ignored when fit_lya=True.
fit_range_A (tuple) – Rest-frame fit window (default 1216–3000 Å, matching Pollock+26 — Lyα to just below the Balmer break).
mask_regions_A (list of (float, float), optional) – Extra rest-frame wavelength ranges to ignore.
fit_lya (bool) – Jointly fit a skewed Lyα emission profile.
lya_params (array of length 4, optional) – Pre-determined Lyα profile (A_peak, μ, σ, α) added to the intrinsic spectrum before attenuation.
prior_log_NHI (tuple) – Uniform-prior bounds for the named parameter (Pollock+26 Sect. 3.2 defaults: log_NHI ∈ [18,24], x_HI ∈ [0,1], β_UV ∈ [-4,0]).
prior_x_HI (tuple) – Uniform-prior bounds for the named parameter (Pollock+26 Sect. 3.2 defaults: log_NHI ∈ [18,24], x_HI ∈ [0,1], β_UV ∈ [-4,0]).
prior_beta_UV (tuple) – Uniform-prior bounds for the named parameter (Pollock+26 Sect. 3.2 defaults: log_NHI ∈ [18,24], x_HI ∈ [0,1], β_UV ∈ [-4,0]).
prior_log_F0 (tuple, optional) – Uniform prior bounds on log10(F0) in the same flux units as the input data. Defaults to ±3 dex around the median continuum flux derived from the input spectrum.
prior_log_A_lya (tuple) – Priors used when fit_lya=True.
prior_sigma_lya (tuple) – Priors used when fit_lya=True.
prior_alpha_lya (tuple) – Priors used when fit_lya=True.
n_live (int) – Number of dynesty live points (default 500).
dlogz (float) – dynesty convergence threshold on log Z (default 0.5).
seed (int) – RNG seed.
emission_lines (list[tuple[float, float, float]] | None)
lam_pivot_A (float | None)

Returns:

Posterior samples for every free parameter, median ± 16/84 percentiles, log-evidence, and a 95th-percentile upper limit on log N_HI when the posterior is unconstrained.

Return type:

DLAResult

Notes

Pollock+26 Eq. 4: F = (F0·(λ/λ_pivot)^β + Lyα) · exp(−τ_DLA) · exp(−τ_IGM).
Lyα emission is added before attenuation (Pollock+26 convention).
Reports an upper limit when the posterior touches the prior lower bound (Pollock+26 do this for 21/48 sources).

jwspecfit.dla.tau_DLA(wave_A, log_NHI, z=0.0, b_kms=30.0)[source]

Compute DLA optical depth at each wavelength.

Implements Eq. 1 of Pollock et al. (2026):: tau_DLA = N_HI * sigma_0 * H(a, u)

where sigma_0 = sqrt(pi) * e^2 * f_alpha / (m_e * c * Delta_nu_D), a = Gamma_alpha / (4 pi Delta_nu_D), and u is the dimensionless frequency offset from Lya.

Parameters:

wave_A (np.ndarray) – Observed-frame wavelength in Angstrom.
log_NHI (float) – log10(N_HI / cm^-2).
z (float) – Source redshift (0 for rest-frame spectra).
b_kms (float) – Doppler parameter in km/s (default 30).

Returns:

Optical depth tau(lambda).

Return type:

np.ndarray

jwspecfit.dla.tau_IGM_DW(wave_A, x_HI, z_source, z_min=5.3, cosmo=None)[source]

IGM damping-wing optical depth from a uniformly neutral IGM.

Implements the closed-form solution of Miralda-Escudé (1998), as written in Totani et al. (2006) Eq. 3 and used in Pollock et al. (2026) Eq. 3. Computes the redward Lyα damping-wing absorption by neutral hydrogen of mean fraction x_HI distributed uniformly between redshift z_min and z_source.

Parameters:

wave_A (np.ndarray) – Observed-frame wavelength array in Angstrom.
x_HI (float) – Volume-averaged neutral hydrogen fraction of the IGM.
z_source (float) – Source redshift (upper bound of the IGM integration).
z_min (float) – Lower bound of the IGM integration; default 5.3 (Bosman+22).
cosmo (astropy.cosmology, optional) – Cosmology used to compute the line-centre Gunn-Peterson optical depth. Defaults to Planck18.

Returns:

Optical depth τ_IGM(λ). Zero for pixels redward of the source’s Lyα (the formula returns zero by construction outside the integration range), and a finite damping-wing contribution between λ_α(1+z_source) and longer wavelengths.

Return type:

np.ndarray

Notes

For pixels blueward of λ_α(1+z_source) the photon would have been absorbed inside the GP forest already; the IGM damping-wing formula is not the right description there. Callers should enforce a separate model = 0 cut blueward of source Lyα.

jwspecfit.dla.voigt_H(a, u)[source]

Voigt-Hjerting function H(a, u) via the Faddeeva function.

Computes the exact H(a, u) = Re[w(u + i*a)] where w(z) is the Faddeeva function. This is the function that Tepper-Garcia (2006) approximates analytically; here we use the full solution.

Parameters:

a (float) – Voigt damping parameter.
u (np.ndarray) – Dimensionless frequency offset from line centre.

Returns:

H(a, u) at each u.

Return type:

np.ndarray