MCFOST Radiative Transfer Code

• MCFOST is a 3D radiative transfer code designed for simulating dusty circumstellar disks and computing both continuum and line emissions.
• It employs Monte Carlo photon-packet transport, ray tracing, and flexible grid discretizations like cylindrical, spherical, and Voronoi meshes to resolve complex geometries.
• The code integrates advanced dust microphysics, detailed scattering, and a non-LTE module (MCFOST-art) for comprehensive disk modeling and coupling with external solvers.

Searching arXiv for relevant MCFOST papers to ground the article in published work. arxiv_search(query="MCFOST radiative transfer code Pinte", max_results=10, sort_by="relevance") arxiv_search(query="MCFOST radiative transfer code", max_results=10, sort_by="relevance") MCFOST—“Monte Carlo radiative transfer of protoplanetary disks”—is a public, three-dimensional continuum and line radiative-transfer code used to model dusty circumstellar disks and related environments. Originally developed for dusty protoplanetary disks, it has been extended to gas and atomic line transfer via the non-LTE module MCFOST-art. Across these configurations, MCFOST combines Monte-Carlo photon-packet transport, ray tracing, flexible spatial discretizations, and dust or gas microphysics to compute dust temperatures, mean radiation fields, spectral energy distributions (SEDs), monochromatic images, polarization, and non-LTE line formation in arbitrary geometries and velocity fields (Grimble et al., 2024, Tessore et al., 2021).

1. Core architecture and spatial discretization

MCFOST is implemented as a modern Fortran-90/2003 Monte-Carlo/ray-tracing radiative-transfer code. In the code description accompanying MCFOST-art, its principal components are organized into grid handling, dust and gas modules, and a radiative-transfer engine. The grid layer supports three spatial discretizations: cylindrical grids, spherical shells, and unstructured Voronoi tessellations. Physical quantities such as density, temperature, velocity, and level populations are defined cell-wise and are taken as constant within each cell (Tessore et al., 2021).

This architectural flexibility underlies the range of geometries used in published applications. Axisymmetric disk calculations commonly use cylindrical or spherical-polar grids, while coupled hydrodynamic applications can construct Voronoi meshes centered on SPH particles. In the on-the-fly PHANTOM–MCFOST coupling for binary-disk simulations, for example, the SPH particle set defines Voronoi cells, with density, composition, and temperature held constant per cell during radiative-transfer updates (Poblete et al., 16 Sep 2025).

The radiative-transfer engine is correspondingly split by physical regime. For dust continuum and molecular lines, MCFOST employs Monte-Carlo sampling. For atomic lines in MCFOST-art, it uses long-characteristic ray tracing, solving the unpolarized transfer equation simultaneously at all wavelengths along each ray (Tessore et al., 2021).

2. Continuum radiative transfer and dust microphysics

At the continuum level, MCFOST solves the monochromatic radiative-transfer equation

$$\frac{dI_\nu}{ds}=-\,\kappa_\nu \rho I_\nu + j_\nu \rho,$$

with dust emission written as

$j_\nu=\kappa_\nu B_\nu(T_d).$

In the Monte-Carlo scheme, photon packets are emitted from the stellar photosphere and, in iterative temperature calculations, from dust thermal emission. Packet propagation is controlled by an optical-depth draw, typically $\tau=-\ln\xi$, after which the packet either scatters or is absorbed and immediately re-emitted, enforcing radiative equilibrium (Grimble et al., 2024, Collaboration et al., 2020).

The dust temperature in each cell is obtained by satisfying the radiative-equilibrium condition

$$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$

This formulation appears across multiple MCFOST studies, including SED fitting for inclined inner disks, protoplanetary disks, and brown dwarf disks. Convergence criteria vary by application, but reported thresholds include relative temperature changes below a few percent, $\Delta T/T\lesssim 1\%$, or $\Delta T/T<10^{-3}$ in large model grids [(Collaboration et al., 2020); (Garufi et al., 2014); (Pinte et al., 2010)].

