Kratos-linerad: GPU-accelerated Monte Carlo radiative transfer of lines with efficient imaging

Kratos-linerad is a GPU-accelerated Monte Carlo radiative transfer code for spectral lines. The code iterates the level populations to statistical equilibrium alongside the escaping spectrum, treating angle-dependent partial frequency redistribution through constant-memory sampling tables, and it extends the two-step imaging scheme to line transfer, similar to the one previously applied to polarized continuum radiative transfer by H. Yang & L. Wang (2025). This method decouples the Monte Carlo sampling of the scattering physics from the imaging geometry by accumulating the scattering emissivity of scattering resolved in Doppler velocity channels during the Monte Carlo pass. A subsequent deterministic ray-tracing pass then synthesizes channel maps at a fixed viewing angle, getting rid of the excessive wastes of photons in direct imaging scheme based on photon counting.

The code treats a spectral line with rest-frame frequency ν0, upper and lower level degeneracies gu and gl, and spontaneous emission rate Aul. Frequency shifts from line center are measured in Doppler units, x ≡ ∆ν/νD = −vlos/b, where b is the thermal Doppler width. The line profile is the Voigt function H(a, x), with damping parameter a = Aul/(4πνD). The line-center cross section per lower-level molecule is σ0 = (gu/gl)(Aul c3)/(8π3/2ν03b), so that the line-center inverse mean free path, or the extinction coefficient due to scattering, is αsca0 = nlσ0.

For frequency redistribution, the code uses the angle-dependent partial frequency redistribution function RIIA (G. B. Rybicki & D. G. Hummer 1992), which correlates the outgoing frequency and direction with the incoming direction and captures the coherent core and the frequency-randomizing wing branches. Rather than evaluating the four-dimensional redistribution kernel on the fly, Kratos-linerad precomputes a sampling table for the conditional distribution P(u∥ xin) of the velocity component u∥ of the scattering atom along the incoming photon direction, whose unnormalized density is the product of the thermal Maxwellian and the atomic Lorentz profile. This tabulated sampler is referred to as the USampler. The table is built once at initialization: at each node of an incoming-frequency grid, comprising eighteen uniformly spaced points below xin = 8 and logarithmically spaced points extending to xin = 300, the density is evaluated on a uniform grid of u∥ ∈ [−6, 6] and is cumulatively summed into a normalized cumulative distribution, which is stored as its logarithm. At propagation time a packet draws u∥ by inversion of the cumulative distribution, with a binary search on the incoming-frequency grid and linear interpolation between adjacent rows; the perpendicular velocity component is drawn from a Gaussian, and the outgoing frequency is reconstructed from the two components and the scattering angle.

The same building procedure also tabulates the three-dimensional RIIA kernel used for the imaging weights. The kernel is the Gaussian-mixture form of the D. G. Hummer (1962) redistribution function. The kernel is pre-tabulated rather than evaluated on the fly, using the discrete form R(∆; xin, g) = Σk pk G[∆ − uk(g − 1); 2−1/2 sin γ], written in terms of the frequency kick ∆ ≡ xout − xin, where pk denote the discrete probabilities of the velocity sampler. The kick distribution is localized, ∆ ≲ 10, so the table spans only ∆ ∈ [−10, 10] at a uniform spacing of 0.1, together with xin ∈ [0, 120] and g ∈ [−1, 1], giving 200 × 200 × 40 nodes. The symmetry R(∆; −xin, g) = R(−∆; xin, g) halves the tabulated incoming-frequency range, lookups are trilinear, and the table returns zero for ∆ > 10. Beyond the tabulated range, xin ≥ 120, the conditional velocity distribution has already converged to the thermal Gaussian, and the kernel is evaluated directly from its asymptotic expression. The tabulation preserves the normalization to better than 0.1 per cent everywhere; in total variation it reproduces broad kernels (g ≲ 0.7) to within 2.5 per cent and near-forward kernels to within 1 per cent.

For statistical equilibrium, the level populations are determined by balancing radiative and collisional transitions. For a two-level system the population ratio is nu/nl = (Γ + klu ncoll)/(Aul + (gl/gu)Γ + kul ncoll), where Γ is the radiative pumping rate, kul and klu are the collisional de-excitation and excitation rate coefficients, and ncoll is the collider density. The collisional destruction probability ϵ = kul ncoll/(Aul + kul ncoll) is added to the absorption component. For multi-level systems Kratos-linerad assembles the full rate matrix from LAMDA-format data and solves it directly. The radiation field and the populations are made self-consistent through a Lambda iteration, beginning from collisional equilibrium and propagating an initial photon population; the accumulated excitation rate updates the populations, which in turn rescale the line opacities for the next propagation. To avoid double-counting the source of radiation, the initial emission photons are frozen across cycles, and only the scattering opacity is updated from the converged lower-level population in each cycle.

