MCRT simulation and sources

Select driver_type: mcrtsimulation for Monte Carlo radiative transfer and place its settings under mcrtsimulation:. See General configuration for global and shared settings, and Output processors for photon streams, spectra, images, and compression.

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
  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

hubble_flow_at_z, omega_matter, omega_lambda, and hubble_mode are accepted but do not affect MCRT calculations. MCRT transport uses hubble_flow directly; output distance conversions use the dataset cosmology.

Output policy is validated before runtime setup. At output initialization, overwrite: false rejects a nonempty destination before removing files. Set overwrite: true deliberately when reusing an output directory. Compression, processor variants and convergence options retain their separate settings.

Simulation restart (mcrtsimulation.restart: true) is currently unsupported and is rejected during preflight. Existing output files remain intact, regardless of their redshift law or whether they contain law metadata.

Source seeds follow mcrtsimulation.seed → emissionmodel.seed → the selected generator’s seed, with each explicit override taking priority. Source seeds must be integers from 0 through 4294967295. Overridden emission seeds and seeds for inactive generators are ignored. The driver seed remains validated for transport.

emissionmodel.mode is required and must be singlesource, listsource, datasetsource, or planesource. Source selection and the effective seed are validated before dataset initialization; YAML key names and MPI source streams are unchanged.

Peel merging requires a single observer. Either nonzero merge threshold enables it, including a negative peel_merge_dist whose square gives the distance threshold. Despite its legacy name, peel_merge_dfreq is in Å; its signed wavelength offset is converted to an absolute frequency difference. Both flat merge keys are validated even when peeling is disabled. Non-finite or unrepresentable thresholds and incompatible observer counts fail preflight.

use_atau0nl: false initializes the nonlocal opacity cache to zero unless an explicit populate mapping supplies AtimesTau0_NL. Such a mapping also takes precedence when computation is enabled.

Transport

Stepping

Spatial propagation keys remain flat under mcrtsimulation. Their values are validated before initialization. After conversion to simulation precision, dist_max must be finite and a positive value must retain a finite, nonzero square. Zero remains an explicit zero cap. Dataset support for automatic limits and periodic-ray stopping conditions are checked during runtime setup.

Hubble flow

The Hubble settings also remain flat under mcrtsimulation. Their values must be finite and representable in simulation precision; nonzero values must remain nonzero after conversion. The rate is converted from km/s/Mpc to s⁻¹ in double precision before narrowing. hubble_flow_discover: true conflicts with an explicit hubble_flow, including zero. Discovery currently has no dataset resolver and fails during initialization; supply the rate explicitly.

With hubble_flow enabled, MCRT defaults to the linear law ν(s) = ν_emit − ν₀ H s/c, where s is accumulated physical path length. To opt into fractional redshifting, ν(s) = ν_emit × exp(−H s/c), build with:

./build.sh --no-clean --cmake-arg -DTHOR_MCRT_FRACTIONAL_HUBBLE=ON

Use OFF to return an existing build to the linear law. This is a CMake option, not a YAML setting. The same selection applies to opacity, scattering, final-flight output, continuum targeting and the max_dlambda distance limit. Output files record mcrt_hubble_redshift_law as linear or fractional.

Both laws keep H, the physical box scale, gas properties and stellar emission fixed. Fractional redshifting defines a frozen-background approximation; it does not evolve the cosmological scale factor or the snapshot.

The stored photon frequency offset is referenced to the local cosmological frame at its most recent emission or scattering. For doublets, the primary line defines this stored offset; both components see the same absolute photon frequency and its selected Hubble shift in opacity and scattering. Peculiar and thermal Doppler shifts still use the interacting component's wavelength. Each scattering incorporates the redshift over the entire incoming flight, including travel across periodic boundaries, before starting a new flight. Peeling photons start in the same frame; use dlambda_peel_shifted to include the redshift over their final flight when analysing their output.

Termination

Distance limits

max_dlambda derives a flight-distance cutoff from the redshift of a line-centre photon; it is not a per-photon wavelength termination condition. Zero disables this proxy. Signed offsets retain the legacy absolute-distance interpretation and must resolve a positive wavelength and absolute frequency for the selected line. Wavelength inputs are converted from Å to cm in double precision before narrowing to the simulation precision.