Dust scattering is treated with detailed grain optics. MCFOST computes extinction, absorption, scattering opacities, phase functions, and, where required, Mueller matrices from Mie theory; some studies also employ a Distribution of Hollow Spheres (DHS) prescription. Full Mueller-matrix scattering and polarized radiative transfer are used in polarimetric applications, while anisotropic phase functions are standard in continuum calculations (Grimble et al., 2024, Chen et al., 2024). This makes MCFOST suitable not only for SEDs and total-intensity images, but also for polarized-intensity maps and derived quantities such as scattering phase functions and polarized fractions.

3. Density parameterizations, grain populations, and disk structure

Published MCFOST applications commonly adopt an axisymmetric density field with a Gaussian vertical profile,

$$\rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}\,h(r)} \exp\!\Bigl[-\frac{z^2}{2h(r)^2}\Bigr],$$

together with a scale-height law

$h(r)=h_0\,(r/r_0)^\beta$

and a surface-density power law such as $\Sigma(r)\propto r^{-p}$ or $\Sigma(r)\propto r^{-\epsilon}$. These forms recur in brown dwarf disk grids, FT Tauri, PDS 453, and the EaRTH Disk Model, with modifications such as tapered outer edges or multi-zone prescriptions when a single radial component is inadequate [(Oliveira et al., 2013); (Garufi et al., 2014); (Martinien et al., 2024); (Grimble et al., 2024)].

Grain populations are usually parameterized by a power-law size distribution,

$j_\nu=\kappa_\nu B_\nu(T_d).$0

or equivalently $j_\nu=\kappa_\nu B_\nu(T_d).$1, between specified $j_\nu=\kappa_\nu B_\nu(T_d).$2 and $j_\nu=\kappa_\nu B_\nu(T_d).$3. The literature surveyed here includes astronomical silicate, amorphous Mg$j_\nu=\kappa_\nu B_\nu(T_d).$4SiO$j_\nu=\kappa_\nu B_\nu(T_d).$5, amorphous MgFeSiO$j_\nu=\kappa_\nu B_\nu(T_d).$6, silicate–carbon mixtures, porous grains, PAHs, and water-ice-bearing mixtures. Applications may also split the dust into multiple populations, as in HD 34700 A, or introduce compositionally distinct radial zones, as in PDS 453 [(Thi et al., 2013); (Chen et al., 2024); (Martinien et al., 2024)].

Vertical dust stratification is likewise supported. In FT Tauri, grains larger than a threshold $j_\nu=\kappa_\nu B_\nu(T_d).$7 follow

$j_\nu=\kappa_\nu B_\nu(T_d).$8

while the brown dwarf disk analysis of GY92 204 quotes a dust-settling exponent $j_\nu=\kappa_\nu B_\nu(T_d).$9. These prescriptions allow MCFOST models to represent settling, flaring, ring-like surface-brightness structures, and discontinuities in density or scale height when required by the data [(Garufi et al., 2014); (Oliveira et al., 2013)].

4. MCFOST-art and non-LTE atomic line transfer

MCFOST-art is the non-LTE multilevel atomic solver embedded in MCFOST. It is designed for multilevel atomic systems and models the close environment of stars in 3D. Its formalism couples the radiative-transfer equation,

$\tau=-\ln\xi$0

to the statistical-equilibrium equations for the level populations $\tau=-\ln\xi$1,

$\tau=-\ln\xi$2

The solver uses Multilevel Accelerated Lambda Iteration (MALI) following the Rybicki and Hummer formalism, with operator splitting, a diagonal approximate operator, full preconditioning including line overlap and cross-coupling terms, and Ng acceleration to reduce the number of non-LTE iterations (Tessore et al., 2021).

The module includes bound-bound and bound-free processes, TOPbase or hydrogenic Kramers cross-sections with Gaunt corrections, Thomson scattering, H free–free, H$\tau=-\ln\xi$3 bound–free and free–free opacity, and multiple broadening terms. Line profiles incorporate thermal plus microturbulent Doppler broadening, natural damping, pressure broadening via the Unsöld approximation and Lindholm theory, and linear and quadratic Stark recipes. Because MCFOST is intrinsically multidimensional, arbitrary prescribed velocity fields $\tau=-\ln\xi$4 are projected along each ray, with Doppler shifts handled explicitly. To avoid artifacts from large velocity gradients, cells can be subdivided whenever the projected gradient exceeds a user-set threshold (Tessore et al., 2021).