During the Monte Carlo pass, each path integration segment contributes to a velocity-resolved scattering emissivity field, recorded on each cell on a mesh. Towards the fixed camera direction ncam on the kth velocity channel in a cell index i, the effective scattering emissivity reads jsca,k(i) = Σpp (Fpp/4π) αsca,0 H(a, xin) RIIA(xout; xin, g) (1 − e−δτe)/δτe, where the sum runs over photon packets traversing the cell, Fpp is the equivalent photon number flux in the current cell calculated with the proper photon weight of packet, αsca,0 the scattering extinction coefficient at the line center, and δτe the extinction optical depth of the segment. The channels vk ≡ vmin + (k + 1/2)∆v bin the escaped photons by Doppler velocity. In media with a bulk velocity field, photon packets are sampled in velocity space so that the Doppler shift vobs = n · vbulk is projected onto each channel.

The propagation offers three photon modes. ph mode= 1 and ph mode= 2 both implement the exact redistribution described above and differ only in where the tables reside: ph mode= 1 keeps the sampling table and the two-dimensional Voigt table in global device memory, whereas ph mode= 2 extracts the one-dimensional slice at the simulation’s damping parameter, retabulated on a one-dimensional logarithmic frequency points, and places it in GPU constant memory together with the constant-memory USampler form. ph mode= 3 retains the exact redistribution sampling but replaces the tabulated Voigt opacity profile with an analytic approximation that blends the Doppler core and the damping wing with a hyperbolic-tangent switch; this mode is offered as a demonstration for the necessity of consistent Voigt profile, and is not for production use.

Following the Monte Carlo sampling, with the effective scattering emissivity ready, Kratos-linerad performs a deterministic ray-tracing pass along lines of sight toward the camera, integrating the line radiative transfer equation cell by cell. Within a cell, the formal solution reads, in each Doppler velocity channel indexed k, Iout,k = Iin,k exp(−δτt,k) + Sk[1 − exp(−δτt,k)], where δτt,k = αt,k δl is the total extinction optical depth across the ray segment of length δl through the cell, and the total opacity combines the scattering and absorption contributions, αt,k = αsca,0 H(a, xk) + αabs; Sk = jsca,k/αt,k + Semis,k, where αt,k is the total extinction coefficient in the kth channel, Sk the effective scattering source function, and Semis,k is the channel-resolved emission source function that sums the contribution from the line transition itself, and other possible sources including thermal emission. Such imaging methods enable significantly higher utilization of photons compared to direct photon-counting imaging: the ray tracing step is based on the already-sampled scattering emissivity and does not require propagating additional photon packets whose majority would be discarded. The imaging cost is therefore decoupled from the scattering optical depth and from the number of Monte Carlo packets, which is precisely what makes the method tractable when combined with the population iteration: this cost is incurred only once, after the population iteration has converged.

Kratos-linerad builds on the heterogeneous computing framework of Kratos and runs the photon propagation on GPUs. The design emphasizes the flexibility of selecting precision, including a full double precision mode, and a mixed-precision mode: floating-point computation is carried out in single precision wherever it does not compromise the accumulated statistics; double precision is reserved for the host-side table construction, the accumulation reductions, and the few operations that require it. The most frequently accessed quantities in the propagation loop are the redistribution sampler and the Voigt profile. Rather than reading these from global memory at every event, Kratos-linerad places the sampling tables in NVIDIA constant memory, which is served by a hardware cache and a broadcast path that is substantially faster than global-memory access when all threads in a warp request the same or nearby addresses. The two-dimensional USampler conditional distribution, with dimensions 251 × 40 and size roughly 40 kB, and a one-dimensional Voigt profile tabulated at 5000 points in logarithmic frequency space both reside in constant memory. This arrangement could yield an access speedup of roughly ≳ 50× relative to an equivalent global-memory lookup. The three-dimensional RIIA kernel table, with dimensions 200 × 200 × 40 and size roughly 6.4 MB, is nonetheless too large for constant memory and resides in global device memory.

