nebular_em

A Python library for computing nebular emission line luminosities from simulation grid data and deriving physical conditions (T_e, n_e) from line ratios using PyNeb.

Why this library?

• Compute luminosities for lines your radiative transfer code does not produce (e.g. UV OIII], CIII], NIV]).
• Cross-check code-output luminosities with an independent pyneb-based calculation.
• Perform observer-style diagnostics (T_e, n_e) from any set of integrated line luminosities — whether computed here or taken directly from your RT code.

• Python ≥ 3.10
• pyneb ≥ 1.1
• numpy, scipy, h5py

# from the repo root
pip install -e nebular_em/
# or simply add the repo root to PYTHONPATH:
export PYTHONPATH=/path/to/nebular-emission:$PYTHONPATH

COLTSnapshot handles all the file I/O and ion-density conversion, then hands off to the core classes:

from nebular_em import COLTSnapshot, Diagnostics

snap_dir = '/path/to/colt_dir'   # directory with states-teq_NNNN.hdf5 + colt_NNNN.hdf5

with COLTSnapshot(snap_dir, snapnum=150, Zgas_Zsun=0.1) as snap:
    # Build a GridEmission from the snapshot
    grid = snap.to_grid_emission(['OIII', 'OII', 'HII', 'SII'])

    # Compute integrated luminosities over the ionized region
    lines = ['OIII_5007', 'OIII_4363', 'OII_3726', 'OII_3729', 'Ha', 'Hb']
    L = grid.integrated(lines, mask=snap.ionized_mask())

    # Observer-style diagnostics
    Te, ne = Diagnostics(L).iterate_Te_ne(
        Te_lines=('OIII_4363', 'OIII_5007'),
        ne_lines=('OII_3726',  'OII_3729'),
    )
    print(f"T_e = {Te:.0f} K,  n_e = {ne:.1f} cm⁻³")

from nebular_em import GridEmission, Diagnostics, LINE_LIBRARY

# Te, n_e, n_OIII, ... are (N,) arrays; vol_cell is a float or (N,) array
grid = GridEmission(
    Te=Te,
    ne=n_e,
    ion_densities={'OIII': n_OIII, 'OII': n_OII, 'SII': n_SII, 'HII': n_HII},
    vol_cell=vol_cell,
)

lines    = ['OIII_5007', 'OIII_4363', 'OII_3726', 'OII_3729', 'SII_6716', 'SII_6731', 'Ha', 'Hb']
cell_L   = grid.cell_luminosities(lines, mask=ionized)   # dict of (N,) arrays [erg/s]
total_L  = grid.integrated(lines, mask=ionized)          # dict of floats [erg/s]

diag = Diagnostics(total_L)    # also accepts COLT-output luminosity dicts
Te, ne = diag.iterate_Te_ne(('OIII_4363', 'OIII_5007'), ('OII_3726', 'OII_3729'))

COLTSnapshot(snap_dir, snapnum, Zgas_Zsun,
             N=512, X_H=0.76, thermal_eq=True)

Reads states-teq_NNNN.hdf5 and colt_NNNN.hdf5 from snap_dir. All arrays are loaded lazily and cached on first access. Use as a context manager (with COLTSnapshot(...) as snap) to ensure HDF5 file handles are closed.

ion_densities(species) → dict[str, ndarray]

Compute number densities [cm⁻³] for a list of ion species using solar abundances (Grevesse et al. 2010) scaled by Zgas_Zsun. Results are cached. See ION_MAP for all supported keys.

ions = snap.ion_densities(['OIII', 'OII', 'HII', 'SII', 'CIII', 'NIV', 'NeIII'])

ionized_mask(threshold=0.2) → ndarray

Boolean mask of cells with x_HII > threshold.

to_grid_emission(species, **kwargs) → GridEmission

Convenience wrapper: calls ion_densities(species) and constructs a GridEmission. **kwargs are forwarded to GridEmission (e.g. Te_range, ne_range).

GridEmission(Te, ne, ion_densities, vol_cell,
             Te_range=(3e3, 1e5), ne_range=(1e-2, 1e7),
             n_Te=80, n_ne=80)

Emissivity tables are built lazily on first request per line and cached for the lifetime of the object.

