Monte Carlo Radiative Transfer in SPH

• Monte Carlo Radiative Transfer (MCRT) is a stochastic technique that simulates photon packet interactions with matter in complex astrophysical media.
• It leverages SPH’s gridless methodology to evaluate densities, emissivities, and optical depths using kernel-based interpolation, preserving adaptive resolution.
• The method minimizes geometric artifacts and efficiently models synthetic images and spectra, making it vital for studying protoplanetary disks, star formation, and turbulent flows.

Monte Carlo Radiative Transfer (MCRT) is a stochastic numerical technique that simulates the interaction of radiation with matter by following the trajectories of discretized radiation packets (or photon packets) through a model medium composed of gas, dust, or plasma. The approach is widely employed in astrophysics due to its flexibility in handling complex geometries, multi-dimensional and time-dependent problems, and arbitrary microphysical processes such as absorption, scattering, and re-emission. MCRT has been implemented "natively" on meshless data structures, most notably Smoothed Particle Hydrodynamics (SPH) simulations, bypassing the need for grid-based discretization and leveraging the full adaptive resolution of particle hydrodynamics (Forgan et al., 2010). This gridless approach enables the direct construction of synthetic images and spectra with fidelity matching that of the underlying hydrodynamics.

1. Direct MCRT in SPH Density Fields

Conventional MCRT tracks photon packets as they propagate, scatter, and are absorbed in a density field—traditionally represented on a grid. The innovation in (Forgan et al., 2010) is to perform all radiative operations directly on an ungridded SPH dataset, evaluating densities, opacities, and emissivities with the same kernel-based methods as the hydrodynamics:

• Point Source Emissivity: For stellar sources, the emitted luminosity follows a blackbody spectrum,

$$L_{\star} = 4\pi^2 R_{\star}^2 \int_{\nu_{\min}}^{\nu_{\max}} B_{\nu}(T_s)d\nu$$

• Diffuse Gas Emissivity: For optically thin thermal gas, the SPH emission per-particle is

$$L_{\text{gas}} = \sum_{i} 4\pi m_i \int_{\nu_{\min}}^{\nu_{\max}} \kappa_{\nu} B_{\nu}(T_i)d\nu$$

• Kernel-Based Emissivity Interpolation: The local emissivity at position $\mathbf{r}$ is sampled as

$$\varepsilon_{\nu}(\mathbf{r}) = 4\pi \sum_j \frac{\kappa_{\nu} B_{\nu}(T_j) m_j}{\rho_j} W(\mathbf{r} - \mathbf{r}_j, h)$$

where $W$ is the SPH smoothing kernel and $h$ the smoothing length.

Photon paths are generated such that at every step the cumulative optical depth is evaluated using the native SPH density field: $$\tau_{\text{scatter}} = \int_0^L \rho \chi_{\nu} d\ell$$ with the cross-section $\chi_{\nu}$ combining absorption and scattering. Scattering angles are sampled from distributions set by the local microphysics, generally by inverting or sampling the cumulative distribution functions of phase functions or Mueller matrices (for polarization).

Imaging is performed by binning the positions of photons as they exit the domain onto a detector plane, with the resulting statistics (e.g., photon counts per pixel) interpreted as synthetic images or datacubes.

2. Optical Depth and Ray Tracing on Meshless Fields

A critical computational detail is the evaluation of column densities and optical depths along photon trajectories in the meshless SPH formalism. Scalar fields $A(\mathbf{r})$ are interpolated in SPH as: $$A(\mathbf{r}) = \sum_{j} \frac{A_j m_j}{\rho_j} W(\mathbf{r} - \mathbf{r}_j, h)$$ and the density by: $$\rho(\mathbf{r}) = \sum_{j} m_j W(\mathbf{r} - \mathbf{r}_j, h)$$

