Monte Carlo Radiative Transfer

This page explains the options for Monte Carlo radiative transfer (MCRT) simulations. The MCRT module is designed for spatially and spectrally resolved line transport and is well-suited for modeling complex environments where line photons interact with gas and dust via resonant scattering, absorption, and emission.

Mathematical descriptions of resonant-line transport, the radiation moment estimators, and continuous Doppler transport are collected under Algorithms. Contributors should also consult the frame and unit conventions.

General MCRT Parameters

The following parameters are important MCRT parameters:

• n_photons – Number of photon packets used in the simulation. [Default value: 1]

• n_photons_per_star – Number of photon packets per star particle. [Optional]

• n_photons_max – Maximum number of photon packets allowed when n_photons_per_star is used. [Optional]

• output_stars – Flag to output the emission and escape data for each star. [Default value: false]

• output_cells – Flag to output the emission and escape data for each cell. [Default value: false]

• output_photons – Flag to output escaped and absorbed photon packets. [Default value: true]

• photon_file – Output a separate photon file with this base name, e.g. “photons”. [Optional]

• output_n_scat – Flag to output the number of scattering events for photon packets. [Compile-time parameter | Default value: false]

• output_kick_history – Flag to record the initial state and every scattering kick for all photon packets in a separate rank-local HDF5 file. [Compile-time parameter | Default value: false | Requires: doppler_frequency = true]

  Each MPI rank independently writes <output_base>_kick_history_rankNNNN.<output_ext>. Events are grouped by the OpenMP thread that processed them and stored under /streams/thread_NNNN/events. Each row in the rank-local /photons/metadata table identifies the corresponding thread and half-open event range through event_start and event_count. The first event in each range is the emitted state, and subsequent events record scattering states. A bounded 32,768-row buffer per active thread ensures long histories can be flushed before their photons complete. The /header group describes the stored data layout, reconstruction rules, units, and run configuration. This output can be much larger than typical photon output, so enable it only for convergence studies with an appropriate packet count and storage allocation.

• output_path_length – Flag to output the cumulative path length for photon packets. [Compile-time parameter | Default value: false | Output units: cm]

• output_freq_length – Flag to output the cumulative frequency path length for photon packets. [Compile-time parameter | Default value: false]

• output_collisions – Flag to output collisional excitation data for photons and cameras. Collision-only images and cubes retain the physical brightness of that emission component. Collisional escape fractions are normalized to its own luminosity and are zero when that luminosity is zero. [Default value: false]

• output_source_frequency – Flag to output the emission frequency for photon packets. [Default value: false]

• output_source_position – Flag to output the emission position for photon packets. [Default value: false]

Line Parameters

The following parameters are related to line specifications:

• line – Name of the line to be modeled. [Default value: Lyman-alpha | Other options: Balmer-alpha, Balmer-beta, Balmer-gamma, Balmer-delta, Balmer-epsilon, Paschen-alpha, Paschen-beta, Brackett-alpha, CII-1577409, CIII-1907-1909, CIV-1548-1550, NII-6550, NII-6585, NIII-1747-1749, NIV-1486, OI-6302, OI-6366, OII-3727-3730, OIII-1661-1666, OIII-4364, OIII-4933, OIII-4960, OIII-5008, OIII-883564, NeIII-3969, MgII-2796-2803, SiII-1190-1193, SiII-1260, SiII-1304, SiII-1527, SiIII-1206, SiIV-1394-1403, SII-6718, SII-6733]

  The infrared fine-structure lines use configuration names based on their rounded vacuum wavelengths in angstroms: CII-1577409 and OIII-883564. Computing their collisional emission requires collisions: true, solver_based_emission: true, and chianti.hdf5 supplied through data_dir. Alternatively, supply j_line with read_line_emissivity: true to read emissivities from a file. Stellar and gas continuum can be transported around any configured line without computing gas-line emission.

• recoil – Include recoil induced Doppler shifting with line scattering. [Default value: true]

• x_crit – Constant critical frequency for core-skipping. [Compile-time parameter | Default value: 0 | Valid range: x_crit >= 0]

• dynamical_core_skipping – Non-local core-skipping for resonance lines. [Compile-time parameter | Default value: false]

• n_exp_atau – HEALPix exponent for the non-local \(a \tau_0\) core-skipping estimates, resulting in \(12 \times 4^{n_\text{exp,atau}}\) directions. [Default value: 0 | Valid range: n_exp_atau >= 0]

• n_side_atau – HEALPix number of base pixel subdivisions for the non-local \(a \tau_0\) core-skipping estimates, resulting in \(12 \times n_\text{side,atau}^2\) directions. [Overrides n_exp_atau | Default value: 1 | Valid range: n_side_atau >= 1]

• p_dest – Line destruction probability between absorption and re-emission. [Default value: 0 | Valid range: 0 < p_dest < 1]

• doppler_frequency – Output frequency in units of Doppler widths, otherwise as a velocity offset from line center (in units of km/s). Doppler units should only be used for idealized simulations. Frequency-resolved internal diagnostics enabled by track_freq_bin require this option, including internal spectra, radial pressure and force spectra, trapping spectra, scattering-count spectra, and kick histories. [Default value: false]

• T_exit – Exit temperature for standardized Doppler frequency. [Default value: 1e4 | Units: K]

• v_exit – Exit velocity for shifted observer frames. [Default value: (0,0,0) | Units: cm/s]

• exit_wrt_com – Exit with respect to the center of mass velocity frame. [Default value: false]

Note

In what follows the units of “freq” refer to km/s unless using doppler_frequency.

Dust Parameters

The following parameters are related to dust specifications:

• dust_model – Dust model for the line radiative transfer setting the dust opacity, albedo, and anisotropy parameters. [Options: SMC, MW, LAURSEN_SMC, filename]

  A tabulated model requires at least two finite, positive, increasing wavelengths. The line wavelength must lie within the table, including either endpoint.

• kappa_dust – Dust opacity at the wavelength of the line, scaling the absorption coefficient by the dust mass: \(k_\text{dust} = \kappa_\text{dust} \rho_\text{dust}\) (also works for graphite kappa_dust_graphite, silicate kappa_dust_silicate, and PAH kappa_dust_PAH). [Default value: Based on dust_model and line | Units: cm^2/g of dust]

• sigma_dust – Dust cross section (per metallicity) at the wavelength of the line, scaling the absorption coefficient by the hydrogen number density and metallicity: \(k_\text{dust} = \sigma_\text{dust} n_\text{H} Z\). It is generally recommended to use kappa_dust instead. [Default value: Based on dust_model and line | Units: cm^2/Z/hydrogen atom]

• albedo – Dust scattering albedo defined as \(A = k_\text{scat} / (k_\text{abs} + k_\text{scat})\) (also works for graphite albedo_graphite, silicate albedo_silicate, and PAH albedo_PAH). [Default value: Based on dust_model and line | Valid range: 0 < albedo < 1]

• g_dust – Anisotropy parameter \(\langle \mu \rangle\) for dust scattering (also works for graphite g_dust_graphite, silicate g_dust_silicate, and PAH g_dust_PAH). [Default value: Based on dust_model and line | Valid range: -1 < g_dust < 1]

• f_ion – Survival fraction of dust in \({\rm H\,{\small II}}\) regions. [Default value: 1 | Valid range: 0 <= f_ion <= 1]

• f_dust – Fraction of metals locked in dust grains. [Default value: 0 | Valid range: 0 <= f_dust <= 1]

• f_PAH – Fraction of carbonaceous dust in polycyclic aromatic hydrocarbons (PAHs). [Requires graphite_scaled_PAH | Default value: 0.01 | Valid range: 0 <= f_PAH <= 1]

• T_sputter – Thermal sputtering temperature cutoff, such that no dust exists within gas hotter than this value. [Default value: 1e6 | Units: K]

Radiation Source Parameters

The following parameters are related to line radiative transfer options:

• j_exp – Luminosity boosting exponent for power-law emission biasing, i.e. probabilities are rescaled according to \(\propto \mathcal{L}_\text{cell}^{j_\text{exp}}\) such that \(j_\text{exp} < 1\) can better sample lower emissivity regions. [Default value: 1 | Valid range: 0 < j_exp < 1]