cell_luminosities(lines, mask=None) → dict[str, ndarray]

Per-cell luminosity [erg/s]. lines is a list of LINE_LIBRARY keys or a custom {name: LineSpec} dict. Cells outside mask are zero.

integrated(lines, mask=None) → dict[str, float]

Sum of cell_luminosities over all (masked) cells.

Diagnostics(luminosities: dict[str, float])

Accepts any {line_name: L_total [erg/s]} dict — from GridEmission.integrated() or directly from an external code.

electron_density(line1, line2, Te=1e4) → float

Derive n_e [cm⁻³] from L[line1]/L[line2] at the assumed T_e. Both lines must be from the same ion (same pyneb atom).

electron_temperature(line1, line2, ne=100.) → float

Derive T_e [K] from L[line1]/L[line2] at the assumed n_e.

iterate_Te_ne(Te_lines, ne_lines, Te0=1e4, ne0=100., n_iter=5) → (float, float)

Iterates the T_e and n_e diagnostics to self-consistency. Returns (Te, ne).

LineSpec(atom, lev_i, lev_j, ion_key, wav, rec=False)

Frozen dataclass describing one emission line. Used to extend LINE_LIBRARY with custom lines (see Adding a custom line).

The UV OIII lines (1661, 1666 Å) use Aggarwal & Keenan (1999) collision strengths via a separate pyneb atom instance.

All ionization stages tracked in the COLT states file are available:

Number densities are computed as:

• H: nH × x
• He: nH × x × X_He / (4 X_H) where X_He = 1 − X_H
• Metals: nH × x × (n_X/n_H)_solar × Zgas_Zsun

with solar abundances from Grevesse et al. (2010).

Use LineSpec directly and pass a {name: LineSpec} dict to cell_luminosities / integrated:

import pyneb as pn
from nebular_em import GridEmission, LineSpec

# Example: [Ne III] 3869 Å (lev_i=3, lev_j=1)
ne3 = pn.Atom('Ne', 3)
custom_lines = {
    'NeIII_3869': LineSpec(atom=ne3, lev_i=3, lev_j=1, ion_key='NeIII', wav=3869.),
}
# ion_densities must include 'NeIII'
grid = GridEmission(Te, ne, ion_densities={'NeIII': n_NeIII}, vol_cell=vol_cell)
cell_L = grid.cell_luminosities(custom_lines)

With COLT data, ion_densities is already provided by COLTSnapshot:

with COLTSnapshot(snap_dir, 150, Zgas_Zsun=0.1) as snap:
    grid = snap.to_grid_emission(['NeIII'])
    cell_L = grid.cell_luminosities(custom_lines, mask=snap.ionized_mask())

To use a custom line in Diagnostics, call pyneb directly:

ratio = total_L['NeIII_3869'] / total_L['NeIII_3967']
ne3.getTemDen(ratio, tem=1e4, to_eval="I(3, 1) / I(2, 1)")

• Tables are log-spaced in both T_e and n_e and interpolated bilinearly in log–log space.
• fill_value=None means the interpolator extrapolates outside the grid bounds; set Te_range/ne_range to cover your data.
• Some pyneb recombination tables (e.g. Storey & Hummer 1995 for H I) only cover T_e ≲ 30 000 K. Cells outside this range are filled by nearest-neighbour extrapolation from the table boundary; errors are small for lines whose emissivity is a smooth power law in T_e.
• Processing is chunked at 5 × 10⁶ cells to limit peak memory.
• Emissivity tables are cached per LineSpec instance — instantiate GridEmission once and call cell_luminosities / integrated multiple times at no extra cost.

If you are not using COLT, convert simulation ion fractions manually:

PROTONMASS = 1.6726e-24   # g
X_H = 0.76

nH     = rho * X_H / PROTONMASS             # hydrogen nucleon density [cm⁻³]
O_solar = 0.00573805 * 1.008 / (X_H * 15.999)   # n_O/n_H at solar Z (Grevesse+2010)
n_OIII  = x_OIII * nH * O_solar * (Z_gas / Z_sun)

The solar_number_abundance(element, X_H) helper from nebular_em performs this calculation for any element in SOLAR_MASS_FRAC.