Gridless SPH-MCRT Implementation

• Gridless SPH-MCRT is a computational technique that applies Monte Carlo radiative transfer directly to ungridded SPH density fields, preserving adaptive resolution and simulation accuracy.
• It leverages an octree-based acceleration structure to efficiently compute ray–sphere intersections and trace photon packet interactions through SPH smoothing volumes.
• The method synthesizes realistic telescope images and datacubes by mapping photon packet weights with support for scattering, absorption, and thermal equilibrium.

A gridless SPH-MCRT implementation is a computational technique that applies Monte Carlo Radiative Transfer (MCRT) directly to Smoothed Particle Hydrodynamics (SPH) density fields without gridding, thus maintaining the fully adaptive nature of SPH and preserving the original simulation’s accuracy. Photon packets are propagated, absorbed, or scattered within the continuous, particle-based fluid representation, generating synthetic telescope images and datacubes that directly reflect the physical and spatial characteristics of the original SPH dataset (Forgan et al., 2010).

1. Algorithm Structure and Workflow

The gridless SPH-MCRT workflow comprises five fundamental phases:

1. Initialization:
   • SPH simulation snapshots are ingested, providing particle positions $r_j$, masses $m_j$, densities $\rho_j$, smoothing lengths $h_j$, and temperatures $T_j$.
   • An acceleration structure (typically an octree) is constructed over all SPH particles. Each node stores an axis-aligned bounding box (AABB) that encapsulates the smoothing volumes (radius $2h_j$) of all its particles.
   • Emission sources (point sources such as stars and/or diffuse emission from SPH particles) are specified. The total photon packet count $N_{\gamma,\mathrm{tot}}$, viewing angles $(\theta_v,\varphi_v)$, image plane size, and resolution are defined.
2. Photon Packet Emission:
   • The photon budget $N_{\gamma,\mathrm{object}}$ per source is allocated according to $$N_{\gamma,\mathrm{object}} = N_{\gamma,\mathrm{tot}} (L_{\mathrm{object}} / L_{\mathrm{tot}})$$.
   • Each packet’s wavelength $\lambda$ is sampled from the emission spectrum (e.g., blackbody $B_\nu(T)$). The initial direction is sampled isotropically unless otherwise specified. For diffuse emission, packets are emitted from $r_j$ with probability proportional to $m_j \kappa_\nu B_\nu(T_j)/\rho_j$.
3. Packet Propagation and Interaction:
   • Each photon packet is traced along its path until its weight drops below a threshold or it escapes.
   • The interaction optical depth $\tau_\mathrm{int}$ is sampled as $-\ln(1-\xi)$ with $\xi$ uniform in $[0,1]$.
   • Along the ray, all SPH particles intersected (as determined by the octree) are identified and traversed in order.
   • The incremental optical depth $\tau_j = \kappa_\nu \int_{\ell_j} \rho(r) d\ell$ is accumulated for each intersected smoothing sphere $\ell_j$.
   • When the total optical depth meets or exceeds $\tau_\mathrm{int}$, the precise interaction point is located by bisection to solve $$\int \kappa_\nu \rho\,d\ell = \tau_\mathrm{int}-\tau$$.
   • At the interaction, packets are scattered (probability $$a = \sigma_\mathrm{scat}/(\sigma_\mathrm{abs}+\sigma_\mathrm{scat})$$) or absorbed and re-emitted.
4. Image-Plane Capture:
   • Exiting photon packets have their rays intersected with a 2D image plane, using the specified geometry.
   • The spatial intersection is mapped to pixel indices, and the photon weight $w$ is accumulated in intensity arrays $I(x, y, \lambda)$, forming a 2D image or 3D datacube depending on spectral binning.
5. Postprocessing:
   • The final image can be smoothed to simulate telescope point spreads with a Gaussian PSF matching the instrument beam.

This entire workflow operates without transforming SPH data onto any grid, thus eliminating re-gridding errors and capturing the full adaptivity of the original hydrodynamical simulation (Forgan et al., 2010).

2. SPH Density Field Representation

SPH density at any spatial location $r$ is evaluated as a continuous sum over neighboring particles: $\rho(r) = \sum_j m_j W(|r - r_j|, h_j)$ where $W$ is the cubic-spline kernel with compact support: $$W(q) = \begin{cases} (1 - \frac{3}{2}q^2 + \frac{3}{4}q^3)/(\pi h^3), & 0 \le q < 1\ \frac{1}{4}(2-q)^3/(\pi h^3), & 1 \le q < 2\ 0, & q \ge 2 \end{cases}$$ with $q \equiv |r - r_j|/h_j$ and each $h_j$ adaptively set so its sphere of radius $2h_j$ contains $N_\text{neigh} \approx 50 \pm 10$ neighbors. Neighbor search and ray–sphere intersection both use the octree's AABBs for efficient spatial queries. This maintains the full accuracy and adaptivity of the particle distribution for radiative transfer and imaging.

3. Monte Carlo Radiative Transfer in SPH Fields

Photon packet propagation uses MCRT principles applied natively to the SPH field:

• Optical Depth Integration: Along a path $s$, the cumulative optical depth is

$\tau(s) = \int_0^s \kappa(\lambda) \rho(r(t)) dt$