When this proxy is active, original photons use the smaller of its distance limit and an enabled dist_max. Peeling photons apply peeling.dist_max independently in the same way. Negative distance caps are disabled; zero remains an explicit zero cap.

Wavelength cutoff

For a cutoff at the packet's actual wavelength, use:

mcrtsimulation:
  pbc: true
  hubble_flow: 70  # km/s/Mpc
  termination:
    redward_dlambda_angstrom: 10  # illustrative margin; check convergence

Omitting termination.redward_dlambda_angstrom disables the cutoff. The offset and Hubble flow must both be positive and finite. Source packets finish at λ₀ + the configured offset in the local cosmological frame, including the current flight's Hubble shift. For doublets, λ₀ is the primary line's wavelength. Scattering before the endpoint starts a new flight with the outgoing frequency. Packets already beyond the cutoff finish without further transport, but still contribute emission and direct peels. The cutoff is independent of the injected wavelength range and spatial termination.mode. Choose the margin by checking that further redward transport does not appreciably change the result.

THOR clips the last step before evaluating opacity and marks completed packets with FREQUENCY_LIMIT. Their stored dfreq references the last emission or scattering; use shifted output fields to include the final flight. The cutoff uses the selected linear or fractional Hubble law. A non-finite cutoff flight distance aborts the run with INVALID_FREQUENCY_LIMIT. Other enabled stopping conditions also apply, including dist_max and the legacy max_dlambda distance proxy. Omit these distance caps to stop on frequency alone.

Set peeling.redward_dlambda_angstrom independently to integrate emission and scattering peels into the red wing, using the source cutoff's frame and Hubble law:

mcrtsimulation:
  pbc: true
  hubble_flow: 70  # km/s/Mpc
  termination:
    redward_dlambda_angstrom: 60
  use_peeling: true
  peeling:
    redward_dlambda_angstrom: 60  # illustrative margin; check convergence
    dist_max: -1                # disable the peel distance cap

The offset and Hubble flow must be positive and finite. Omit max_dlambda and disable peeling.dist_max to avoid stopping before resonance; peeling.tau_max still allows early removal. Increase the redward margin until the result converges.

For periodic IFUs, use spectral.frame: box_center to bin the launch event. Binning the endpoint wavelength would pile up flux at the cutoff.

Spatial termination

mcrtsimulation.termination limits transport to a cube or sphere, in box coordinates.

Geometry must be finite and remain nondegenerate in the simulation's floating-point precision. Fields for an unselected mode are ignored with a warning.

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 peel-ray travel distance (box widths)
    tau_max: 100.0       # maximum optical depth before discarding peel

Omitting peeling, using an empty mapping, or using a null section uses the defaults. A negative dist_max disables the distance limit. An empty observer list is an error. When use_peeling: false, a supplied peeling section is ignored with a warning.

max_peel_in_flight and drain_transport_skip apply only with fanout_at_scatter: false. Drain thresholds additionally require drain_transport_skip: true. Supplying these settings outside their active mode produces a warning and their values are ignored.

Stop conditions and watchdogs

MCRT uses separate stop reasons for ordinary work limits, navigation recovery, and watchdog events.

The two configurable policies accept off, warn, finish, and error. warn logs and permits continuation; finish terminates affected photons and continues the run; error aborts the run. Fatal transport failures always abort, even when a coincident boundary marked the packet finished. These checks run before configurable dispositions or host boundary validation.

mcrtsimulation:
  stuck_photon_mode: error
  max_steps_without_scatter: 1000000
  debug:
    step_budget_mode: off

The no-scatter counter resets on scattering and at each kernel launch. It can count positive-distance steps, so reaching the cap is a watchdog event rather than proof of zero spatial progress. Peel photons never scatter, making this a consecutive-step cap for them. Budget exhaustion also occurs in legitimate thick scatterers: finish can bias results, and error can abort valid runs.

Set max_steps_without_scatter: 0 or stuck_photon_mode: off to disable the no-scatter detector. A cap at or above nsteps_per_photon_max clamps to one less than that budget; if the budget is at most one, detection is disabled. These adjustments and inactive settings produce warnings. This does not disable explicit navigation failures or the persistent MCRT zero-length detector. The latter reports GEOMETRY_STALLED after 64 consecutive zero-length transport steps across kernel relaunches and MPI transfers. Positive travel or an actual interaction resets its per-photon counter; peel-buffer backpressure does not count as a transport step. This guard always aborts and is independent of the two configurable policies.

