Configuration Reference

Thor reads a single YAML configuration file. This page documents the available options.

File Structure

# Global settings
## Required
dataset_type: 'shellmodel'
driver_type: 'mcrtsimulation'

## Optional
device: "cpu"          # options: "cpu", "cpu-openmp", "gpu", "multigpu"
log_level: info        # options: trace, debug, info, warn, error
log_dir: "./logs"

# Dataset section (name must match dataset_type)
shellmodel:
  # ...

Global Settings

Device Options

Dataset Types

• shellmodel - Spherical shell
• powerlaw_shellmodel - Power-law shell profiles
• unigrid - Uniform Cartesian grid
• infiniteslab - Infinite slab
• pointcloud_voronoi - Voronoi tessellation (experimental)
• pointcloud_voronoi_tree - Tree-accelerated Voronoi (incomplete)
• pointcloud_sph - SPH interpolation (incomplete)

Driver Types

• mcrtsimulation - Monte Carlo Radiative Transfer
• raytracer - Forward ray tracing (projections, absorption spectra, tracers, etc)

MCRT Simulation Settings

mcrtsimulation:
  outputpath: output
  overwrite: true
  linename: Lya

  # Photon control
  nphotons_max: 100000
  nphotons_step_max: 100000
  nsteps_per_photon_max: 10000000

  # Physics
  xcrit: 3.0
  acc_scheme: Smith15
  interactor: ResonantWithDustInteractor

  # Stepping
  max_step: 1.0
  min_step: 0.0

  # Boundary conditions
  pbc: false
  hubble_flow: 300.0       # km/s/Mpc (for cosmological simulations)

  # Peeling
  use_peeling: true
  peeling:
    dist_max: 0.6
    tau_max: 100.0

  # Emission model
  emissionmodel:
    # ...

  # Output processors
  outputprocessors:
    # ...

MCRT Parameters

Line Names

The linename selects the resonant line to simulate. It determines the scattering cross-section, the ion density field that drives opacity, and which Cloudy emission/absorption fields to load. Use one of the following abbreviations:

Acceleration Schemes

Only use values other than none with Lyman-alpha.

Interactors

Dust

When interactor: ResonantWithDustInteractor, the dust absorption coefficient field DustAbsorptionCoefficient [cm⁻¹] is built from a top-level dust: block:

dust:
  model: constant           # constant | laursen | colt
  variant: LAURSEN_SMC      # SMC | LAURSEN_SMC       (laursen/colt only)
  k_dust_value: 1.0e-21     # [cm^-1]                   (constant only)
  z_ref: 0.005              # reference metallicity     (laursen/colt)
  f_ion: 0.01               # HII-region dust survival  (laursen/colt)
  t_sputter: 1.0e6          # [K], thermal sputtering   (colt only)
  dust_to_gas_field: D      # optional input field name (colt only)

Byrohl & Nelson 2025 validates only the constant model (infinite slab, Fig. 4). The laursen and colt models implement literature formulae (rascas, COLT) but are not yet independently validated in THOR.

Dataset support: laursen and colt require recipe-system field access and are currently PointCloud-only. Shell, uniform-grid, and infinite-slab datasets support only constant. For infinite slab, the legacy infiniteslab.tau_dust key still works and takes precedence over dust.k_dust_value.

Debug Settings

THOR has two stuck-photon detectors, sharing one disposition vocabulary (off / warn / finish / error):

mcrtsimulation:
  stuck_photon_mode: error          # real stall guard; finish for unattended sweeps
  max_steps_without_scatter: 1000000
  debug:
    step_budget_mode: off           # debug only; also fires on thick scatterers

The no-progress counter resets on every scatter, so thick scatterers never trip it — hence its error default. The step-budget signal fires on any budget-exhausting photon, including legitimately thick ones, so keep it off/warn: finish biases results and error aborts valid runs. (No-progress detection compiles out with THOR_NSTEPS_SINCE_SCATTER_ATTR, default ON. The old debug.stuck_photon_print/stuck_photon_finish keys are deprecated, mapped to step_budget_mode.)

Peeling Settings

When use_peeling: true, peeling photons are generated at each scattering event and propagated towards the observer:

mcrtsimulation:
  use_peeling: true
  peeling:
    dist_max: 0.6       # maximum observer distance (box widths)
    tau_max: 100.0       # maximum optical depth before discarding peel

Emission Models

Single Source

emissionmodel:
  mode: "singlesource"
  singlesource:
    nphotons: 1000
    lum_total: 1e42
    position: [0.5, 0.5, 0.5]
    forced_weight: 1.0  # optional

List Source

emissionmodel:
  mode: "listsource"
  listsource:
    loader: gadget
    nphotons_per_source_min: 1
    nphotons_per_source_max: 1000
    spatial_filter:               # optional, also supported by datasetsource
      mode: sphere                # none | cube | sphere
      center: [0.5, 0.5, 0.5]
      radius: 0.25

Multi Source Options

The listsource and datasetsource generators share these options:

Dataset Source

emissionmodel:
  mode: "datasetsource"
  datasetsource:
    linename: HalphaRecomb        # required: which line/recipe drives the Emissivity field
    emission_position: centroid   # optional: 'centroid' (default) or 'uniform_in_cell'
    spatial_filter:               # optional, also supported by listsource
      mode: sphere                # none | cube | sphere
      center: [0.5, 0.5, 0.5]
      radius: 0.25
    lum_threshold: 1.0e-6              # optional absolute luminosity cut
    lum_threshold_fraction: 0.5        # optional brightest-emitter cumulative luminosity fraction
    lum_threshold_policy: "stricter"   # "stricter" = both cuts, "looser" = either cut

Use uniform_in_cell when the cell linear size is comparable to or larger than the radiative transfer length scale you care about resolving — otherwise the per-cell centroid is fine and avoids the per-photon sampling cost.

Output Processors

Completed MCRT output files store line_name and line_rest_wavelength_angstrom as root HDF5 attributes in each enabled photon stream. The wavelength is the reference used for photon frequency offsets. Python postprocessing uses this value, so adding a line does not require a separate Python wavelength entry. Older files use the configured line.name or linename and the legacy wavelength table; an unknown line without metadata raises an error.

Output processors transform or aggregate photon data. The photon processor (raw photon dump) is the most common.

Common parameters:

Stream Filtering

THOR runs up to three photon streams: original (escaped photons), peel (peeling-off photons), and input (initial photon state). By default, each output processor is active for all streams. The streams parameter restricts a processor to specific streams:

outputprocessors:
  - type: "photon"
    active: true

  - type: "sbmap"
    streams: ["peel"]        # only process peeling photons (default for sbmap)
    nbins_x: 64
    nbins_y: 64
    left_x: 0.0
    right_x: 1.0
    left_y: 0.0
    right_y: 1.0
    weight: "weight"

  - type: "histogram1d"
    # no streams → active for all streams
    field: "dlambda"
    weight: "weight"
    left: -1.5
    right: 1.5
    nbins: 30
    assess_convergence: true
    convergence_threshold: 0.1

Some processors have built-in defaults. The frame-based processors — sbmap, apertureSpectrum, and the histogram1d aperture — need a fixed observer frame: without an explicit frame they default to streams: ["peel"] (the line of sight is inferred from the peeling observer directions), but with an explicit line_of_sight or lines_of_sight they default to streams: ["input"], since a fixed frame is meant for the intrinsic source photons. An explicit streams in the config always overrides these defaults.

HDF5 Compression

By default, compressed output datasets use HDF5's built-in shuffle + gzip filters. Both are part of libhdf5 itself — no filter plugin is needed, and h5py, h5ls, and any other stock HDF5 reader decompress transparently.

To use lossless LZ4 compression instead, set:

mcrtsimulation:
  compression:
    codec: lz4
    level: 5
    shuffle: bit

The same block can be set on an individual output processor to override the inherited settings. Available codecs are lz4, zlib, and raw (uncompressed). Shuffle can be none or byte; LZ4 also supports bit. Existing compress_* settings still work; use either the new block or compress / compress_level / compress_shuffle at each configuration level. Chunk sizing and compress_float32 remain separate options.