• V_exp – Additional volume bias for cell-source sampling. The sampling probability is proportional to \(\mathcal{L}_\text{cell}^{j_\text{exp}} V_\text{cell}^{V_\text{exp}-1}\). A value of 1 leaves the luminosity-based distribution unchanged, while values below 1 favor smaller cells. Packet weights compensate for the changed sampling probabilities. Nonemitting cells retain zero sampling probability, including boundary cells with zero stored volume. [Default value: 1]

Gas Emission Parameters

These parameters control gas emission processes:

• recombinations – Include recombination emission with luminosity \(\mathcal{L}_\text{rec} = h \nu_0 \int \alpha_\text{B,eff}(T) n_e n_\text{HII} \text{d}V\). [Default value: false]

• T_floor_rec – Apply a temperature floor to recombination line-emission probabilities. Total recombination coefficients still use the gas temperature. The floor uses physical gas temperature for every species. Hydrogenic HeII table lookups then reduce the floored temperature to \(T/4\). The built-in PEQ line probabilities retain the nearest fit value outside \(0.004 \leq 10^{-4}T/Z^2 \leq 2\), where \(Z\) is the recombining ion charge. This avoids extrapolating the fitted probability and does not cap the total rate. [Default value: 0 | Units: K]

• collisions – Include collisional excitation emission. For a transition of frequency \(\nu\), the luminosity is \(\mathcal{L}_\text{col} = h \nu \int q_\text{col}(T, n_e) n_e n_i \text{d}V\), where \(n_i\) is the lower-state ion density. Each component of a doublet uses its own transition frequency and excitation coefficient. This flag works for all implemented hydrogen lines and many collisionally excited metal lines. [Default value: false]

• solver_based_emission – Compute collisional line emission from CHIANTI statistical-equilibrium level populations instead of the existing analytic or tabulated rate coefficients. Requires collisions: true and a supported metal line. Leaving this option disabled preserves the existing rate selection. [Default value: false]

• data_dir – Directory containing chianti.hdf5 for solver-based emission. An empty value uses the process working directory. The file is needed only for line components whose emissivities are computed, which can be circumvented by reading from the input file. See Packaging CHIANTI Atomic Tables for the file format and packaging utility. [Default value: “” | Used with solver_based_emission: true]

Solver-based emission supports the configured CII, CIII, CIV, NII, NIII, NIV, OI, OII, OIII, NeIII, and SII lines, provided the requested transitions are present among the retained CHIANTI levels. The retained level counts are set by ion_n_levels in source/vars.cc. Requests for ions without retained levels, missing transitions, or indistinguishable doublet matches report a setup error instead of silently falling back to another rate model.

The current level solver includes electron collisions and spontaneous radiative decays, but not excitation by neutral hydrogen or molecules, or radiative pumping. This distinction is important when interpreting infrared emission from predominantly neutral gas. The transition probabilities used for the emissivity come from the supplied CHIANTI file, while the built-in line constants above set the photon energy and transfer parameters.

The solver uses each cell’s temperature and electron density to determine the level populations, which are self-consistent with the thermal equilibrium solver in the ionization module. For an upper-level fraction \(p_u\) and radiative rate \(A_{ul}\), it returns \(q_\text{col} = p_u A_{ul}/n_e\), so the existing luminosity normalization gives \(h\nu n_i p_u A_{ul} V\). Cells without electrons produce no collisional emission. Both doublet components use the same carrier-ion data, but their own matched transitions and photon energies. Recombination emission and the line-transfer parameters are unchanged. Direct level solves are more expensive than table interpolation, so this option remains opt-in.

• collisions_limiter – Limit collisional excitation emission by the photoheating rate. Note: Requires collisions to be true and reads the photoheating rate from the initial conditions file. [Default value: None | Valid range: > 1]

• emission_n_min – Minimum number density of emitting cells. [Default value: 0 | Units: cm^-3]

• emission_n_max – Maximum number density of emitting cells. [Default value: infinity | Units: cm^-3]

• emission_ne_min – Minimum electron number density of emitting cells. [Default value: 0 | Units: cm^-3]

• emission_ne_max – Maximum electron number density of emitting cells. [Default value: infinity | Units: cm^-3]

Active total- or electron-density limits load the electron fraction even when recombinations and collisions are disabled. A configured electron_fraction takes precedence over the input x_e field, which defaults to 1 when absent.

• emission_T_min – Minimum temperature of cells eligible for gas-line emission. Cells at this limit are included. [Default value: 0 | Units: K]

• emission_T_max – Maximum temperature of cells eligible for gas-line emission. Cells at this limit are included. [Default value: infinity | Units: K]

• cell_center_emission – Use the cell center for emission instead of uniformly sampling the cell volume. [Default value: false]

Stellar Source Parameters

These parameters define the stellar sources in the simulation:

• continuum – Include stellar continuum emission around the line. This option cannot be combined with unresolved_HII. [Default value: false]

• free_free, free_bound, two_photon – Independently add gas continuum from H/He bremsstrahlung, radiative recombination, or HI/HeII two-photon decay. These options use the same coefficients as the ionization module, evaluated once per cell at the selected line wavelength. They act independently from stellar continuum and gas-line emission. [Default value: false for each process]

• helium_ions – Include helium in the gas continuum. Enabling a gas continuum process requests electron and hydrogen ion fractions, and helium additionally requires the helium abundance and HeI/HeII fractions. HeIII is inferred as the remaining helium fraction. [Default value: true | Used with gas continuum]

Gas continuum is approximated as flat in \(L_\lambda\) over continuum_range (velocity offsets in km/s). Its luminosity is the central spectral emissivity multiplied by cell volume and the corresponding wavelength width. Keep this range narrow enough that continuum slopes and recombination edges are unimportant. Free-free and free-bound coefficients use the ionization tables’ temperature limits, and free-bound emission is zero beyond the wavelength coverage of the existing coefficient tables, including these far-infrared lines. No extrapolated high-level recombination tail is added.

Gas continuum photons use the existing cell-source sampling and transport, including packet-based images/cubes and output_mcrt_emission / output_mcrt_attenuation. The sources/L_gas_cont attribute records their band-integrated luminosity in erg/s, and sources/L_cont includes all continuum components. MCRT’s direct output_proj_* line projections cannot be combined with gas continuum; use flat continuum projections for deterministic continuum images without frequency sampling.

File-backed stellar continuum spectra use logarithmic interpolation. Before interpolation, COLT replaces nonpositive or nonfinite spectral samples in the loaded table with its smallest finite positive value. At least one such value must exist. This also applies to monochromatic continuum tables in the escape-continuum module and does not modify the input file.

• continuum_range – Continuum velocity offset range for frequency sampling. [Default value: (-2000,2000) | Units: km/s]

• unresolved_HII – Generate Lyman-alpha, Balmer-alpha, or Balmer-beta emission directly from the ionizing-photon rates of stellar sources, using conversion factors of 0.68, 0.45, and 0.12, respectively. For example, the Lyman-alpha luminosity is \(\mathcal{L}_{\text{Ly}\alpha} = 0.68 \, h \nu_0 (1 - f_\text{esc}^\text{HII}) \dot{N}_\text{ion}\). This option cannot be combined with continuum. When cell-based line emission is also enabled, use avoid_SFR: true to avoid double-counting emission from star-forming cells. [Default value: false]

• T_unresolved_HII – Set a common source temperature for the initial line profile of unresolved HII regions. If omitted, each source uses the temperature of its containing gas cell. This parameter changes the emitted line width without changing the gas temperature used during radiative transfer. [Optional | Units: K | Valid range: > 0]

• use_stellar_frame – Use the stellar velocity frame, rather than the containing gas cell’s frame, for stellar continuum and unresolved-HII emission. Enabling this option requires stellar velocities in the initial conditions. It changes only the initial emission frequency, and the gas velocity still determines the subsequent radiative transfer. [Default value: false]

• focus_groups_on_emission, focus_subhalos_on_emission – Center group or subhalo positions and velocities on their selected emission. Both options require stellar velocities when stellar emission is enabled. Their halo centroids include the selected gas and stellar luminosities, with one normalization after both contributions are accumulated. Both doublet components contribute to the selected gas luminosity. Effective halo cell counts use the squared sum of those components in each cell. [Default value: false]