The propagation is parallelized over photon packets. In the default server-worker mode, worker threads fetch packets from a global work queue rather than being bound to a fixed packet, which balances the load when the per-packet propagation cost varies strongly across the volume; this improves the throughput by about a factor of two relative to the classic one-thread-per-packet scheduling. The random number generator is a portable generator with per-thread state, which avoids contention. The scattering emissivity, the excitation flux, and the escaping photon flux are accumulated with atomic operations into global memory. The optional imaging ray-tracing pass is parallelized with one thread per image pixel. By design the code does not employ core-skipping, and every photon is propagated through the full optical depth of the medium. The pipeline is organized as a Python orchestrator that manages the population iteration, the field preparation, and the input and output, together with a set of GPU kernels, with the Monte Carlo transport and the imaging implemented as separate kernels. The kernels are compiled through CUDA, but the code is not tied to NVIDIA hardware: it builds equally against AMD’s HIP and can also run on CPUs through HIP-CPU.

A suite of tests validates the proper implementation. The escaped spectrum is first validated against the analytical scaling of D. A. Neufeld (1990). A plane-parallel slab with 128 × 2 × 2 cells, periodic boundaries transverse to the slab normal, a central isotropic source of line-center photons, and a damping parameter a = 0.149, spanning mean optical depths τ0 ∈ 200, 500, 2000, 8000, 32000 is propagated with 105 photon packets. The core statistics measured are the peak locations of the double-peaked profile, whose analytic scaling is xpeak = 0.881(aτ0)1/3. The exact treatments agree with the analytic scaling to within a few percent across the full range of optical depths, and the two memory variants are consistent at the one to ∼ 2% level.

The two-step scheme decouples absorption from scattering: the Monte Carlo pass samples the scattering emissivity without absorption, and the imaging pass applies the absorption in the formal solution. This separation is validated against the dusty-slab test of A. Verhamme et al. (2006), in which a uniform static slab at T = 10 K (a = 0.015) with line-center optical depth τ0 = 105 and a central plane of monochromatic line-center photons is embedded with dust at seven absorption depths. The escape fraction fesc is measured from 2 × 105 photon packets per point. The Kratos-linerad measurements track the analytic reference curve of A. Verhamme et al. (2006) to within ≃ 12 per cent across nearly three decades in fesc, from the nearly unabsorbed limit fesc ≃ 1 down to fesc ∼ 10−3.

The imaging pass is first validated in the optically thin regime, where the formal solution reduces to a direct volume integration of the emissivity. A uniform scattering slab of thickness Lslab and line-center optical depth τ0 = 0.01, with a near-Doppler profile (a = 0.01, b = 1 km s−1), is placed in a cubic domain with free boundaries, illuminated by a plane-parallel beam of line-center photons propagating in the plane of the slab, and viewed face-on by a camera looking along the slab normal, so that the scattering angle between the incoming photon direction and the line of sight is 90◦ (g = 0). The slab moves at vz = b along the line of sight. The measured channel spectrum agrees with the analytic shape to within ≃ 0.4 per cent and the integrated intensity to within ≃ 2 per cent with 2 × 105 photon packets, and the line peaks at v = −vz as expected.

The imaging module is also validated by comparing the peak of the velocity-resolved channel map against the same Neufeld scaling. A camera is placed on the +x face with θ = π/2 and ϕ = 0, using 64 uniformly spaced velocity channels whose range adapts to the expected escape width, spanning mean optical depths τ0 ∈ 300, 1000, 3000, 10000, 30000, 100000 with 105 photon packets per point. The imaging peaks reproduce the Neufeld prediction to within roughly 10% across the tested range of optical depths, with the imaging peak lying slightly inside the escape peak at the lowest optical depths and the two converging at the highest. This offset is physical rather than numerical.

The wall-clock performance of the two stages is summarized as follows. Sampling the scattering emissivity with 105 photon packets and up to 107 propagation steps, with exact redistribution and no core-skipping, completes in roughly ≲ 2 s on a single NVIDIA RTX 3090 GPU. This is about 15× faster than the eight-core CPU reference (AMD Ryzen 7 5800X) running the same propagation. The subsequent ray-tracing imaging stage synthesizes a 256×256-pixel map with ten velocity channels in ≲ 3 s, and this cost is incurred only after the population iteration has converged.