LZ4 output requires a Blosc-enabled reader. For Python, install hdf5plugin and import hdf5plugin before reading with h5py. If a selected codec is unavailable in your build, THOR warns and writes uncompressed output. Invalid settings are errors. The defaults below are unchanged.

The default depends on the processor class, because the two write patterns have very different costs:

An explicit compress key beats the class default, and a per-processor value beats a driver-level one.

Without compress_float32, compression is a storage change only: values round-trip bit-for-bit. A raw writer left at its default keeps its historical 1024-row chunks, so the streamed columns' layout and dtype are exactly as before. (That is a per-dataset guarantee. Whole HDF5 files are not reproducible byte for byte between runs in any case — MCRT photon ordering varies with thread scheduling.)

mcrtsimulation:
  compress_float32: true       # inherited by every processor below
  outputprocessors:
    - type: "photon"
      compress: true           # opt the streamed columns in, against their class default
      compress_level: 0        # narrow but do not deflate — the efficient point (see below)
    - type: "sbmap"
      compress_float32: false  # override one key; still compressed by class default
      nbins_x: 512
      nbins_y: 512
      left_x: 0.0
      right_x: 1.0
      left_y: 0.0
      right_y: 1.0

Chunk shape. Chunks are sized by target bytes, trimming the slowest-varying axis first, so a [nlos, nx, ny] map chunks into whole images or contiguous row bands of one image. Sizing matters for the ratio, not just for the cache: deflate compresses each chunk independently, so a chunk far below its window wastes compressibility. On a real 512×512 sbmap, one-row 4 KiB chunks reach 27.0× where 1 MiB chunks reach 31.9× — and 48.2× vs 65.0× once compress_float32 halves the element size.

Measured on a 5000-photon shell-model peel run (512² sbmap, ~194k peel photons):

Ratios are strongly data-dependent. Sparse, smooth arrays such as SB maps compress enormously; a float64 column of noisy per-photon physics values (wavelength offsets, optical depths) has effectively random mantissa bits and compresses only ~1.2×. That asymmetry is why compress_float32 buys much more on the photon columns than on the maps.

Opting the streamed columns in

The class default keeps them uncompressed because gzip there runs inside the writer mutex on the transport hot path. Measured on a deliberately write-bound run — optically thin shell (τ ≈ 1), 2M photons, peeling off, 4M rows / 422 MB of photon columns. Medians of 3 runs; the writer column is time spent inside the HDF5 mutex, which is 22% of wall at baseline:

Reading this:

• gzip on photon columns is a poor trade in isolation. It roughly doubles the wall time of a write-bound run for 2.35×. Dropping to compress_level: 1 saves almost nothing and gives up no size — level is not a useful escape hatch here.
• compress_float32 is the efficient lever: 1.71× for an 11% wall cost. The columns that dominate volume are noisy float64 physics values whose mantissa bits are near-random — deflate does full work on them for ~1.2×, whereas narrowing gives a guaranteed 2× for the price of a type conversion.
• Combine them when size matters most: 4.82× for a 1.57× wall cost.

Both runs above had IO effectively free (tmpfs, and a local disk whose writes fit in page cache), so they isolate CPU cost. On a bandwidth-limited filesystem — cluster scratch, NFS, a shared parallel FS — writing 4.8× fewer bytes shifts the balance back toward compression by an amount these numbers do not capture. Measure your own case with log_level: trace and compare the Writing N entries to disk took ... MB/s lines.

Photon Output

Writes escaped photon data to HDF5.

Histogram 1D

Bins photon field values into a histogram (e.g., spectrum).

The optional aperture restricts the spectrum to a circular region of the image plane, using the same centered projection convention as sbmap (box center → image (0.5, 0.5)). With multiple peeling observers it applies per-observer (each photon is masked against its own line-of-sight frame) and the output is stacked [nlos, nbins], matching the full-field spectrum. This produces the aperture-integrated line profile during the run, without dumping raw photons for offline masking.

The single aperture_radius key is unit-aware:

kpc is accepted as an alias for pkpc. An unknown unit, or arcsec on a non-cosmological dataset, is a config error.

Aperture Spectrum (apertureSpectrum)

