  We propose in this work a Monte Carlo method for three dimensional scalar
radiative transfer equations with non-integrable, space-dependent scattering
kernels. Such kernels typically account for long-range statistical features,
and arise for instance in the context of wave propagation in turbulent
atmosphere, geophysics, and medical imaging in the peaked-forward regime. In
contrast to the classical case where the scattering cross section is
integrable, which results in a non-zero mean free time, the latter here
vanishes. This creates numerical difficulties as standard Monte Carlo methods
based on a naive regularization exhibit large jump intensities and an increased
computational cost. We propose a method inspired by the finance literature
based on a small jumps - large jumps decomposition, allowing us to treat the
small jumps efficiently and reduce the computational burden. We demonstrate the
performance of the approach with numerical simulations and provide a complete
error analysis. The multifractional terminology refers to the fact that the
high frequency contribution of the scattering operator is a fractional
Laplace-Beltrami operator on the unit sphere with space-dependent index.