Benchmarking against TURBOspectrum and RH was carried out on MARCS stellar photospheres and the FAL-C solar atmosphere. In LTE continuum tests, agreement between MCFOST and RH is reported at the $\tau=-\ln\xi$5 level over most of the optical and near-IR, while TURBOspectrum differs by up to $\tau=-\ln\xi$6 around the H$\tau=-\ln\xi$7 minimum at $\tau=-\ln\xi$8 in the coolest model and can differ by $\tau=-\ln\xi$9 at the Balmer jump in the hottest model. For solar non-LTE hydrogen lines, departure-coefficient stratifications from MCFOST-art and RH agree to better than $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$0 everywhere; emergent disk-integrated profiles agree to $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$1 in the line wings, while Ly $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$2 core discrepancies of $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$3–$$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$4 fall below $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$5 after increasing the angular sampling. Across the LTE and non-LTE benchmarks, $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$6 in line wings, $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$7 in strong-line cores, and $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$8, with typical RMS differences of $$\int_0^\infty \kappa_\nu J_\nu\,d\nu = \int_0^\infty \kappa_\nu B_\nu(T_d)\,d\nu.$$9 between MCFOST-art and RH (Tessore et al., 2021).

5. Coupling to other frameworks and representative applications

A major mode of MCFOST usage is coupling to external thermochemical or dynamical solvers. In the DENT framework, MCFOST is interfaced with ProDiMo to compute a grid of 300 000 disc models for Herschel/GASPS. MCFOST supplies dust temperatures $\Delta T/T\lesssim 1\%$0, mean radiation intensities $\Delta T/T\lesssim 1\%$1, SEDs, and line source functions; ProDiMo then solves chemistry and gas thermal balance, returns level populations for key coolants, and MCFOST performs the subsequent line transfer. The implementation summary for the DENT grid reports 323 020 distinct dust and continuum MCFOST runs and 1 610 150 line-transfer calculations for 29 lines and five inclinations [(Pinte et al., 2010); (Woitke et al., 2010)].

In dust-continuum studies, MCFOST is routinely used to fit broadband SEDs and synthesize resolved images. For the optically thin disk around HD 141569A, standard public MCFOST was used to compute the 2D dust temperature structure, derive the emergent SED from UV to mm wavelengths, and synthesize the $\Delta T/T\lesssim 1\%$2 image used for comparison with VISIR; those dust results then served as input to ProDiMo for gas-line modeling (Thi et al., 2013). In the FT Tauri study, MCFOST delivered a dust model with $\Delta T/T\lesssim 1\%$3, $\Delta T/T\lesssim 1\%$4 AU, $\Delta T/T\lesssim 1\%$5 AU, $\Delta T/T\lesssim 1\%$6, and $\Delta T/T\lesssim 1\%$7, reproducing the SED from the optical edge through the far-IR PACS photometry (Garufi et al., 2014). For brown dwarf disks in $\Delta T/T\lesssim 1\%$8 Ophiuchi, grids of synthetic MCFOST SEDs were compared to Herschel and ancillary photometry through Bayesian inference, constraining inner radii, flaring indices, and scale heights (Oliveira et al., 2013).

MCFOST is also used in interferometric and spectro-interferometric modeling. In the RY Lup analysis, MCFOST-based image synthesis and Fourier transformation were used to reproduce PIONIER and GRAVITY visibilities and closure phases, yielding a final model with $\Delta T/T\lesssim 1\%$9, $\Delta T/T<10^{-3}$0 au, $\Delta T/T<10^{-3}$1, and $\Delta T/T<10^{-3}$2 mag (Collaboration et al., 2020). In non-axisymmetric magnetospheric accretion, the 2021.06 branch of MCFOST was configured with a 20-level hydrogen atom, accelerated $\Delta T/T<10^{-3}$3-iteration, and monochromatic image ray tracing to compute Br$\Delta T/T<10^{-3}$4 line profiles, synthetic channel maps, visibilities, phases, and photo-centre shifts. That study found rotationally modulated inverse P Cygni profiles, line-emitting radii ranging from $\Delta T/T<10^{-3}$5 to $\Delta T/T<10^{-3}$6 of the truncation radius, and velocity-dependent photo-centre tracks that recover the 3D tilt of the flow (Tessore et al., 2023).