The remaining mcrtsimulation.debug controls are:

An explicit step_budget_mode overrides the deprecated stuck_photon_print and stuck_photon_finish keys. Otherwise, stuck_photon_finish: true selects finish; if it is false or absent, stuck_photon_print: true selects warn. Photon tracing remains configured separately under mcrtsimulation.trace.

Photon tracing

Tracing controls live under mcrtsimulation.trace and require a build with THOR_DEBUG_PHOTONTRACE=ON.

Active block sizes must fit the device index type. The overhead factor must be finite and nonnegative. It is ignored with a warning when blocksize: 0; both values are ignored with warnings when tracing was not compiled in. Allocation limits that depend on the launch size are checked at runtime.

Timing

These options belong under mcrtsimulation. Durations are in seconds.

The kernel target must be finite and nonnegative. Zero retains the minimum adaptive step budget; it does not disable transport. Probe counts are ignored with a warning when probing is disabled; the scheduler margin is likewise ignored when use_slurm_timelimit: false.

When enabled and present, TIMELIMIT overrides runtime_max. It must contain whole seconds and exceed the safety margin. Runtime setup computes the effective limit as TIMELIMIT - runtime_slurm_safetymargin; --validate-config checks YAML inputs only. Run limits are checked between steps, so they are not hard deadlines.

Memory budgets

These keys belong under mcrtsimulation:

Fractions must be finite and within [0, 1]; allocation counts are nonnegative integers. The peel factor must fit in uint32_t. With fanout_at_scatter: true, it is raised during setup to fit all observers, including when configured as zero. With fanout_at_scatter: false, it must be positive. Zero remains valid for the dynamic factor and minimum allocation. Runtime checks enforce the allocation's index limits once the launch size is known.

The four peel-only options are ignored with a warning when use_peeling: false.

Resonant physics

Line names

Use line.name or the compact linename spelling. Both may be present with the same value; conflicting names are rejected. A partial line: section can accompany linename to set phase or deuterium options.

line.phase_function: auto uses the existing atomic default. An explicit lya phase requires a Lyα transition. The deprecated line.isotropic_phasefunc cannot accompany phase_function; its historical false behavior selects the Lyα phase only for the exact name Lya and isotropic scattering for other names.

A positive line.deuterium_fraction enables coupled H/D transfer and disables core skipping. It requires hydrogen Lyα and a single-line or doublet resonant interactor. The ratio must be finite, nonnegative and representable in the simulation precision. With deuterium disabled, ResonantDoubletInteractor requires the exact primary name MgII. These checks run during preflight; atomic constants and scattering algorithms are unchanged.

The line name selects the atomic properties used for transport. Density and emissivity recipes are selected separately. Supported short names include:

Higher Lyman lines use their table labels: Lyb, Lyc, Lyd, Lye, then Ly6 to Ly39 for levels 7p to 40p. Their upper level np also decays to lower s and d levels, so the Voigt damping uses its total radiative width Γ_np. A scattered packet returns in the same line with probability A(np → 1s)/Γ_np, for example 0.882 for Lyβ; otherwise it leaves the transported line, and the cascade products are not followed. Widths and per-absorption cascade outcomes come from src/LymanCascadeTableGenerated.h, generated by scripts/gen_lyman_cascade_header.py. Sightline absorption from tau: blocks and from thor.absorption uses the same widths.

Ion density

By default, Lyman lines use HIdensity; the listed metal ions use their existing cloudydensity_<ion> fields. Existing field aliases still apply. Loading another provider does not change these defaults.

Select a different recipe explicitly, for example with the AtomDB tables:

mcrtsimulation:
  line:
    name: FeXXV_w
    density_field: atomdbdensity_FeXXV

The dataset must configure the selected recipe. Gadget, Tipsy and octree loaders request it during setup. Explicit selections require matching species metadata and physical number-density units, \(\mathrm{cm}^{-3}\). Cloudy, AtomDB and derived HI recipes provide this metadata; raw snapshot fields without it cannot be selected explicitly. Old field caches may need regeneration. Emission fields remain independently selected by the source configuration.

Full line aliases such as H 1 1215.67A and registered NIST labels such as CIV 1548 select the same density defaults as their species. Each NIST entry retains its own atomic constants. Species without a default require an explicit field. Arbitrary wavelength spellings and Cloudy emissivity blends such as Blnd 2798.00A are not transport-line selectors.