• f_esc_HII – Ionizing escape fraction from unresolved HII regions, e.g. \(\mathcal{L}_{\text{Ly}\alpha} = 0.68 \, h \nu_0 (1 - f_\text{esc}^\text{HII}) \dot{N}_\text{ion}\). [Default value: 0 | Valid range: 0 <= f_esc_HII <= 1]

• f_esc_HII_group – Use a model for the line escape fraction from unresolved HII regions based on group metallicities, e.g. \(L_\text{esc}^\text{HII} = L_\text{int} \exp(-\tau_\text{group})\). [Default value: false]

• f_esc_HII_subhalo – Use a model for the line escape fraction from unresolved HII regions based on subhalo metallicities, e.g. \(L_\text{esc}^\text{HII} = L_\text{int} \exp(-\tau_\text{subhalo})\). [Default value: false]

• f_esc_HII_Z0 – Normalization factor \(Z_0\) for the line escape fraction from unresolved HII regions using a power-law model of \(f_\text{esc}^\text{HII} = \exp(-\tau)\) where \(\tau = (Z/Z_0)^\beta\). [Required if f_esc_HII_group or f_esc_HII_subhalo is true | Valid range: 0 < f_esc_HII_Z0 <= 1]

• f_esc_HII_beta – Power-law index \(\beta\) for the line escape fraction from unresolved HII regions using a power-law model of \(f_\text{esc}^\text{HII} = \exp(-\tau)\) where \(\tau = (Z/Z_0)^\beta\). [Default value: 0 | Valid range: f_esc_HII_beta >= 0]

• source_model – Stellar population model used by continuum and unresolved_HII. These models require the initial stellar masses, metallicities, and ages. Currently, unresolved_HII supports only BPASS-CHAB-100. [Default value: “” | Options: BC03-IMF, BPASS-IMF-135-100, BPASS-CHAB-100]

Point Source Parameters

These parameters define point sources in the simulation:

• point – Point source position in [x y z] format. [Optional | Units: cm]

• x_point, y_point, z_point – Point source \((x, y, z)\) coordinates. [Optional | Units: cm]

• L_point – Point source line luminosity. [Default value: 1 | Units: erg/s]

Plane Source Parameters

These parameters define plane sources in the simulation:

• plane_direction – Direction of the plane source. [Options: ‘+x’, ‘-x’, ‘+y’, ‘-y’, ‘+z’, ‘-z’]

• S_plane_cont – Plane continuum surface flux. [Optional | Units: erg/s/cm^2/angstrom]

• F_plane_cont – Plane continuum source flux. [Optional | Units: erg/s/angstrom]

• L_plane_cont – Plane continuum source luminosity. [Default value: 1 | Units: erg/s]

• plane_beam – Use an ellipsoidal beam instead of a rectangular plane. [Default value: false]

• plane_center_x_bbox, plane_center_y_bbox, plane_center_z_bbox – Center of the plane source in bounding box units. For each axis, the supplied fraction \(f\) weights the lower boundary, so the center is \(f x_\mathrm{min} + (1-f) x_\mathrm{max}\). Thus 0.25 places the center three quarters of the way from the lower to the upper boundary, while the default 0.5 selects the midpoint. [Valid range: 0 < value < 1]

• plane_radius_x_bbox, plane_radius_y_bbox, plane_radius_z_bbox – Radii of the plane source in bounding box units. [Valid range: > 0]

Sphere Source Parameters

These parameters define sphere sources in the simulation:

• sphere_center_cont – Sphere center position (continuum) in [x, y, z] format. [Optional | Units: cm]

• sphere_radius_cont – Sphere radius (continuum). [Optional | Units: cm]

• L_sphere_cont – Sphere continuum source luminosity. [Default value: 1 | Units: erg/s]

Shell Source Parameters

These parameters define sphere sources in the simulation:

• r_shell_line – Shell source radius (line emission). Note: A shell source is a spherical shell with a uniform emissivity and outward initial flux. [Optional | Units: cm]

• L_shell_line – Shell line source luminosity. [Default value: 1 | Units: erg/s]

• r_shell_cont – Shell radius (continuum emission). [Optional | Units: cm]

• L_shell_cont – Shell continuum source luminosity. [Default value: 1 | Units: erg/s]

• shell_blackbody – Set the shell luminosity based on a blackbody spectrum. [Default value: false]

• T_shell – Shell blackbody temperature. [Required if shell_blackbody is true | Units: K]

• r_shock – Shock source radius (continuum emission). Note: A shock source is a spherical shell with a uniform emissivity and isotropic initial flux. [Optional | Units: cm]

• shock_blackbody – Set the shock luminosity based on a blackbody spectrum. [Default value: false]

• T_shock – Shock blackbody temperature. [Required if shock_blackbody is true | Units: K]

• L_shock – Shock source luminosity. [Default value: 1 | Units: erg/s]

Note

Stellar continuum and unresolved-HII sources cannot be combined with plane, sphere, shell, or shock continuum sources.

Observed Output Parameters

This section describes the output options for MCRT simulations, including angle-averages, peel-off cameras, and directional binning parameters.

Angle-averaged Parameters

The following parameters are related to angular-averaged escape options:

• output_flux_avg – Output the angle-averaged flux. [Default value: true]

• output_radial_avg – Output the angle-averaged radial surface brightness in radial_avg/data. [Default value: false | Dimensions: (n_radial_avg_pixels) | Output units: erg/s/cm^2/arcsec^2]

• output_radial_cube_avg – Output the angle-averaged radial spectral profile in radial_cube_avg/data. [Default value: false | Dimensions: (n_radial_cube_avg_pixels, n_radial_cube_avg_bins) | Output units: erg/s/cm^2/arcsec^2/angstrom]

• n_radial_cube_avg_bins – Number of frequency bins in the angle-averaged radial spectral profile. The optional radial_cube_avg_freq_range, radial_cube_avg_freq_min, and radial_cube_avg_freq_max parameters configure its frequency range independently of the main spectrum. [Default value: n_bins]

Camera Parameters

The following parameters are related to cameras options:

• n_bins – Number of frequency bins. A positive freq_bin_width determines the integer bin count instead. The resolved width then spans the full frequency range with that count, and output metadata and normalization use this resolved width. The same rule applies to spectral maps, cubes, slits, and halo spectra. [Default value: 100]

• n_slit_bins – Number of frequency bins for both spectral slits. The optional slit_freq_range, slit_freq_min, and slit_freq_max parameters configure their frequency range independently of the main spectrum. [Default value: n_bins]

• n_cube_bins – Number of frequency bins for data cubes. [Default value: n_bins]

• n_radial_cube_bins – Number of frequency bins for radial data cubes. [Default value: n_bins]

• freq_range – Frequency extrema as an exactly two-element [min max] sequence. All product-specific frequency ranges, including group and subhalo variants, require the same shape. Note: Similar setup for slit_freq_range, cube_freq_range, and radial_cube_freq_range. [Default value: (-1000,1000) | Units: freq]

• freq_min – Frequency range minimum. Note: Similar setup for slit_freq_min, cube_freq_min, and radial_cube_freq_min. [Default value: -1000 | Units: freq]

• freq_max – Frequency range maximum. Note: Similar setup for slit_freq_max, cube_freq_max, and radial_cube_freq_max. [Default value: 1000 | Units: freq]

• output_fluxes – Output spectral fluxes. [Default value: true | Dimensions: (n_cameras, n_bins) | Output units: erg/s/cm^2/angstrom]

• output_images – Output the images/data surface-brightness dataset. The image axes follow the configured camera orientation, with the x-pixel index before the y-pixel index. See the group and dataset reference. [Default value: true | Dimensions: (n_cameras, nx_pixels, ny_pixels) | Output units: erg/s/cm^2/arcsec^2]

• output_slits – Output slits/x and slits/y, one spectral slit along each image axis. In both datasets the spatial-pixel index precedes the frequency-bin index. [Default value: false | Dimensions of each dataset: (n_cameras, n_slit_pixels, n_slit_bins) | Output units: erg/s/cm^2/arcsec^2/angstrom]