For a photon path, the integrated column density can be recast as: $$\Sigma = \int_0^L \rho(\mathbf{r}) d\ell = \sum_{j} \int_0^L m_j W(\mathbf{r} - \mathbf{r}_j, h_j) d\ell$$ This reduces to calculating the overlap of the ray with each particle's smoothing kernel, typically (for spline kernels) an analytic or semi-analytic integral dependent on particle–ray impact parameter. Only particles whose smoothing sphere ($\sim2h$ radius) intersects the ray contribute; candidate sets are efficiently identified using hierarchical spatial partitioning (e.g., octrees with Axis-Aligned Bounding Boxes).

The location of the next interaction is found by accumulating optical depth $\tau$ along the path and solving for the point where it matches a randomly sampled value $-\ln(1 - \zeta)$, using, e.g., recursive bisection.

3. Preservation of Adaptive Resolution and Gridless Advantages

The gridless SPH-MCRT approach offers the following fundamental advantages over grid-based radiative methods:

• Preservation of native SPH resolution: No information is lost to grid interpolation; radiative calculations retain all spatial detail present in the original hydrodynamics, including in regions of rapid density variation.
• Adaptive sampling: High-density (thus, more highly resolved) regions inherently receive more detailed radiative transfer treatment, reflecting the greater number of SPH particles present.
• Artifact avoidance: By not imposing grid geometry, the method minimizes geometric artifacts that could bias radiative transport or synthetic image features, which is critical for high-resolution observational comparison.
• Efficient intersection tests: Octree spatial acceleration avoids the high computational cost that would otherwise arise from checking all particles for every photon.

In sum, these features make the method particularly well suited for synthetic imaging and spectral modeling in scenarios with highly variable or complex spatial structure, such as protoplanetary disks, star formation simulations, and the fragmentation of massive structures.

4. Astrophysical Applications

The meshless SPH-MCRT formalism has been explicitly demonstrated in several astrophysical contexts:

Model/Application	Description and Demonstrated Results
Protoplanetary disc (e.g. HL Tau)	Imaging of an SPH model reproducing candidate protoplanet and spiral arms; brightness and detectability functions evaluated for ALMA.
Stellar disc encounter/outbursts	Imaging of protostellar discs under strong tidal encounters; disc distortion and tidal arm production imaged with spatial and temporal fidelity.
Star and planet formation scenarios	Improved realism in synthetic observation compared to gridded approaches; near-1:1 mapping with underlying SPH simulation resolution.

In the HL Tau case, the gridless SPH-MCRT implementation could robustly recover the expected features of a forming planet and its surrounding spiral structure. Matches to expected instrument wavelength ranges and spatial resolution were achieved by post-processing the ungridded SPH data (Forgan et al., 2010).

5. Future Extensions: Radiative Equilibrium and Beyond

While the presented method is applied to post-processing with fixed, precomputed temperature distributions from SPH, the framework is structurally extensible to full radiative equilibrium iterative schemes within the meshless context. Energy absorbed by photon packets can be distributed back to the underlying SPH particles directly using smoothing kernel interpolation: $$\Delta T(r_j) = \sum_i \frac{\Delta T_{ij} m_i}{\rho_i} W(r_i - r_j, h_i)$$ with closure achieved by normalizing temperature increments and selecting flux-decreasing spatial distributions, e.g. $g(r_{ij}) = e^{-\tau_{ij}/4}$.

Such an approach would enable direct coupling of radiative transfer to thermal evolution and, ultimately, to the hydrodynamics, without ever interpolating to a grid. This extension opens the path to self-consistent, on-the-fly, gridless radiation hydrodynamics with full adaptive fidelity, and to future extensions including non-LTE line radiative transfer.

6. Conclusions and Theoretical Significance

The SPH-native MCRT method as developed provides a rigorous, high-fidelity route to synthetic imaging and spectral modeling that avoids the main limitations of traditional grid-based mesh radiative transfer in dynamically adaptive, particle-based simulations. It maintains spatial accuracy, leverages efficient hierarchical methods for ray-particle intersection, faithfully preserves substructure, and is directly extensible to iterative radiative equilibrium. As such, the method is a valuable tool for the interpretation of high-resolution observations of star formation regions, discs, and turbulent hydrodynamic flows captured in modern SPH simulations (Forgan et al., 2010).