• Interaction Sampling: For each photon, the interaction length $L$ is sampled by solving $\tau(L) = \tau_\mathrm{int}$ where $\tau_\mathrm{int} = -\ln \xi$, $\xi \in (0,1)$.
• Scatter vs. Absorption: The single scattering albedo is $$a(\lambda) = \sigma_\mathrm{scat}/(\sigma_\mathrm{abs}+\sigma_\mathrm{scat})$$, dictating the probabilistic outcome at each interaction.
• Scattering Phase Function: The Henyey–Greenstein function,

$$P(\mu) = \frac{1 - g^2}{4\pi (1 + g^2 - 2g\mu)^{3/2}}$$

with anisotropy parameter $g \in [-1,1]$, models the angular distribution of scattered photons.

• Weight Attenuation: The packet’s weight decays between interactions as $$w_\text{out} = w_\text{in} \exp(-\tau_\text{segment})$$.
• Polarization: Packet Stokes vectors are updated via the local Mueller matrix at each scattering.

For absorption, the possibility of radiative equilibrium (immediate re-emission at the updated local temperature) is supported.

4. Ray Traversal and Optical Depth Calculation

Photon–SPH interaction is governed by precise ray traversal:

• Acceleration Structure Usage: Rays traverse the octree hierarchy using the slab method to find intersecting AABBs; only leaf nodes are further checked for SPH particle intersections (spheres of radius $2h_j$).
• Ray–Sphere Intersection: For each intersected sphere, the segment $\ell_j$ within the kernel support is determined, and the incremental optical depth is

$$\Delta \tau_j = \kappa(\lambda)m_j \int_{\ell_j} W(|r - r_j|, h_j) d\ell$$

The cubic-spline kernel's analytic form allows these integrals to be precomputed or reduced to a 1D lookup.

• Interaction Point Localization: Within the final intersected smoothing volume, the path length to reach $\Delta \tau = \tau_\mathrm{int}$ is found by recursively bisecting the segment until the solution is within tolerance.

Photon packets that exit the system without interaction are projected onto the image plane.

5. Image and Datacube Synthesis

Synthetic images are assembled by accumulating escaping photon packet weights onto a user-defined pixel grid. The final image properties are governed by:

• Image Grid: A grid of size $N_x \times N_y$ encompasses the modeled field of view, resolution matching that of the desired observational setup. Photon impacts binned as $I(x,y,\lambda)$ form either 2D images (broadband) or 3D datacubes (if spectral binning is applied).
• Spectral Binning: If included, photon wavelengths are accumulated in discrete bins to generate spectrally resolved maps.
• Postprocessing: Gaussian point spread functions can be applied to match telescope angular resolution (FWHM equal to the instrument beam).

This enables the generation of synthetic observations that faithfully reproduce the spatial and spectral features of the simulated physical system.

6. Performance Characteristics and Resolution

The computational cost $T$ is approximately

$$T \sim N_\gamma \times \langle N_\text{ray} \rangle \times N_\text{scatt}$$

where $$\langle N_\text{ray} \rangle \simeq N_p P_\text{intersect}$$ is the mean number of smoothing volumes intersected per ray, and $N_\text{scatt} \sim \langle \tau \rangle^2$. $P_\text{intersect}$ for a uniform sphere is roughly $1 - \exp[-N_p^{1/3}]$, and in optically thin fields $\langle N_\text{ray} \rangle \sim N_p$.

For typical parameters ($N_p \approx 5 \times 10^5$, $N_\gamma \approx 10^8$), wall times are $\sim 10^2$ CPU-hours. The adaptive nature concentrates computational effort where mass and opacity are high, and no grid artifacts occur. Compared to grid-based MCRT with $N_\text{cells}$ instead of $N_p$, gridless SPH-MCRT preserves or improves accuracy, while incurring a higher per-ray cost if many small smoothing volumes are intersected.

7. Extensions and Radiative Equilibrium

Radiative equilibrium is implemented via immediate re-emission (the Bjorkman–Wood scheme), depositing $\Delta T$ at each absorption site, then re-emitting the absorbed packet at the updated local temperature. In SPH terms, temperature at each particle is interpolated as: $$T(r_j) = \sum_i \frac{T_i m_i}{\rho_i} W(r_j - r_i, h_i)$$ After absorption, particle temperature increments $\Delta T_i$ can be apportioned with

$\Delta T_{ij} = A_{ij} \Delta T g(\tau_{ij}),$

where $g(0) = 1$, $g \to 0$ as $\tau \gg 1$, and the normalization is

$$\sum_i \frac{m_i}{\rho_i} A_{ij} g(\tau_{ij}) W(r_j - r_i, h_i) = 1$$

Alternatively, adopting Lucy’s path-length method, each smoothing volume accumulates the total path length $L_j$ of all photon segments passing through, with the basic constraint

$$\sum_\gamma \kappa L_j / (4\pi V_j) = \sum_i m_i \sigma T_i^4$$

Parameters controlling the physical model include $N_p$ (typically $10^5$–$10^7$), $N_\text{neigh} \sim 50$, smoothing lengths $h_j$ so $2h_j$ encloses $N_\text{neigh}$, $N_\gamma$ ($10^6$–$10^9$), and optical depths per packet sampled from the exponential distribution $e^{-\tau}$.

This approach allows for gridless, fully adaptive radiative transfer imaging with support for absorption, scattering, polarization, and thermal equilibrium, utilizing the SPH field at its native resolution (Forgan et al., 2010).