• output_cubes – Output spectral data cubes in cubes/data. See the spectral product reference. [Default value: false | Dimensions: (n_cameras, nx_cube_pixels, ny_cube_pixels, n_cube_bins) | Output units: erg/s/cm^2/arcsec^2/angstrom]

• output_images2 – Output statistical second-moment images, e.g. to estimate convergence statistics as the relative error is approximately \(\sqrt{\text{images/data2}} / \text{images/data}\). This dataset accumulates squared photon-packet contributions to approximate the variance of the Monte Carlo estimator. Enabling this option also enables output_images. [Default value: false | Dimensions: (n_cameras, nx_pixels, ny_pixels) | Output units: erg^2/s^2/cm^4/arcsec^4]

• output_radial_images – Output radial surface brightness images in radial_images/data. [Default value: false | Dimensions: (n_cameras, n_radial_pixels) | Output units: erg/s/cm^2/arcsec^2]

• output_radial_cubes – Output radial spectral data cubes in radial_cubes/data. [Default value: false | Dimensions: (n_cameras, n_radial_cube_pixels, n_radial_cube_bins) | Output units: erg/s/cm^2/arcsec^2/angstrom]

• output_radial_images2 – Output statistical second-moment radial images, e.g. to estimate convergence statistics as the relative error is approximately \(\sqrt{\text{radial\_images/data2}} / \text{radial\_images/data}\). This dataset accumulates squared photon-packet contributions to approximate the variance of the Monte Carlo estimator. Enabling this option also enables output_radial_images. [Default value: false | Dimensions: (n_cameras, n_radial_pixels) | Output units: erg^2/s^2/cm^4/arcsec^4]

The following parameters produce intrinsic or attenuation only images. Sampled planar variants are stored as images/data_int and images/data_ext; ray-integrated planar variants use proj_images/data_int and proj_images/data_ext. The corresponding cubes and radial products use the same data_int and data_ext child names. Each product group carries its own geometry:

• output_mcrt_emission – Output intrinsic emission without transport based on MCRT sampling (also works for radial images). [Default value: false | Dimensions: (n_cameras, nx_pixels, ny_pixels) | Output units: erg/s/cm^2/arcsec^2]

• output_mcrt_attenuation – Output attenuated emission without scattering based on MCRT sampling (also works for radial images). [Default value: false | Dimensions: (n_cameras, nx_pixels, ny_pixels) | Output units: erg/s/cm^2/arcsec^2]

• output_proj_emission – Output intrinsic emission without transport with ray-based imaging. [Default value: false | Dimensions: (n_cameras, nx_pixels, ny_pixels) | Output units: erg/s/cm^2/arcsec^2]

• output_proj_attenuation – Output attenuated emission without scattering with ray-based imaging. [Default value: false | Dimensions: (n_cameras, nx_pixels, ny_pixels) | Output units: erg/s/cm^2/arcsec^2]

• output_proj_cube_emission – Output intrinsic emission without transport with ray-based imaging. [Default value: false | Dimensions: (n_cameras, nx_cube_pixels, ny_cube_pixels, n_cube_bins) | Output units: erg/s/cm^2/arcsec^2/angstrom]

• output_proj_cube_attenuation – Output ray-integrated line emission with grey dust extinction in proj_cubes/data_ext. This product excludes resonant-line extinction. [Default value: false | Dimensions: (n_cameras, nx_cube_pixels, ny_cube_pixels, n_cube_bins) | Output units: erg/s/cm^2/arcsec^2/angstrom]

• output_proj_cube_direct – Output unscattered line emission in proj_cubes/data_direct, including dust, enabled resonant-line extinction, and enabled electron scattering. The profile is integrated through the full foreground column within each frequency band. Sobolev configurations are not supported. [Default value: false | Dimensions: (n_cameras, nx_cube_pixels, ny_cube_pixels, n_cube_bins) | Output units: erg/s/cm^2/arcsec^2/angstrom]

• proj_emission_profile – Select delta, gaussian, or voigt emission for deterministic spectral cubes, including unresolved-HII sources. Voigt extinction is independent of this choice. Select gaussian to reproduce the previous emission model. See Camera integration for band integration and Hubble flow and continuous Doppler shifts for Hubble-flow support. [Default value: voigt]

• proj_spectral_rtol – Relative spectral integration tolerance, independent of spatial pixel_rtol. Must be finite and strictly between zero and one. [Default value: 1e-6]

• proj_spectral_atol – Absolute spectral tolerance relative to total unattenuated ray luminosity within the requested frequency window. Must be finite and strictly between zero and one. [Default value: 1e-12]

Note

The intrinsic and attenuation-only images do not require resonant-line scattering. Therefore, the projection outputs provide a noiseless ray-based solution and may be requested without enabling the corresponding packet-sampled output_images or output_cubes product. With output_collisions, collision-only projections are computed alongside the total emission. Requesting these images does not change the emissivity used by the corresponding cubes.

If unresolved_HII is enabled, its stellar emission is co-added into both the intrinsic and attenuation-only images and cubes. The attenuated product traces the line-wavelength dust extinction from each stellar position to the camera. There is intentionally no deterministic projection counterpart to the escaped MCRT image because a meaningful escaped line image must include resonant scattering.

Deterministic cubes integrate normalized Gaussian line profiles analytically between each pair of frequency-bin edges. Bin values are spectral densities in the configured output units, so their frequency integral gives the luminosity within the selected range independently of spectral binning. Gas and secondary-line profiles use their local Doppler widths, while unresolved HII sources also respect T_unresolved_HII. Radial grids also integrate the changing projected gas velocity along each cell segment, including configured spherical velocity gradients. The Doppler transformation follows the packet transport convention. These Gaussian projections omit natural line wings and resonant scattering.

MCRT projection outputs default to adaptive spatial convergence for cell emission. Set adaptive: false to use one nine-point Simpson estimate per pixel instead, or set adaptive: true and configure pixel_rtol explicitly. Adaptive integration does not affect unresolved point sources, which are deposited directly into their projected pixels. See Adaptive pixel quadrature for the cube convergence criterion and Hubble flow and continuous Doppler shifts for continuous cosmological shifts in gas and stellar cube emission when hubble_flow is enabled.

The following parameters produce spatially integrated camera summaries [Dimensions: n_cameras]:

• output_escape_fractions – Output full-camera escape fractions and isotropic-equivalent luminosities in erg/s. When aperture_radius is set, the aperture group adds selected luminosities. Intrinsic L_ints, attenuated L_exts, and collisional _col luminosities are included when those products are enabled. [Default value: true]

• output_freq_avgs – Output full-camera mean frequencies. An active aperture_radius adds selected means under aperture, including intrinsic, attenuated, and collisional variants. Each mean is normalized by the light contributing to its product. [Default value: false | Output units: freq]

• output_freq_stds – Output frequency standard deviations. [Default value: false | Output units: freq]

• output_freq_skews – Output frequency skewnesses. [Default value: false]

• output_freq_kurts – Output frequency kurtoses. [Default value: false]

Global and escaped camera frequency summaries use weighted central moments, with zero for the standard deviation, skewness, and kurtosis of an empty or monochromatic sample. Accumulating relative to a local frequency origin preserves narrow widths when a large common offset would cause cancellation in raw moments. Thread and MPI partitions are merged with their weights and means. The aperture luminosities and frequency summaries include the entire sampled spectrum, independently of freq_min and freq_max. Image, map, and mu products retain normalized raw moments. The aperture group contains integrated luminosities, camera spectra, and frequency summaries. Angular and spatially resolved products retain their full outputs. See Packet frequency moments for definitions and aperture normalization.

The following parameters produce flux-weighted raw frequency moments for each image pixel [Dimensions: (n_cameras, nx_pixels, ny_pixels)]. The dataset of order \(n\) contains \(\langle \nu^n \rangle\), rather than a central moment. Requesting an order automatically enables every lower order needed to calculate it and also enables output_images. For example, output_freq4_images: true produces freq, freq2, freq3, freq4, and data inside the images group.

• output_freq_images – Output flux-weighted mean-frequency images. [Default value: false | Output units: freq]

• output_freq2_images – Output flux-weighted second raw frequency-moment images. [Default value: false | Output units: freq^2]