A missing required density stops initialization. Synthetic setups can retain density_check_lenient: true when Density already contains the absorbing population. This flag never replaces a missing or incompatible explicitly selected field. Do not combine line.density_field with populate.Density.

Outputs record the resolved line, atomic-data source, population approximation, and selected density field. Species and recipe provenance are also recorded when available. Legacy fields without species metadata retain their existing contract with a warning. Transport keeps its existing approximation of using the ion density as the absorbing population.

Acceleration schemes

Only use values other than none with Lyman-alpha.

xcrit sets a minimum only for Laursen09; supplying it with none or Smith15 produces an inactive-setting warning. It must be finite. Negative values clamp to zero with a warning. The deprecated xcrit: -1 also selects Smith15 when acc_scheme is absent or empty; use acc_scheme: Smith15 and omit xcrit instead.

Core skipping is disabled for destructive lines and the doublet interactor; explicit acceleration settings then produce a warning during preflight.

Interactors

Lyman-series transport

LymanSeriesInteractor transports Lyα together with the higher H I Lyman lines. It requires linename: Lya: photon frequencies keep the Lyα reference, and each higher line sits at its fixed offset. A packet blueward of Lyβ can pass through the higher resonances, scatter in them, or be converted to a lower Lyman line.

mcrtsimulation:
  linename: Lya
  interactor: LymanSeriesInteractor
  lyman_series:
    n_max: 20       # highest upper level np; 2 transports Lyα only
    tau_tol: 1.0e-8 # skip a line whose per-step optical-depth bound is lower
    eps_rel: 1.0e-4 # skip a non-nearest line below this fraction of the nearest line

Each step evaluates Lyα and the lines that can matter: the line nearest in velocity, if the photon is in its core or its Lorentzian-wing bound exceeds tau_tol, and other lines that also exceed eps_rel times the nearest line's optical depth. Paths that only meet Lyα therefore reproduce SimpleInteractor and its random-number sequence.

An interaction selects the absorbing line by its share of the opacity. Lyα events behave as with SimpleInteractor, including core skipping. An absorption into np with n ≥ 3 draws one outcome of the radiative cascade from src/LymanCascadeTableGenerated.h:

• resonant return in the same line, coherent in the atom frame and isotropic;
• conversion to a lower Lyman line m, emitted isotropically at line centre in the frame of the absorbing atom;
• decay through 2s by two-photon emission, which absorbs the packet.

Balmer and other non-Lyman cascade photons are not transported. Peels of higher-line events follow the sampled outcome, as isotropic emission from the absorbing atom.

Build with -DTHOR_MCRT_FRACTIONAL_HUBBLE=ON for Hubble flow. The linear law redshifts every photon at the Lyα rate, which biases higher-line optical depths high by ν_n/ν_Lyα; THOR warns at startup in that case. Resolve each resonance crossing within a cell by limiting max_step so that the Hubble shift per step is a fraction of the Doppler parameter. continuum_advance is not supported.

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.

Deuterium scattering

Set mcrtsimulation.line.deuterium_fraction to include neutral deuterium in Lyα transfer. The value is the number ratio nDI / nHI; it defaults to 0 (disabled). For example:

mcrtsimulation:
  interactor: SimpleInteractor
  line:
    name: Lya
    deuterium_fraction: 3e-5

THOR multiplies the neutral hydrogen density by this ratio. Deuterium shares hydrogen's temperature, bulk velocity, and turbulent broadening; its larger mass gives a narrower thermal profile at 1215.34 Å (about 82 km/s blueward of H Lyα). Scattering and peeling use the selected isotope. Phase overrides apply to both.