Produces a spectrum for several concentric circular apertures at once, in a single binning pass. It is the multi-aperture generalization of the histogram1d aperture: instead of one aperture_radius, you give a list of aperture_radii and get one enclosed spectrum per radius (all flux within that radius), written to /aperturespectrum/.

Each photon is projected onto its observer's image plane and assigned to the radial bucket between the sorted radii (a single upper_bound); a photon beyond the largest radius is dropped. At finalize the buckets are summed inner→outer, so aperture k holds the flux enclosed within radius r_k — the same "inside radius" meaning as the single histogram1d aperture, for every radius simultaneously. Doing the bucket sort once avoids re-masking the photon stream per aperture.

Like the histogram1d aperture, it needs a fixed observer frame. Without an explicit frame it defaults to the peel stream (peeling.lines_of_sight); with an explicit line_of_sight or a lines_of_sight list it defaults to the input stream (intrinsic source photons). An explicit streams overrides either default. Output is [n_apertures, nbins] for a single observer, or stacked [nlos, n_apertures, nbins] with N>1 observers or any explicit lines_of_sight broadcast list (even a one-entry list → [1, n_apertures, nbins]); per-observer frames are in /aperturespectrum/line_of_sight, /e_x, /e_y. The sorted radii and projection geometry are stored as dataset attributes (aperture_radii, n_apertures, aperture_center, center, north).

Physical units are available for peel and input. The peel stream describes radiation directed toward an observer after transport. Input projections describe intrinsic emission, assuming isotropic emission without attenuation or scattering. Each input view receives the full projected source luminosity; THOR does not divide it among viewing directions. Physical units are unavailable for original, whose packets sample different escape directions. The radiation_origin attribute records how to interpret the output.

With units: erg/s, each bin contains the isotropic-equivalent luminosity enclosed by the aperture. Flux modes divide this luminosity by \(4\pi D_L^2\). Mean surface-brightness modes divide it by \(4\pi A(1+z_\mathrm{obs})^4(206265)^2\), where \(A\) is the full circular aperture area in physical cm², including any area outside the dataset. Because brightness is averaged over the whole aperture, increasing its radius need not increase the brightness.

Per-Å modes also divide by the observed wavelength-bin width, \(\Delta\lambda_\mathrm{obs}=(1+z_\mathrm{obs})(\mathrm{right}-\mathrm{left})/\mathrm{nbins}\). Use quantity: dlambda, or dlambda_peel_shifted on the peel stream. The left and right attributes describe the transport wavelength offsets in Å. Flux and surface-brightness modes preserve the absolute normalization; neither divides by the spectrum sum.

THOR treats \(|z_\mathrm{obs}|<10^{-10}\) as zero. Luminosity and mean surface brightness remain valid at zero redshift. Flux requires an explicit luminosity_distance_cm. At positive redshift, THOR can calculate the distance from a valid flat matter-plus-Λ cosmology. You can also supply a distance for local or non-cosmological datasets.

Angular apertures in physical modes use the selected observer distance, with \(D_A=D_L/(1+z)^2\). Thus arcsec radii also require an explicit distance at zero redshift. Apertures in pkpc or box fractions need no distance for luminosity or surface brightness. Observer overrides change output calibration and angular aperture sizes. They leave photon transport and the gas state unchanged.

outputprocessors:
  - type: apertureSpectrum
    units: erg/s/cm2/arcsec2/Angstrom
    aperture_radii: [10 pkpc, 30 pkpc]
    left: -5
    right: 5
    nbins: 80

THOR writes one spectrum per configured processor. Its dataset, such as /aperturespectrum/dlambda, contains values in the selected units and records units, observer_redshift, luminosity_unit_erg_s, and weight_to_spectrum_factors. It also records the available physical scale, distance, aperture areas, and observed bin width.

THOR computes physical bins in float64, including with FP32 transport. With compress_float32: true, it stores the spectrum as float32 if every value fits that numeric range. Otherwise, it warns and stores the entire spectrum as float64. Lossless compression remains available. Raw mode uses its existing storage format and records units: weight.