• output_freq3_images – Output flux-weighted third raw frequency-moment images. [Default value: false | Output units: freq^3]

• output_freq4_images – Output flux-weighted fourth raw frequency-moment images. [Default value: false | Output units: freq^4]

The radial-image options follow the same raw-moment and dependency rules and store freq through freq4 in radial_images, with one value per radial image pixel [Dimensions: (n_cameras, n_radial_pixels)]:

• output_freq_radial_images – Output average frequency radial images. [Default value: false | Output units: freq]

• output_freq2_radial_images – Output frequency^2 radial images. [Default value: false | Output units: freq^2]

• output_freq3_radial_images – Output frequency^3 radial images. [Default value: false | Output units: freq^3]

• output_freq4_radial_images – Output frequency^4 radial images. [Default value: false | Output units: freq^4]

The following parameters are useful for on-the-fly false color imaging:

• output_rgb_images – Output flux-weighted frequency RGB images and automatically enable output_images. The normalized image log intensity can be used as the image opacity to illustrate the corresponding spectral data cube. [Default value: false | Dimensions: (n_cameras, nx_pixels, ny_pixels, 3) | Output units: RGB color]

• output_rgb_radial_images – Output flux-weighted frequency RGB radial images and automatically enable output_radial_images. [Default value: false | Dimensions: (n_cameras, n_radial_pixels, 3) | Output units: RGB color]

• output_rgb_map – Output a flux-weighted frequency RGB HEALPix map and automatically enable output_map. [Default value: false | Dimensions: (12 * n_side_map^2, 3) | Output units: RGB color]

• rgb_freq_range – Frequency extrema in [min max] format. [Default value: (-1000,1000) | Units: freq]

• rgb_freq_min – Frequency range minimum. [Default value: -1000 | Units: freq]

• rgb_freq_max – Frequency range maximum. [Default value: 1000 | Units: freq]

• adjust_camera_frequency – Adjust each camera using offsets read from the freq_avgs HDF5 dataset. The dataset must contain one value per camera and use Doppler units when doppler_frequency is enabled, or km/s otherwise. This option was introduced to center Lyα RGB and other camera outputs on the Hα line center. [Default value: false]

• freq_offset_file – Frequency offset file (optional). [Default value: init_file]

• freq_offset_dir – Frequency offset directory (optional). [Default value: init_dir]

• freq_offset_base – Frequency offset file base (optional). [Default value: init_base]

• freq_offset_ext – Frequency offset file extension (optional). [Default value: hdf5]

HEALPix Binning Parameters

The following parameters are related to line-of-sight HEALPix binning options:

• n_exp_map – HEALPix exponent for the line-of-sight map, resulting in \(12 \times 4^{n_\text{exp,map}}\) directions. Note: Similar setup for n_exp_flux (flux maps), n_exp_radial (radial maps), and n_exp_cube (radial cube maps). [Default value: 2 | Valid range: 0 <= n_exp_map <= 13]

• n_side_map – HEALPix number of base pixel subdivisions for the line-of-sight map, resulting in \(12 \times n_\text{side,map}^2\) directions. Note: Similar setup for n_side_flux (flux maps), n_side_radial (radial maps), and n_side_cube (radial cube maps). [Mutually exclusive with n_exp_map | Default value: 4 | Valid range: 1 <= n_side_map <= 13377]

All HEALPix camera and map resolutions must have a pixel count that fits in a 32-bit signed integer. This includes group and subhalo maps and nonlocal optical-depth directions. The same exponent and subdivision limits apply to their corresponding options. Available memory can impose a lower practical limit.

• n_map_bins – Number of frequency bins for the line-of-sight flux map. Note: Similar setup for n_cube_map_bins. [Default value: 100]

• map_freq_range – Frequency extrema in [min max] format for the line-of-sight flux map. Note: Similar setup for cube_map_freq_range. [Default value: (-1000,1000) | Units: freq]

• map_freq_min – Frequency range minimum for the line-of-sight flux map. Note: Similar setup for cube_map_freq_min. [Default value: -1000 | Units: freq]

• map_freq_max – Frequency range maximum for the line-of-sight flux map. Note: Similar setup for cube_map_freq_max. [Default value: 1000 | Units: freq]

• map_radius – Outer radius of the line-of-sight HEALPix radial map. A bare number is interpreted in centimeters. The value may instead include a cm, pc, kpc, Mpc, cMpc, arcsec, bbox, or Rvir suffix, for example map_radius: 1 kpc. The bbox unit is relative to the distance from the camera center to the nearest bounding-box face, while Rvir is relative to the selected halo virial radius \(R_{\rm vir}\). Observed angular sizes require a cosmological simulation. Legacy aliases are allowed as suffixes, e.g., map_radius_bbox: <value> is equivalent to map_radius: <value> bbox. Note: cube_map_radius accepts the same syntax. [Default value: full simulation domain | Units: cm, or as given by the suffix | Valid range: map_radius > 0]

• n_map_pixels – Number of radial bins in the line-of-sight HEALPix map. A configured map_pixel_width overrides this value. With logarithmic spacing, this specifies the number of annuli outside the central disk, and COLT adds the central bin. Note: n_cube_map_pixels has the same behavior. [Default value: 100 | Valid range: n_map_pixels >= 1]

• map_pixel_width – Width of each radial bin in the line-of-sight HEALPix map. It accepts the same unit suffixes as map_radius; for example, map_pixel_width: 10 pc is equivalent to the legacy key map_pixel_width_pc: 10. The legacy key map_pixel_arcsec also remains accepted. A bare number remains a length in centimeters. Specifying a width adjusts n_map_pixels and the outer radius so that the map contains complete bins. Note: cube_map_pixel_width has the same behavior. [Units: cm, or as given by the suffix | Valid range: map_pixel_width > 0]

• min_map_radius_pc – Minimum radius for logarithmic radial bins in parsecs. A positive value makes bin 0 the disk inside this radius; the remaining bins are evenly spaced in \(\log(r)\) up to map_radius. This option is incompatible with a configured bin width. Note: min_cube_map_radius_pc has the same behavior. [Units: pc | Valid range: min_map_radius_pc >= 0]

• min_map_radius_kpc – Minimum radius for logarithmic radial bins in kiloparsecs, with the same behavior as min_map_radius_pc. Note: min_cube_map_radius_kpc has the same behavior. [Units: kpc | Valid range: min_map_radius_kpc >= 0]

• output_map – Output a line-of-sight HEALPix map of escape fractions in map/data. [Default value: false | Dimensions: (12 * n_side_map^2)]

• output_map2 – Output a statistical second-moment line-of-sight HEALPix map, e.g. to estimate convergence statistics as the relative error is approximately \(\sqrt{\text{map/data2}} / \text{map/data}\). This dataset accumulates squared photon-packet contributions to approximate the variance of the Monte Carlo estimator. Enabling this option also enables output_map. [Default value: false | Dimensions: (12 * n_side_map^2)]

• output_flux_map – Output a line-of-sight HEALPix map of spectral fluxes in flux_map/data. [Default value: false | Dimensions: (12 * n_side_flux^2, n_map_bins) | Output units: erg/s/cm^2/angstrom]

• output_radial_map – Output a line-of-sight HEALPix map of radial surface brightness images in /radial_map/data. [Default value: false | Dimensions: (12 * n_side_radial^2, n_map_pixels) | Output units: erg/s/cm^2/arcsec^2]

• output_cube_map – Output a line-of-sight HEALPix map of radial spectral data cubes in /cube_map/data. [Default value: false | Dimensions: (12 * n_side_cube^2, n_cube_map_pixels, n_cube_map_bins) | Output units: erg/s/cm^2/arcsec^2/angstrom for cosmological runs and erg/s/cm^2/Doppler or erg/s/cm^2/(km/s) otherwise]

Angular and radial geometry, including the radial-bin area arrays, is stored as attributes on each map group. flux_map and cube_map also carry their own complete frequency attributes; units remain on each dataset. See HEALPix radial maps for field names.