A controlled, decomposed benchmark was constructed against SKIRT 9, operated in its Lyα transfer mode with all acceleration schemes turned off. The setup is a uniform cubic medium on a 323 Cartesian grid at T = 100 K (b = 1.285 km s−1, damping parameter a = 4.73×10−3), with line-center optical depth τ0 = 103 and a uniform volume source emitting a Gaussian line profile of the thermal width. The escaped spectra of the two codes are consistent: both show double peaks at x ≃ 2.4–2.6, matching to 11 per cent. The transport kernel is 6× faster at Nph = 103 and 132× faster at Nph = 106, saturating at 155× at Nph = 107, reflecting the throughput contrast between a single RTX 3090 GPU and a 16-thread CPU. The end-to-end ratio climbs from 0.5× at Nph = 103 to 112× at Nph = 107, the fixed overhead being significant only at the lowest photon counts.

As an application example, the pipeline is applied to a snapshot of the fiducial turbulent diffuse molecular cloud simulation of N. Yue et al. (2024). The simulation follows three-dimensional, non-ideal magnetohydrodynamics coevolved with nonequilibrium thermochemistry (a 23-species network) with Athena++ in a periodic (0.04 pc)3 box on a uniform 1283 grid. The snapshot analyzed here has a mean gas temperature ⟨T⟩ ≃ 47 K (peak ≃ 610 K), mean densities ⟨nCO⟩ ≃ 2.8×10−5 cm−3 and ⟨nOH⟩ ≃ 5.8×10−6 cm−3, and a line-of-sight velocity dispersion σv ≃ 0.5 km s−1. Two transitions are synthesized: the OH 18 cm ground-state Λ-doublet line at 1665.402 MHz, with a mean line-center optical depth ⟨τ⟩ ≃ 0.3 rising to τ ≃ 1 along the densest sightlines, and the CO J = 1 → 0 line at 115.27 GHz, which remains optically thin (τ ≃ 0.05). Both are modeled as two-level systems with LTE level populations at the local gas temperature, exact RIIA redistribution with the constant-memory kernel table, and no continuum opacity. Channel maps are generated with the two-step imaging scheme on 1282 pixels with five velocity channels of width ∆v = 1 km s−1 spanning v ∈ [−2, 2] km s−1, viewed along the x-axis of the simulation box; the Monte Carlo pass used ≃ 2.1 × 106 photon packets per transition. At relatively low optical depths, the velocity-integrated CO and OH intensities correlate with the true column densities with correlation coefficients r ≃ 0.90 and 0.97 in logarithmic space, respectively. The channel maps also recover the cloud kinematics, in which the emission centroid sweeps by ∼ 0.01 pc across the box between the −1 and +1 km s−1 channels. The cost of obtaining these channel maps remains modest: each deck of images, one Monte Carlo pass followed by the deterministic ray-tracing pass through all channels, took ≃ 12 s of wall-clock time per transition on a single NVIDIA RTX 3090 GPU including the IO overheads.

The principal current limitation is the split between the GPU-resident Monte Carlo transport and the Python-side orchestration. The population iteration, the photon generation, and the file IO are handled by the Python pipeline, and their combined overhead can exceed the GPU transport kernel time by one to two orders of magnitude. The absence of core-skipping makes extremely optically thick media, such as resonant lines with τ0 ≳ 106, relatively expensive, while the planned discrete-diffusion Monte Carlo extension addresses this regime. Additionally, the current implementation targets single-species line transfer, and the treatment of multiple interlacing transitions is a natural extension. The two-step architecture is well suited to forward-model fitting: once the scattering emissivity has been sampled, re-synthesizing images for new channel definitions or spatial binning is computationally cheap, while a new viewing geometry requires re-accumulating the camera-directed scattering emissivity in a new Monte Carlo pass.

Improvements to AI Systems:

1. GPU-Accelerated Radiative Transfer for AI Training Data Generation

• Integrate Kratos-linerad as a fast, GPU-based forward model to generate synthetic spectral line observations (e.g., CO, OH, Lyα) for training AI models in astrophysical inference.

• The improved AI system can rapidly produce large, physically accurate datasets of channel maps and spectra across diverse optical depths, velocity fields, and damping parameters, enabling robust training of deep learning models for tasks like gas column density mapping, kinematic reconstruction, and line profile classification.

2. Hybrid AI-Physics Surrogate Modeling

• Use Kratos-linerad’s two-step imaging (Monte Carlo scattering emissivity + deterministic ray tracing) to decouple expensive scattering sampling from cheap imaging, allowing AI surrogates to learn the scattering emissivity field directly from input conditions (density, temperature, velocity).

