Add dust to galaxy Lyα transfer

Dust can absorb Lyα photons or scatter them into new directions. Resonant scattering changes the paths photons take through dusty gas, so attenuation can depend on source location, wavelength, and aperture. Here we add dust to the three-component emission example and compare the escaping light at the same absolute luminosity scale.

Dust model

We estimate dust opacity from the cutout's gas metallicity and ionization fractions using THOR's laursen prescription, based on Laursen, Sommer-Larsen & Andersen (2009), with local interaction coefficient:

THOR uses \(\sigma_{d,\rm SMC}=3.95\times10^{-22}\,\mathrm{cm^2}\) at Lyα. We choose z_ref: 0.005 and f_ion: 0.01, the fraction of dust surviving in ionized gas. Both \(Z\), supplied by GFM_Metallicity, and \(Z_{\rm ref}\) are metallicity mass fractions.

dust:
  model: laursen
  variant: LAURSEN_SMC
  z_ref: 0.005
  f_ion: 0.01

The dust-enabled configurations also load DustAbsorptionCoefficient and set mcrtsimulation.interactor: ResonantWithDustInteractor. The coefficient includes both absorption and scattering. THOR uses albedo 0.32 and a Henyey–Greenstein scattering phase function with asymmetry parameter 0.73.

The laursen model has no explicit temperature cutoff for sputtering. See the dust reference for the colt alternative, which includes a cutoff.

Run the comparison

First run the dust-free emission example. Use the same demo cutout and download these three complete configurations:

OMP_NUM_THREADS=8 thor lya-recombination-dust.yaml
OMP_NUM_THREADS=8 thor lya-collisional-dust.yaml
OMP_NUM_THREADS=8 thor lya-stellar-dust.yaml

For GPU execution, use the same device and backend settings as in the emission example.

The emissivities, stellar luminosities, packet allocation, and viewing direction match the preceding example. Star-forming gas remains excluded from diffuse emission but contributes dust opacity and resonant scattering. The stellar factor stays at lya_sfr_frac: 1.0: applying an additional empirical Lyα escape factor could count attenuation twice. These runs do not model absorption of ionizing photons before stellar Lyα is produced or an unresolved birth-cloud dust screen.

With plot_lya_components.py and thorplot.py in the working directory, run:

uv run --python 3.11 --with h5py --with numpy --with matplotlib plot_lya_components.py --dust

Inspect the result

The upper row shows the dust-free escaping light; the lower row adds dust. A shared color scale preserves attenuation and the relative component strengths. The stellar component is more strongly attenuated: its all-direction photon escape fraction is 38%, compared with 80–84% for diffuse gas emission.

The spectra sum the three components within 10 and 30 kpc. Their normalization is unchanged between runs, so a lower curve indicates less escaping light in that direction and spectral interval. Scattering can redirect light into a particular pixel, aperture, or wavelength bin; attenuation need not be uniform there.

The script calculates the all-direction photon escape fraction as \(f_{\rm esc}=\sum w_{\rm escaped}/\sum w_{\rm emitted}\), using photon-rate weights. The sums include all packets in the selected component, or all three components for the combined fraction.

These runs give approximately:

Increase packet counts before interpreting small spectral differences. Vary f_ion to test dust survival in ionized gas, or z_ref to change opacity at fixed snapshot metallicity; increasing z_ref decreases the opacity. For quantitative escape fractions, also check convergence against reduced or disabled core skipping (acc_scheme: none, use_atau0nl: false), since changes to photon paths can matter for absorption.