Group and subhalo radial outputs require the corresponding output_groups or output_subhalos compile-time option. Because these are observer-frame products, requesting one automatically enables output_groups_obs or output_subhalos_obs. The separate virial- and galaxy-aperture statistics are controlled by the options described under Halo Catalog Outputs. Each profile is centered on its catalog halo, with every halo of the same type using the same radial geometry.

• output_radial_map_groups – Output one HEALPix map of radial surface brightness profiles per group in /group/radial_map/data. Note: Similar configuration to output_radial_map but for groups, and likewise for subhalos in /subhalo/radial_map/data. [Default value: false | Dimensions: (n_groups, 12 * n_side_radial_groups^2, n_map_pixels_groups) | Output units: erg/s/cm^2/arcsec^2 for cosmological runs and erg/s/cm^2 otherwise]

• output_radial_avg_groups – Output one angle-averaged radial surface brightness profile per group in /group/radial_avg/data. Note: Similar configuration to output_radial_avg but for groups, and likewise for subhalos. [Default value: false | Dimensions: (n_groups, n_pixels_groups) | Output units: erg/s/cm^2/arcsec^2 for cosmological runs and erg/s/cm^2 otherwise]

The map direction keys are n_exp_radial_groups and n_side_radial_groups. The radial map keys are map_radius_groups, n_map_pixels_groups, and map_pixel_width_groups; the angle-averaged profile keys are radius_groups, n_pixels_groups, and pixel_width_groups. Replace groups with subhalos for subhalo outputs. The radius and width keys accept the same units as map_radius and remain centimeters when given a bare number. A positive min_map_radius_pc_groups or min_map_radius_kpc_groups selects logarithmic map bins, while min_radius_pc_groups or min_radius_kpc_groups does the same for the angle-averaged profile. In each case, the configured bin count describes the logarithmic annuli outside the central disk.

Legacy keys remain supported. These include map_radius_bbox_groups, map_radius_kpc_groups, map_pixel_arcsec_groups, map_pixel_width_pc_groups, radius_bbox_groups, radius_kpc_groups, pixel_arcsec_groups, and pixel_width_pc_groups, together with their subhalos counterparts. For example, map_radius_kpc_groups: 5 is equivalent to map_radius_groups: 5 kpc.

The following parameters produce flux-weighted raw frequency moments for each HEALPix map pixel [Dimensions: 12 * n_side_map^2]. As with image moments, requesting a higher order automatically enables all lower orders and the underlying output_map.

• output_freq_map – Output the flux-weighted mean-frequency map. [Default value: false | Output units: freq]

• output_freq2_map – Output the flux-weighted second raw frequency-moment map. [Default value: false | Output units: freq^2]

• output_freq3_map – Output the flux-weighted third raw frequency-moment map. [Default value: false | Output units: freq^3]

• output_freq4_map – Output the flux-weighted fourth raw frequency-moment map. [Default value: false | Output units: freq^4]

Group and subhalo HEALPix maps use the same option families with _groups or _subhalos appended, such as output_map2_groups, output_freq_map_groups, and output_freq4_map_subhalos. These variants add a leading halo dimension and require the corresponding output_groups or output_subhalos compile-time option. Frequency-moment variants are observer-frame products, so requesting one also enables output_groups_obs or output_subhalos_obs. Their datasets are written inside /group/map or /subhalo/map with the same child names as the global map group.

Angular Binning Parameters

The following parameters are related to line-of-sight angular cosine (\(\mu=\cos\theta\)) binning options:

• n_mu – Number of directional bins for the line-of-sight \(\mu\) map. [Default value: 100]

• n_mu_flux – Number of directional bins for the spectral \(\mu\) map. Its solid-angle normalization uses this resolution independently of n_mu. [Default value: 100]

• n_mu_bins – Number of frequency bins for the line-of-sight flux \(\mu\) map. [Default value: 100]

• mu_freq_range – Frequency extrema in [min max] format for the line-of-sight flux \(\mu\) map. [Default value: (-1000,1000) | Units: freq]

• mu_freq_min – Frequency range minimum for the line-of-sight flux \(\mu\) map. [Default value: -1000 | Units: freq]

• mu_freq_max – Frequency range maximum for the line-of-sight flux \(\mu\) map. [Default value: 1000 | Units: freq]

• output_mu – Output a line-of-sight \(\mu\) map of escape fractions in mu/data. [Default value: false | Dimensions: n_mu]

• output_mu2 – Output a statistical second-moment line-of-sight \(\mu\) map, e.g. to estimate convergence statistics as the relative error is approximately \(\sqrt{\text{mu/data2}} / \text{mu/data}\). This dataset accumulates squared photon-packet contributions to approximate the variance of the Monte Carlo estimator. Enabling this option also enables output_mu. [Default value: false | Dimensions: n_mu]

• output_flux_mu – Output a line-of-sight \(\mu\) map of spectral fluxes in flux_mu/data. [Default value: false | Dimensions: (n_mu, n_mu_bins) | Output units: erg/s/cm^2/angstrom]

The mu and flux_mu groups each carry their own n_mu attribute, and flux_mu also has its own frequency attributes. See Directional maps.

The following parameters produce flux-weighted raw frequency moments in mu/freq through mu/freq4 for the \(\mu\) map [Dimensions: n_mu]. Requesting a higher order automatically enables all lower orders and output_mu.

• output_freq_mu – Output the flux-weighted mean frequency in each \(\mu\) bin. [Default value: false | Output units: freq]

• output_freq2_mu – Output the flux-weighted second raw frequency moment in each \(\mu\) bin. [Default value: false | Output units: freq^2]

• output_freq3_mu – Output the flux-weighted third raw frequency moment in each \(\mu\) bin. [Default value: false | Output units: freq^3]

• output_freq4_mu – Output the flux-weighted fourth raw frequency moment in each \(\mu\) bin. [Default value: false | Output units: freq^4]

Radiation Moment Parameters

These compile-time options in defines.yaml enable estimators of the internal radiation field. All default to false. Momentum and pressure estimators include continuous absorption and the changing radial basis in spherical and cylindrical coordinates; see Radial radiation moments for the integrals and their absorption-free limits. Field names and units are listed under Field Names and Descriptions.

• output_energy_density – Output the path-based radiation energy density.

• output_internal_spectrum – Output the path-based energy density in each frequency bin and cell, with integrated contributions below and above the frequency window.

• output_acceleration_path – Output the path-based radiation acceleration. Spherical geometry records the radial component, cylindrical geometry records radial and vertical components, and Cartesian-like geometries record the vector components.

• output_acceleration_scat – Output radiation acceleration from scattering kicks and continuous absorption, with the same geometry-dependent components.

• output_force_path, output_force_scat – Output the corresponding radial forces in spherical geometry.

• output_pressure – Output radial pressure in spherical geometry; radial, vertical, and radial-vertical shear pressure in cylindrical geometry; or the six independent tensor components in Cartesian-like geometries.

• output_pressure_x – Output radial pressure resolved by frequency in spherical geometry.

• output_P_u_cum – Output the cumulative pressure-to-energy-density ratio in spherical geometry.

• output_trap_x – Output the frequency-dependent trapping ratio and its integrated tails in spherical geometry.

Internal spectra require doppler_frequency: true and use the local gas-frame frequency in local thermal Doppler widths. The bins cover [freq_min, freq_max). The accompanying u_rad_x_tails, P_rad_x_tails, and trap_x_tails arrays contain the integrated contributions below freq_min and at or above freq_max, in that order. To recover the full frequency integral, multiply the spectral bins by freq_bin_width and add both tails. These estimators assign each path segment to a single frequency bin, before the radial velocity-basis update.

Force and acceleration spectra include contributions outside the frequency range in the nearest edge bin. Jump histograms do the same, with zero jumps included in the first bin. Excursion and cumulative frequency diagnostics use absolute Doppler frequencies and require freq_min <= 0 < freq_max. Excursion maxima beyond the range are included in the last bin. Kick histories also record frequencies outside the configured range.

The frequency-resolved momentum options, force multipliers, and pressure-divergence diagnostics use the same segment integrals. The available _e2 options square individual event contributions, while _cm options form cross moments after accumulating a packet’s complete trajectory.

Radial Flow Parameters