• The improved AI system can predict emissivity fields and then apply the deterministic ray-tracing step on-the-fly, reducing inference time by orders of magnitude while retaining physical accuracy for arbitrary viewing angles and velocity channels.

3. AI-Driven Adaptive Sampling and Optimization

• Leverage the code’s constant-memory tabulated samplers (USampler, RIIA kernel) to train AI models that predict optimal photon packet distributions or adaptive frequency grids for Monte Carlo transport, minimizing variance and computational cost.

• The improved AI system can dynamically allocate computational resources (e.g., photon counts, grid resolution) based on predicted optical depth or scattering complexity, accelerating simulations for high-τ regimes without sacrificing accuracy.

4. Real-Time Bayesian Inference and Parameter Estimation

• Combine Kratos-linerad’s fast GPU transport (≤2 s for 105 packets) with AI-based emulators to enable real-time Markov Chain Monte Carlo (MCMC) or nested sampling for fitting observed line profiles.

• The improved AI system can infer physical parameters (temperature, density, velocity dispersion, dust content) from multi-line observations in seconds, rather than hours, by replacing full radiative transfer calls with AI-predicted spectra trained on Kratos-linerad outputs.

5. Multi-Species and Multi-Transition Analysis

• Extend the code’s single-species capability using AI to model interlacing transitions (e.g., CO + OH + [CII]) by learning cross-species coupling effects from Kratos-linerad simulations.

• The improved AI system can jointly analyze multiple spectral lines to disentangle temperature, density, and chemical abundance gradients in turbulent clouds, improving 3D tomographic reconstructions of interstellar medium.

6. AI-Accelerated Population Iteration

• Replace the Lambda iteration for statistical equilibrium with a neural network that predicts converged level populations from local conditions (density, radiation field, collider rates), trained on Kratos-linerad’s exact rate-matrix solutions.

• The improved AI system can skip iterative cycles, directly outputting populations for multi-level systems, reducing computational cost by 10–100× for complex molecules and enabling large-scale parameter surveys.

7. Automatic Core-Skipping and Photon Management

• Train a reinforcement learning agent using Kratos-linerad’s performance metrics to decide when to apply core-skipping or discrete-diffusion Monte Carlo for optically thick media (τ0 ≳ 106).

• The improved AI system can autonomously switch between transport algorithms based on local optical depth, balancing accuracy and speed, and extending feasibility to extremely thick resonant lines.

8. AI-Based Image Synthesis and Super-Resolution

• Use Kratos-linerad’s channel maps as ground truth to train generative models (e.g., diffusion or GANs) that produce high-resolution synthetic observations from low-resolution or partial emissivity data.

• The improved AI system can upscale 1282 pixel maps to 10242 with physical consistency, or generate new viewing angles from a single Monte Carlo pass, enabling rapid exploration of observational geometries for telescope planning.

9. Uncertainty Quantification for AI Predictions

• Use Kratos-linerad’s exact solutions to calibrate AI models’ predictive uncertainty, particularly in optically thick or multi-scattering regimes where approximations fail.

• The improved AI system can output confidence intervals for inferred parameters, flagging regions where the AI emulator deviates from full radiative transfer, thus guiding follow-up high-fidelity simulations.

10. Portable AI-Hardware Co-Design

• Leverage Kratos-linerad’s HIP-CPU portability to train AI models that optimize kernel execution across GPUs (NVIDIA, AMD) and CPUs, learning hardware-specific scheduling policies.

• The improved AI system can automatically select precision modes (single vs. double) and memory layouts (constant vs. global) to maximize throughput on any given hardware, reducing energy consumption and wall-clock time for large surveys.

Spectral lines encode the velocity, the temperature, and the chemical structure of astrophysical gas; interpreting them requires radiative transfer that is accurate at high optical depth, consistent with the local excitation, and efficient enough to synthesize velocity-resolved images. We present Kratos-linerad, a GPU-accelerated Monte Carlo radiative transfer code for spectral lines. The level populations can be iterated to statistical equilibrium together with the escaping photon distribution, treating angle-dependent partial frequency redistribution through constant-memory sampling tables. The code adopts the two-step imaging scheme to line transfer, in which a Monte Carlo pass samples a velocity-resolved scattering emissivity, and a deterministic ray-tracing pass synthesizes channel maps decoupled from the scattering geometry. Validation reproduces the analytic scaling of escaped spectra and the imaging double-peak profiles. GPU parallelism enables both stages sufficiently fast for