Polarimetry and ice spectroscopy form another major application class. For HD 34700 A, MCFOST was used to generate total-intensity and polarized-intensity images across six near-IR wavelengths, enabling qualitative constraints on grain-size mixtures, porosity, and Mie versus DHS grain optics (Chen et al., 2024). For PDS 453, MCFOST models with a two-zone radial structure, mixed silicate-and-ice grains, and $\Delta T/T<10^{-3}$7 reproduced the two-nebula morphology, the ring-like brightening at $\Delta T/T<10^{-3}$8 au, polarization levels, and a $\Delta T/T<10^{-3}$9 water-ice band depth of $$\rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}\,h(r)} \exp\!\Bigl[-\frac{z^2}{2h(r)^2}\Bigr],$$0 (Martinien et al., 2024). In the spatially resolved ice-band study, MCFOST was used to extract spectra from different disk locations and inclinations, showing that apparent H$$\rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}\,h(r)} \exp\!\Bigl[-\frac{z^2}{2h(r)^2}\Bigr],$$1O-ice band minima depend on the balance between absorption and scattering (Martinien et al., 20 Mar 2025).

6. Validation, limitations, and interpretive issues

Across the literature, MCFOST’s strengths are consistently associated with self-consistent dust radiative transfer, flexible geometries, detailed scattering physics, and compatibility with line, chemistry, and hydrodynamic workflows. At the same time, the same literature identifies persistent modeling limits. For MCFOST-art, current limitations explicitly include the coherent-scattering approximation for electrons in extreme winds, approximate Stark broadening limited to the first series, and memory overhead for large atoms in 3D when Ng acceleration is used (Tessore et al., 2021).

In continuum and polarimetric applications, limitations often arise from the adopted dust optics and geometry rather than from the transport solver itself. The HD 34700 A study reports that no model matched all observed properties of the disk, and attributes this in part to the limits of Mie-theory calculations for irregular aggregate grains and to the simplified axisymmetric disk setup (Chen et al., 2024). The PDS 453 analysis similarly required a two-zone geometry because no one-zone model could reproduce the observed ring (Martinien et al., 2024). In the $$\rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}\,h(r)} \exp\!\Bigl[-\frac{z^2}{2h(r)^2}\Bigr],$$2 Eridani debris-disk simulations, the authors report that MCFOST was used in a largely black-box fashion, with ring geometry, dust mass, astrosilicate composition, and $$\rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}\,h(r)} \exp\!\Bigl[-\frac{z^2}{2h(r)^2}\Bigr],$$3 specified, but without reporting power-law slopes, scale-height functions, photon-packet counts, or explicit scattering settings (Bao et al., 20 Sep 2025).

Interpretation of observables can also remain degenerate even when the radiative-transfer calculation is well converged. Within the DENT grid, the analysis of far-IR lines is described as complex, and the inversion of line fluxes into physical quantities as difficult; the recommended strategy is to combine Herschel fine-structure lines with continuum data and rotational CO lines in the (sub-)millimetre regime (Pinte et al., 2010). For ice spectroscopy, the band minimum of H$$\rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}\,h(r)} \exp\!\Bigl[-\frac{z^2}{2h(r)^2}\Bigr],$$4O at the central source can shift by up to $$\rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}\,h(r)} \exp\!\Bigl[-\frac{z^2}{2h(r)^2}\Bigr],$$5, comparable to the difference expected between amorphous and crystalline ices, showing that apparent spectral shifts need not map one-to-one onto composition without tailored radiative-transfer modeling (Martinien et al., 20 Mar 2025).

These published results place MCFOST in a specific methodological niche. It is not only a dusty SED code, nor only a Monte-Carlo image generator, but a radiative-transfer framework that has been extended from passive dust-heating problems to non-LTE atomic line formation, interferometric forward modeling, thermochemical pipelines, and on-the-fly coupling to hydrodynamics. The applications explicitly listed for MCFOST-art include stellar photospheres, chromospheres, winds, protoplanetary-disc plus magnetosphere interaction zones, molecular-line disks, polarized line transfer as future work, and planet-forming regions (Tessore et al., 2021).