The radial flow data is a 3D probe of the radiation field. The output data is in the radial_flow group with an attribute n_radial_bins and fields radial_edges [cm] (radial bin edges), L_src (emission as a function of source radii), L_esc (escape as a function of source radii), L_abs (absorption as a function of absorption radii), L_pass (first passage by source radii and shell boundaries), and L_flow (net flux by source radii and shell boundaries).

The following parameters are compile-time options and need to be set in the defines.yaml file:

• output_radial_flow – Output radial flow data. [Compile-time parameter | Default value: false]

• output_group_flows – Output group flow data. Requires at least one selected group. [Compile-time parameter | Default value: false]

• output_subhalo_flows – Output subhalo flow data. Requires at least one selected subhalo. [Compile-time parameter | Default value: false]

The associated runtime parameters are:

• radial_linspace_pc – Radial linear-spaced binning in parsecs for radial flow data as [min max n_edges]. [Default value: 1 | Units: pc]

• radial_logspace_pc – Radial log-spaced binning in parsecs for radial flow data as [log_min log_max n_edges]. [Default value: 1 | Units: pc]

• radial_linspace_Rvir – Radial linear-spaced binning in units of the selected halo virial radius \(R_{\rm vir}\) for radial flow data as [min max n_edges]. [Default value: 1 | Units: Virial radius]

• radial_logspace_Rvir – Radial log-spaced binning in units of the selected halo virial radius \(R_{\rm vir}\) for radial flow data as [log_min log_max n_edges]. [Default value: 1 | Units: Virial radius]

• radial_edges_pc – Radial edges in parsecs for radial flow data. At least two entries are required. Note: Similar config for radial_edges_Rvir and radial_edges_cm. [Default value: 1 | Units: pc]

Advanced Parameters

The following parameters are related to advanced options:

• velocity_gradients – Integrate line optical depths with a continuous Doppler shift from the local line-of-sight velocity gradient, using the Gaussian profile contribution. See Continuous Doppler transport for the formulation and its scope. [Compile-time parameter | Default value: false | Requires spherical geometry]

• voigt_velocity_gradients – Include the terms through second order in the Voigt damping parameter in the integrated profile. This also enables velocity_gradients. [Compile-time parameter | Default value: false | Requires spherical geometry]

• two_level_atom – Line radiative transfer incorporates occupation number densities from both lower and upper transitions. If this is false then we assume that bound electrons are always in the ground state. Note: This currently also switches from a Voigt to a Gaussian line profile, and assumes complete redistribution for scattering. [Default value: false]

• spontaneous – Include spontaneous emission in two_level_atom line transfer as \(\mathcal{L}_\text{sp} = h \nu_0 A_{21} \int n_2 \text{d}V\). [Default value: false]

• resonant_scattering – Turn resonant scattering on/off. Note: It is not recommended to override this option unless you know exactly what you are doing.

  With a positive compile-time deuterium_fraction, the built-in isotope data support resonant Lyman-alpha only. Other hydrogen lines require isotope parameters in their line implementation. Nonresonant lines without these parameters have zero deuterium opacity. [Default value: Based on line]

• voigt_method – Select the packet Voigt evaluator: second_order uses the expansion through second order in the damping parameter, while faddeeva evaluates the full Voigt profile. Selection occurs once during initialization. It applies to primary and secondary resonant-line opacities and packet camera attenuation. The two-level-atom mode retains its Gaussian profile, and disabling resonant scattering retains zero line opacity. Faddeeva is incompatible with velocity_gradients when a Voigt profile is active because that solver uses separate approximate profile integrals. Hubble-flow builds without that solver support both choices. This option does not change deterministic cube profiles or scattering redistribution. [Default value: second_order]

• anisotropic_scattering – Flag to turn anisotropic scattering on or off (in the core of the line only). [Default value: true for Lyman-alpha, false for other lines]

• load_balancing – Use MPI load balancing algorithms. Note: Due to the OpenMP load balancing this is not recommended in general as the number of MPI tasks is usually much smaller than the number of photon packets. [Default value: false]

Halo spectral controls

With output_groups or output_subhalos compiled, output_flux_avg_groups and output_flux_avg_subhalos request angle-averaged halo spectra. output_flux_map_groups and output_flux_map_subhalos request their HEALPix-resolved equivalents. These flags default to false and enable the corresponding observer escape statistics. output_radial_avg_subhalos and output_radial_map_subhalos request the subhalo radial products described under Halo Catalog Outputs.

Angle-averaged group spectra accept n_bins_groups and either freq_range_groups: [minimum, maximum] or freq_min_groups / freq_max_groups. Subhalos use n_bins_subhalos, freq_range_subhalos, freq_min_subhalos, and freq_max_subhalos. Each defaults to the global bin count and frequency limits. All counts must be positive, and ranges require exactly two increasing values in the configured frequency units.

Flux maps instead use n_map_bins_groups, map_freq_range_groups, map_freq_min_groups, and map_freq_max_groups, with the corresponding _subhalos forms for subhalo maps. Their defaults also come from the global n_bins, freq_min, and freq_max. Angular resolution is set independently by n_side_flux_groups / n_exp_flux_groups or n_side_flux_subhalos / n_exp_flux_subhalos. Changing angular resolution does not change the angle-integrated normalization.

Lyman-alpha transfer through a static slab

This calculation transports line-center photons through a single-cell static neutral hydrogen slab with \(a\tau_0=10^4\) at \(T=100\) K. Here \(a\) is the Voigt damping parameter and \(\tau_0\) is the line-center optical depth from the midplane to either surface. Photons escape after diffusing in frequency, producing a symmetric double-peaked spectrum that can be compared with the analytic static-slab solution.

Escaping packet spectrum and the analytic static-slab solution. Monte Carlo sampling produces small differences between the two peaks.

Choose the optical depth and source

The files in examples/lya_double_peak_slab set the slab half-width to 1 kpc. For the chosen temperature and target \(a\tau_0\), the generator calculates the required neutral hydrogen density:

DnuD = nu0 * np.sqrt(2. * kB * T / mH) / c  # Doppler width [Hz]
a = 0.5 * dnu_L / DnuD                # Voigt damping parameter
sigma0 = np.sqrt(np.pi / 2.) * f12 * ee**2 * np.sqrt(mH / (kB * T)) / (me * nu0)  # Cross section [cm^2]
tau0 = atau0 / a                      # Line center optical depth
n_HI = tau0 / (sigma0 * half_length)  # Neutral hydrogen density [cm^-3]
rho = mH * n_HI / (X * (1. - Z))      # Gas density [g/cm^3]

The generator writes one gas cell between -half_length and half_length, with zero velocity and neutral fraction one. One cell is sufficient because the gas is uniform. The slab geometry extends indefinitely in the transverse directions.

The build settings select this geometry and enable dynamical core-skipping to accelerate diffusion through the line core:

# Compile-time definitions

GEOMETRY: slab                 # Mesh geometry
dynamical_core_skipping: true  # Turn on dynamical core skipping
optimize_temperature: true     # Optimize constant temperature simulations
NO_MPI: true  # Single-rank example with OpenMP

The runtime configuration selects Lyman-alpha, fixes the temperature at 100 K, and disables dust and recoil to match the analytic solution. A point source at the midplane emits at line center:

# Central point source at line center
point: [0, 0, 0]                     # Point source position [cm]
L_point: 1.0e40                      # Point source luminosity [erg/s]

n_photons: 20000                     # Number of photon packets
n_bins: 60                           # Number of frequency bins
freq_range: [-60, 60]                # Frequency range [Doppler widths]
output_photons: true                 # Output escaped photon packets

doppler_frequency: true expresses output frequencies in thermal Doppler widths. T_exit: 100 uses the same temperature for this conversion as the slab gas. The complete configuration sets these options. The packet count controls Monte Carlo sampling, while the number of gas cells controls spatial resolution. They are independent choices.

Run the transfer calculation

From the repository root, generate the slab, copy the configuration, compile the slab executable, and run it with four OpenMP threads:

example_dir="$(pwd)/examples/lya_double_peak_slab/run"
python3 examples/lya_double_peak_slab/generate.py "$example_dir/ics/lya_static_slab.hdf5"
cp examples/lya_double_peak_slab/config.yaml "$example_dir/config.yaml"
make DEFS=examples/lya_double_peak_slab/defines.yaml BUILD="$example_dir/build" EXE="$example_dir/colt" -j4
(cd "$example_dir" && OMP_NUM_THREADS=4 ./colt config.yaml)