In MPI runs, aperture spectra sum contributions from all ranks during finalization. Only rank 0 writes the complete spectra and observer frames, in data.0.h5. Each spectrum records mpi_aggregation: sum and mpi_rank_count. The Python read_aperture_spectrum() reader reads these complete products directly and still sums older rank-local spectrum files. It rejects mixed global and rank-local products to prevent double-counting.

For an input-stream aperture spectrum, quantity: dlambda is evaluated in the requested observer direction rather than blindly reusing the stored source-photon direction. THOR caches the emitting cell's bulk velocity once per photon batch and converts the frequency separately for every LOS. With write_inputphotons_restframe: true (the default), it evaluates dfreq_los = dfreq_gas + nu0 (v_cell . los) / c. With write_inputphotons_restframe: false, stored dfreq already contains the Doppler shift along the photon's original direction k, so it evaluates dfreq_los = dfreq + nu0 [v_cell . (los - k)] / c; this replaces the old directional shift instead of applying the velocity twice. This local peculiar-velocity transform does not add the peel photon's final-leg Hubble shift represented by dlambda_peel_shifted.

outputprocessors:
  - type: "apertureSpectrum"
    quantity: "dlambda"
    lines_of_sight: [[1, 0, 0], [0, 1, 0], [-1, 0, 0]]
    aperture_radii: ["1 pkpc", "5 pkpc"]
    nbins: 200
    left: -10.0
    right: 10.0

outputprocessors:
  - type: "apertureSpectrum"
    quantity: "dlambda_peel_shifted"
    aperture_radii: ["1 arcsec", "2 arcsec", "5 arcsec"]
    nbins: 200
    left: -10.0
    right: 10.0

Surface Brightness Map (sbmap)

Projects photon positions onto a 2D image plane perpendicular to the line of sight. Produces an HDF5 dataset at /sbmap/image.

The image frame is computed from the line of sight direction and the north hint vector. By default (peel stream), the LOS comes from the peeling observer directions. Providing an explicit line_of_sight (or a lines_of_sight list) images the intrinsic source photons instead: it routes the map to the input stream automatically — no streams override needed — and a lines_of_sight list broadcasts each source photon into every frame, giving a stacked [nlos, nbins_x, nbins_y] output. An explicit frame combined with a forced streams: ["peel"] is rejected, since the peel observer directions would otherwise silently override the configured frame.

Image coordinates use a centered convention: image_x = (pos - center)·e_x + 0.5 (and likewise for y), so the projection center — the box center by default — lands at image coordinate (0.5, 0.5) for every observer direction. For an axis-aligned observer the unit box projects exactly onto [0, 1] in each image axis, so the default [0, 1] bounds capture the full cross-section. For an off-axis observer the projected box footprint is larger (e.g. a [1, 1, 1] line of sight spans about [-0.207, 1.207]), so photons near the projected corners fall outside [0, 1] and are clipped — widen left/right (or accept the central crop) for diagonal lines of sight. Zooming means narrowing left/right symmetrically around 0.5, or moving center onto the target. For the canonical +z observer this is identical to projecting raw box coordinates. If every photon an observer received falls outside the binning range, a warning is emitted at finalization (the image would otherwise be silently empty).

With multiple peeling observer directions (peeling.lines_of_sight, N>1), one image slab per observer is accumulated and written as a stacked [nlos, nbins_x, nbins_y] dataset, with per-observer frames stored as [nlos, 3] datasets (/sbmap/line_of_sight, /sbmap/e_x, /sbmap/e_y).

To write absolute surface brightness, add these settings to the sbmap processor:

units: erg/s/cm2/arcsec2
observer_redshift: 0.0  # optional; defaults to the snapshot redshift

THOR needs the physical domain length to convert pixel areas. Physical units support two streams:

• peel: radiation directed toward an observer after transport.
• input: intrinsic emission, assuming isotropic emission without attenuation or scattering. Each view receives the full projected source luminosity.

The radiation_origin attribute records which interpretation applies. Physical units are unavailable for original, whose packets sample different escape directions.