Supported interactors are SimpleInteractor, ResonantDoubletInteractor, and ResonantWithDustInteractor (to include dust). The fraction must be finite, nonnegative, and representable in the build's precision. Enabling deuterium disables core skipping, increasing runtime at high optical depth. Output records deuterium_fraction; frequency offsets remain relative to H Lyα.

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], finite and > 0       (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: Gadget and Tipsy point clouds and RAMSES/ENZO/ART octrees use the same top-level model selection and parameters. All three recipes honor field_filters.DustAbsorptionCoefficient on the final opacity. Explicit model requests cannot fall back to the legacy nH × Z recipe; legacy sigma_per_H_solar/z_solar keys remain supported but cannot be mixed with modern model parameters in one dust block. Ordinary shellmodel accepts constant, with opacity confined to the shell; laursen and colt are rejected. infiniteslab also accepts constant; its legacy infiniteslab.tau_dust key still works and takes precedence over dust.k_dust_value.

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.

For conservative lines without a dust interactor, these settings have no effect; explicitly supplying either produces a warning during preflight.

Emission models

Single source

mcrtsimulation:
  emissionmodel:
    mode: "singlesource"
    singlesource:
      nphotons: 1000
      lum_total: 1e42
      position: [0.5, 0.5, 0.5]

nphotons, lum_total, and position are required. The photon count is a positive global budget distributed across MPI ranks. Luminosity is finite and nonnegative, in erg/s; positions use dataset user coordinates and may lie outside the unit cube. A positive forced_weight overrides the luminosity-derived packet weight; zero or negative values leave it active.

Omit direction, or set it to [0, 0, 0], for isotropic emission. A nonzero three-vector selects directed emission and is normalized during setup. With that direction, cone_half_angle_deg selects a cone from 0 to 180 degrees; omitting it or using a negative value gives a collimated beam. Cone settings without a direction are ignored with a warning.

Preflight validates these inputs before opening output files or initializing sources. It checks finite values, array sizes, and whether luminosity, packet weights and direction normalization can be represented in the build's floating precision. Dataset geometry and MPI-dependent photon-ID capacity remain runtime checks.

List source

This example uses Lya_StellarParticle_Luminosity for stellar Lyα emission from a cosmological Gadget snapshot with stellar formation times and masses. Replace /data/snapshot with your snapshot path.

mcrtsimulation:
  linename: Lya
  emissionmodel:
    mode: listsource
    listsource:
      loader: gadget
      gadget:
        path: /data/snapshot
        particle_type: PartType4
        luminosity_field: Lya_StellarParticle_Luminosity
      nphotons_per_unit_lum: 100
      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

Plane source

mode: planesource places equal-luminosity emitters at the cell centers of a rectangular grid. Each emitter radiates isotropically. Set these options under emissionmodel.planesource:

Negative extents reverse sampling order; zero extents collapse that grid axis. The generator clips emitters to its ownership box and applies the shared source filters below. MCRT checks geometry and luminosity conversion before allocation; coordinates must remain finite in the build's precision, grid counts must be positive without overflowing their product, and positive luminosity must remain nonzero per emitter.

Multi source options

For listsource, datasetsource, and planesource, nphotons_per_unit_lum sets packets per luminosity unit (normally \(10^{42}\) erg/s). Each positive source receives its luminosity times this rate, rounded up and bounded by nphotons_per_source_min and nphotons_per_source_max. The rate and minimum default to 1. target_nphotons is no longer accepted; ionizing RT uses ionizing_rt.nphotons.

MCRT validates budgets and luminosity cuts before loading the catalog. Ratios must be finite, positive, and within the build's normal floating-point range. Cuts must be finite; a fractional cut must be at most 1, with negative values disabling it. Only an explicitly supplied absolute cut participates in the looser policy.

Spatial filters are also validated before loading the catalog. Cube bounds must remain finite and satisfy min < max in simulation precision; the minimum is included and the maximum excluded. A sphere includes its boundary and defaults to center [0.5, 0.5, 0.5]. Its radius must be positive with a finite, nonzero square in simulation precision. Geometry fields for other modes are ignored with a warning; unknown keys inside spatial_filter are errors.

Restart paths are checked as strings during preflight; checkpoint files are opened at runtime. Each {rank} in either path is replaced by the MPI rank. With multiple ranks and no placeholder, .rank<N> is appended. Empty paths disable the corresponding operation. These checkpoints resume source emission; they do not enable full simulation restart.

Dataset source

datasetsource creates emitters from the dataset's Emissivity field and cell volumes. Load this field with the dataset or select it through mcrtsimulation.populate.Emissivity.

mcrtsimulation:
  emissionmodel:
    mode: "datasetsource"
    datasetsource:
      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.

MCRT validates the emission_position choice before dataset loading. Inactive dataset-source positioning is ignored with a warning. The concrete dataset's support for uniform_in_cell is checked during source preparation, before allocating emitters.

Analytic emission spectra

emissionmodel.spectrum accepts the shapes delta, gaussian, doublegaussian, uniform, and gaussian_with_continuum. Omission selects a zero-shift delta. The default constant recipe reads parameters from spectrum.constant; thermal requires a Gaussian and reads temperature in Kelvin.

The default global scope requires the selected recipe subsection. For a Gaussian mixture, mus and sigmas must each contain two entries. The optional weights defaults to [0.5, 0.5]; a global empty weight list also means equal weights. local scope allows an omitted recipe subsection, uses zero-valued fallback parameters and equal mixture weights, and applies the recipe per emitter. Local thermal recipes read emitter temperatures at runtime.

Velocities are in km/s, and equivalent width is in Angstrom. Preflight checks selectors, array sizes, finite values, nonnegative widths and temperatures, ordered uniform edges, and nonnegative mixture weights with a finite positive sum in the sampler's precision. The existing velocity limit is c/2, counting Gaussian tails to ten widths. Global continuum spectra require positive equivalent width; nonpositive local equivalent width disables the continuum.

The legacy scalar syntax, such as spectrum: gaussian with a sibling gaussian_spectrum: {mu: 0, sigma: 100}, remains supported. Unknown active parameter keys are errors; supplied inactive settings produce warnings.

Packet weights

An MCRT packet weight is a photon rate in line units. One unit is the photon rate of \(10^{42}\) erg/s at the line rest frequency ν0, i.e. \(10^{42}\) erg/s / (hν0). Frequency changes from scattering and Hubble flow preserve this photon rate. Absorption and survival weighting change it; peel contributions also include angular weighting and attenuation toward the observer.

Emitter luminosities and spectral shapes describe energy: a spectrum is L_ν in \(10^{42}\) erg/s, for lines and continua alike, as for BPASS L_λ. At launch, a packet emitted at ν_rest receives its luminosity share times ν0/ν_rest, so its rest-frame energy is that share. A line at ν0 keeps its share exactly; a Doppler-broadened line changes by about σ/c (1e-4 at 30 km/s). Spectral velocities are in km/s and may reach at most c/2 (150,000 km/s, counting Gaussians to 10σ); setup rejects larger offsets, which usually mean a unit slip.

Energy outputs multiply each contribution by ν_recorded/ν0. Binned outputs already include this conversion, even with units: weight; that setting omits physical calibration, not frequency conversion. The recorded frequency depends on the output:

• histogram1d, histogram2d and apertureSpectrum use the endpoint frequency for dlambda_peel_shifted, the box-centre frequency for dlambda_boxcenter, and the stored frequency for any other quantity;
• sbmap uses its cube's frame. Without a cube it uses the box-centre frequency, which does not depend on how far peels fly; on the original stream it uses the stored frequency. The image attribute energy_frequency names the choice;
• intrinsic apertureSpectrum and sbmap outputs on the input stream use each observer's frequency. Intrinsic sbmap outputs need a configured unit line_of_sight or peeling.

A histogram2d whose two axes record different frequencies is rejected, for example dlambda_peel_shifted against dlambda_boxcenter or stored dlambda. The photon attributes energy_factor, energy_factor_peel_shifted and energy_factor_boxcenter give these factors. Raw photon output stores the unscaled weight. Convert it with photons/s = weight × \(10^{42}\)/(hν0) and erg/s = weight × \(10^{42}\) × (ν0 + dfreq)/ν0. Output files record mcrt_packet_weight: photon_rate, mcrt_weight_unit_erg_s and mcrt_weight_reference_frequency_hz. Outputs from earlier THOR versions lack these attributes; their weights kept each packet's launch energy. The Python analysis guide shows how to read both conventions and calculate a photon escape fraction or luminosity.

Native BPASS stellar continuum settings

Sample a stellar continuum directly from the native BPASS spectra:

mcrtsimulation:
  emissionmodel:
    mode: listsource
    bpass_continuum:
      enabled: true                              # default: false
      spectra_path: /data/bpass-spectra.h5        # required when enabled
      wavelength_range_angstrom: [1025.72, 1215.67]  # default: direct Lyα pumping band
    listsource:
      loader: gadget
      nphotons_per_unit_lum: 100
      gadget:
        path: /data/snapshot
        particle_type: PartType4
        mass_has_h_factor: true                  # TNG, including z = 0

This opt-in mode requires listsource. Omit emissionmodel.spectrum and luminosity_field: the native bpass-spectra-v1 file supplies both. It is separate from the derived-quantity bpass_table_path.

Wavelengths are in the star's rest frame, before velocity and Hubble shifts. The finite bounds must satisfy 0 < minimum < maximum and lie inside the native grid. Configuration validation keeps them in double precision without opening the file. Disabled settings are reported and ignored; unknown keys are errors.

Use the non-legacy Gadget stellar loader (PartType4, or PartType1 where stars are stored there). It requires formation scale factors, total metallicity, initial masses and stellar Velocities; missing initial mass is an error. Wind particles are excluded. Formation times rounded just past the snapshot at their stored precision give age zero; larger overshoots are rejected. Ages require a flat matter-plus-Λ cosmology: abs(Omega0 + OmegaLambda - 1) <= 1e-6. The loader converts masses and physical peculiar velocities, including Gadget's square-root-of-scale-factor convention. For TNG at z = 0, explicitly set mass_has_h_factor: true (its z = 0 default is false).

The reader loads only the band and its neighboring samples. The sampler clips boundary intervals to the exact limits and shares template CDFs across stars. It interpolates linear luminosity with weights from log age and log metallicity, then selects a neighboring template in proportion to interpolation weight × photon rate. Per-star caches retain float64 physical rates through catalog cuts; band luminosities use the generator's 1e42 erg/s units.

Packets sample photon number, proportional to λ Lλ. The allocator's luminosity weight is multiplied by \(h\nu_0 / (L/Q)\) to obtain the stored photon-rate weight, where \(\nu_0\) is the reference line frequency, L is the band luminosity and Q is the photon rate. With N packets from a source, each represents Q/N photons per second, independent of its sampled wavelength.

Energy outputs multiply the stored weight by \(\nu_\mathrm{recorded}/\nu_0\), as described in Packet weights. At emission in the source rest frame, this gives the allocator's luminosity weight multiplied by E / (L/Q), where E = hν is the sampled photon energy. The resulting luminosities sum to L statistically. Set nphotons_per_unit_lum to control sampling, subject to per-source rounding and limits.

Launch, emission peeling and intrinsic aperture spectra/IFU cubes use stellar velocities; scattering uses gas velocities. Stellar emission and intrinsic records are completed on the source rank before MPI handoff. Received packets resume transport without stellar lookups. Emission Doppler shifts use each sampled frequency to first order in v/c. SimpleInteractor treats higher Lyman resonances as transparent; LymanSeriesInteractor transports them with cascades (see Lyman-series transport). BPASS sources can opt into the stellar continuum launch shortcut described below.

Stellar continuum launch shortcut

For a stellar continuum near Lyα, an optional preprocessing kernel skips the initial flight to the entry face of the periodic image containing the configured blueward margin. It runs once per fresh emission batch, after the source Doppler shift and emission records have been captured.

mcrtsimulation:
  continuum_advance:
    enabled: true       # default: false
    margin_kms: 2000    # nonnegative blueward velocity offset; default: 2000

This option treats the configured singlesource or listsource population as stellar emission. It supports BPASS continuum and the analytic uniform and gaussian_with_continuum spectra, applying to packets above the margin. It requires Lya, SimpleInteractor, positive finite Hubble flow, and a fully periodic unit box. Spatial termination surfaces and dust interactors are unsupported. A finite, nonnegative margin is required; packets below it use ordinary transport. For u = c × (ν/ν₀ − 1), the target distance in physical units depends on the build's redshift law:

• Default linear (THOR_MCRT_FRACTIONAL_HUBBLE=OFF): s = (u − margin)/H.
• Fractional (THOR_MCRT_FRACTIONAL_HUBBLE=ON): s = (c/H) × ln[(c + u)/(c + margin)].

Here u, margin and c use the same velocity units. The jump stops at the entering face of the periodic image containing that point.

The packet keeps its direction, weight, optical-depth target and stored frequency. The skipped distance is added to lsp_dist, so normal opacity and scattering account for the Hubble shift once. Source records and emission peel seeds remain at the star. Restored packets and MPI relaunches do not trigger the shortcut. Destinations on another rank enter the existing communication queue before local transport. A skip that reaches the incoming-flight distance cutoff is declined.

The skipped flight is assumed transparent: the margin alone does not bound its Lyα wing optical depth. Check convergence against larger margins and a run with the shortcut disabled before using it for production results.