Plot the escaped spectrum

output_photons: true writes the escaped packet frequencies and weights to run/output/lya_double_peak_slab.hdf5. The plotting script scales the frequencies by \((a\tau_0)^{1/3}\) and forms a weighted histogram:

x = frequency / atau0_13
bins = np.linspace(-2., 2., 61)
flux, _ = np.histogram(x, bins=bins, weights=weight, density=True)

The analytic curve is normalized in the same frequency coordinate. Packets are not mirrored about line center, so the figure preserves differences between the two peaks from finite sampling. Increasing n_photons reduces those differences. Run the plot with:

python3 examples/lya_double_peak_slab/plot.py "$example_dir" "$example_dir/lya_slab.png"

See Output Files for the packet datasets and Scientific Regression Checks for the independent comparison with the analytic solution.

Static emission and frequency bins

A spectral bin contains the average of a line profile over a finite frequency interval. This example compares COLT’s band averages with a static Gaussian profile, using uniform gas emissivity and no dust or resonant scattering.

Set the emission and frequency window

The generator in examples/frequency_binning writes a Cartesian box with eight cells, a temperature of \(10^4\) K, and line emissivity \(j=10^{-25}\) erg/s/cm³:

with h5py.File(path, 'w') as f:
    f.attrs.update(r_box=1e18, nx=2, ny=2, nz=2)
    for name, value, units in (('rho', 1e-24, 'g/cm^3'), ('T', 1e4, 'K'),
                               ('j_line', 1e-25, 'erg/s/cm^3')):
        f[name] = np.full(8, value)
        f[name].attrs['units'] = units.encode()

read_line_emissivity: true uses this supplied j_line field. resonant_scattering: false and zero dust leave only the emission profile. The deterministic cube settings select a Gaussian, two opposing cameras, and 16 frequency bands:

output_proj_cube_emission: true
proj_emission_profile: gaussian
doppler_frequency: true
cameras: [[1, 0, 0], [-1, 0, 0]]
cube_radius: 5e17          # Aperture half-width [cm]
proj_radius: 1e18          # Projection half-depth [cm]
n_cube_pixels: 4
n_cube_bins: 16
cube_freq_range: [-2, 2]
adaptive: false

cube_radius sets the image half-width, and proj_radius sets the integration half-depth. The aperture is smaller than the box, so every pixel crosses the same emitting depth. doppler_frequency: true and T_exit: 1e4 use the gas thermal width as the frequency unit. Although the configuration also sets a small Monte Carlo packet count, the deterministic cube is calculated by ray integration and does not depend on that count.

Generate the input and run COLT from the repository root:

example_dir="$(pwd)/examples/frequency_binning/run"
python3 examples/frequency_binning/generate.py "$example_dir"
cp examples/frequency_binning/config.yaml "$example_dir/config.yaml"
make DEFS=examples/frequency_binning/defines.yaml BUILD="$example_dir/build" EXE="$example_dir/colt" -j4
(cd "$example_dir" && OMP_NUM_THREADS=4 ./colt config.yaml)

Read the band averages

The plotting script selects one camera, averages its spatial pixels, and constructs the band edges from the saved frequency range:

with h5py.File(directory / 'output/cube.hdf5', 'r') as f:
    g = f['proj_cubes']
    data = g['data_int'][0].mean(axis=(0, 1))
    edges = np.linspace(g.attrs['freq_min'], g.attrs['freq_max'], len(data) + 1)

The unscaled values in proj_cubes/data_int have shape (2, 4, 4, 16) and units erg/s/cm^2/Doppler. For uniform emissivity and path length \(2R\), the expected band average is

\[\overline I_b=\frac{Rj}{x_{b+1}-x_b}\left[\operatorname{erf}(x_{b+1})-\operatorname{erf}(x_b)\right].\]

Here \(R=10^{18}\) cm. Summing the band averages times their widths gives \(2Rj\,\operatorname{erf}(2)\), retaining about 99.53% of the full line integral. The missing fraction lies outside the requested frequency window. Increasing the number of bins resolves the profile more finely but does not recover that missing emission without also widening the window.

python3 examples/frequency_binning/plot.py "$example_dir" "$example_dir/gaussian_bands.png"

The curve and steps use the same analytic peak normalization. The highest bin is slightly below one because it averages over a finite interval. Dotted lines mark the requested frequency window. The display normalization leaves the physical HDF5 values unchanged.

Integrate the last dataset axis with freq_bin_width and the image axes with pixel_area. See Spectral cubes and slits for the output layout and Camera integration for band transfer.

P-Cygni profile from an expanding wind

An expanding shell absorbs continuum on the blue side of a line and scatters light into the red side, producing a P-Cygni profile. This example uses a spherical grid with a central photospheric source, homologous velocity \(v(r)=r/t\), and Sobolev optical depth one. The input and analytic comparison follow the existing P-Cygni benchmark.

The generator in examples/p_cygni_sphere places gas between the photospheric and outer radii. The radial velocity increases with radius, and the first cell contains an ionized emitting core with no neutral hydrogen opacity:

# tau_S = pi f e^2 lambda n_HI / (m_e c dv/dr), with dv/dr = 1/time.
n_HI = me * c / (np.pi * 0.41641 * charge**2 * 1215.6701e-8 * time)
edges = np.r_[0., np.linspace(r_ph, r_out, n_cells)]
velocity = np.zeros((n_cells, 3))
velocity[:, 0] = 0.5 * (edges[:-1] + edges[1:]) / time
neutral = np.ones(n_cells)
neutral[0] = 0.                         # The emitting core contributes no line opacity

The build selects spherical geometry, spherical escape, and the Sobolev approximation:

# Compile-time definitions for the Sobolev expanding wind
GEOMETRY: spherical      # Geometry type: slab, spherical, cylindrical, cartesian, octree, voronoi
sobolev: true            # Adopt the Sobolev approximation for line opacities
spherical_escape: true   # Photons escape from a sphere
NO_MPI: true             # One process with OpenMP

The continuum source is a shell at the photospheric radius. Its luminosity is spread over the specified velocity interval:

# Expanding-shell continuum source
line: Lyman-alpha            # Name of the resonant line
continuum_range: [-5000.0, 5000.0] # Continuum velocity-offset range [km/s]
r_shell_cont: 3.496779230112e12 # Shell continuum source radius [cm]
L_shell_cont: 2e39           # Shell continuum source luminosity [erg/s]

The outer escape radius and camera aperture enclose the wind. The output frequency interval matches the injected continuum interval, and 100,000 packets sample the spectrum in 50 bins. Generate the initial conditions and run the example from the repository root:

example_dir="$(pwd)/examples/p_cygni_sphere/run"
python3 examples/p_cygni_sphere/generate.py "$example_dir/ics/p_cygni_wind.hdf5"
cp examples/p_cygni_sphere/config.yaml "$example_dir/config.yaml"
make DEFS=examples/p_cygni_sphere/defines.yaml BUILD="$example_dir/build" EXE="$example_dir/colt" -j4
(cd "$example_dir" && OMP_NUM_THREADS=4 ./colt config.yaml)

The plotting script reads flux_avg or fluxes and divides by the injected continuum luminosity per velocity interval, calculated from the saved source metadata:

sources = output["sources"]
continuum_luminosity = float(sources.attrs["L_shell_cont"])
continuum_min = float(sources.attrs["Dv_cont_min"])
continuum_max = float(sources.attrs["Dv_cont_max"])
continuum = continuum_luminosity / (continuum_max - continuum_min)

This normalization does not depend on noisy samples in the spectrum’s wings. Plot the result with:

python3 examples/p_cygni_sphere/plot.py "$example_dir/output/p_cygni_sphere.hdf5" "$example_dir/p_cygni.png"

The steps show the Monte Carlo spectrum, so some sampling noise is expected. The curve is a semi-analytic Sobolev solution integrated over impact parameter. Increasing n_photons reduces sampling noise without changing the model. See Spectral cubes and slits for the spectral outputs.