/sbmap/image contains surface brightness integrated over the line. The selected units apply to this single image dataset. THOR computes physical pixels in float64, including with FP32 transport. With compress_float32: true, it stores the image as float32 if every value fits that numeric range. Otherwise, it warns and stores the entire image as float64. Lossless compression remains available. Standard luminosity-based emitters assign one weight unit to \(10^{42}\) erg/s. If you change that convention or override packet weights, set luminosity_unit_erg_s to the luminosity represented by one accumulated weight unit.

With the default luminosity scale, a pixel with summed weight \(W\) and projected physical area \(A\) in cm² has surface brightness

This preserves the absolute luminosity without normalizing by the map total. The distance cancels from the conversion, so no cosmological distance integration is needed. At zero observing redshift, the brightness has no cosmological dimming. Changing observer_redshift affects only the output conversion; the gas state and photon transport stay the same.

The image records units, length_unit_cm, observer_redshift, pixel_area_cm2, luminosity_unit_erg_s, and weight_to_image_factor so you can reproduce the conversion. Raw output records units: weight and the physical length when available.

los_indices and los_stride are mutually exclusive, per-processor selectors. They reduce only the SB-map slabs retained in memory and written to disk; peeling photons are still generated and transported for every configured observer, so other processors such as apertureSpectrum continue to receive all observers. Explicit indices must be unique, non-negative, in range, and non-empty. Selection requires peeling observer directions and therefore cannot be used with an inferred non-peel line of sight.

In MPI runs, SB maps sum contributions from all ranks during finalization. Only rank 0 writes /sbmap/image and its observer frames, in data.0.h5. Other rank files retain their raw photon outputs but contain no SB map. The image records mpi_aggregation: sum and mpi_rank_count; do not sum this image again across rank files. Empty ranks contribute zero pixels. Other histograms retain their existing per-rank output behavior.

For a subset selected from a multi-observer context, /sbmap/image keeps the shape [K, nbins_x, nbins_y] and the frame datasets keep [K, 3], including when K is one. The image dataset receives a los_indices attribute whose entries give the original observer index for each output slab.

Dataset-Specific Configuration

See Datasets for detailed configuration of each dataset type.

Resonant scattering phase functions

mcrtsimulation.line.phase_function selects the unpolarized angular distribution for resonant scattering and peeling. The default, auto, selects:

mcrtsimulation:
  line:
    name: FeXXV
    phase_function: auto

Explicit choices are isotropic, rayleigh, half-rayleigh, and lya (Lyα only). An omitted setting and auto both use isotropic scattering when no verified physical default is known, with a warning identifying the fallback. This preserves historical behavior; it does not establish isotropic scattering as the physical default. Explicit choices are respected. THOR warns if you override the physical default or if no physical default is known for the line. For the MgII doublet, an explicit override applies to both components.

MgII follows the scalar interference phase of Seon (2023), Eqs. 8–9, with the line-center limits described by Dijkstra & Loeb (2008), Appendix A. The sampled ion velocity sets the atom-frame frequency for scattering and peeling; polarization is not propagated.

Migrating existing configurations

The default changed from isotropic to auto. Set phase_function: isotropic to preserve isotropic scattering.

The deprecated isotropic_phasefunc override still accepts true for isotropic scattering. With false, it selects Lyα only for the literal Lya name and isotropic scattering otherwise. Do not specify both keys.

Survival estimator

mcrtsimulation.survival_cutoff controls packet survival in dust interactions and resonant lines with radiative branching:

• -1 (default): stochastic destruction. The packet scatters at unchanged weight with the survival probability; otherwise it is removed.
• 0: survival weighting. Multiply the packet weight by the survival probability and scatter.
• A positive finite value: survival weighting with Russian roulette below the cutoff. A packet of weight w survives with probability w/cutoff and takes the cutoff weight; otherwise it is removed. This preserves the expected weight while reducing work on low-weight packets.

For dust, the survival probability is the albedo. For a resonant line, it is the probability of radiative decay back into that transition.

The cutoff is an absolute packet weight. For example, survival_cutoff: 1.0e-5 enables roulette below that weight. Other negative values, infinity, and NaN are rejected.

Alternatively, mcrtsimulation.survival_mode accepts StochasticDestruction (equivalent to -1) or SurvivalWeighting (equivalent to 0). Specify only one of survival_mode and survival_cutoff.