A Framework to Model Stellar Irradiated Disks with Frequency-dependent Absorption and Scattering Opacities in Athena++

Abstract: The frequency dependence of opacity is crucial for determining the thermal structure of protoplanetary disks, which in turn influences disk dynamics and planet formation. Yet many disk models adopt simplified thermodynamics, and common radiation-hydrodynamic approaches often use gray opacities, ignore scattering, and yield inaccurate results in regions with intermediate optical depth. We present a comprehensive framework that models stellar irradiation with frequency-dependent absorption and scattering across all optical depths using the Athena++ finite-volume code, extended with multigroup radiation transport and newly implemented radial rays to more accurately represent the stellar flux. To calibrate this framework, we focus exclusively on hydrostatic disk models, allowing us to isolate radiative effects and evaluate the method without additional dynamical complexity. Because dust opacity increases strongly with frequency, ultraviolet stellar irradiation heats the tenuous disk atmosphere while the optically thick midplane remains cooler. This vertical temperature gradient is captured more accurately when more frequency bands are used or when scattering is included. Our hydrostatic models achieve equilibrium temperatures that differ from Monte Carlo radiative-transfer benchmarks on average by 2--5% with 64 frequency bands and 7--11% with 3 bands. Reducing the number of bands lowers computational cost by at least an order of magnitude while increasing the maximum possible temperature deviation only from 8% to 19%. This calibration demonstrates the accuracy and efficiency of the framework and provides a solid foundation for future self-consistent studies of irradiated protoplanetary disks, including fully dynamical simulations and applications involving chemical processes and time-dependent stellar luminosity.

• The paper introduces a new multigroup radiative transfer framework for accurately modeling irradiated protoplanetary disks using frequency-dependent absorption and scattering opacities in Athena++.
• It employs an enhanced discrete ordinates solver with optimized radial ray integration to capture precise thermal and scattering effects in disk structures.
• Benchmarking against Monte Carlo methods demonstrates low error margins, paving the way for improved simulations of disk thermodynamics, chemistry, and planet formation.

A Framework for Frequency-dependent Absorption and Scattering in Athena++: Modeling Stellar Irradiated Disks

Overview and Motivation

The paper "A Framework to Model Stellar Irradiated Disks with Frequency-dependent Absorption and Scattering Opacities in Athena++" (2606.08859) introduces a robust, extensible methodology for modeling irradiated protoplanetary disks with full frequency-dependent absorption and scattering using the Athena++ magnetohydrodynamics framework. This work is motivated by the critical role of frequency-dependent radiative transfer in determining disk vertical structure, chemistry, and conditions for planet formation, and by the deficiencies of existing methods—gray opacities, moment-based approximations, and Monte Carlo techniques—that fail to consistently yield accurate solutions across the full range of optical depths and discard crucial physical processes (e.g., frequency-dependent scattering, beam crossing). The proposed framework targets high-fidelity coupling between stellar irradiation and disk hydro/thermodynamics, laying essential groundwork for dynamic and chemically rich models in future studies.

Numerical Methodology: Extensions to Multigroup Radiative Transfer in Athena++

The core of the contribution is an overhaul of the nonrelativistic discrete ordinates solver in Athena++ [Jiang 2021, 2022], extending it to support arbitrary multigroup (frequency-binned) absorption and scattering, and implementing radial rays for angular discretization optimized for centrally irradiated, nearly axisymmetric disk structures.

Frequency Discretization and Opacity Treatment

The system divides the radiation field across $N_f$ frequency groups, with frequency-dependent (input tabulated) absorption and scattering coefficients calculated for each band using Mie theory-based grain opacities from the DSHARP database [Birnstiel et al. 2018]. Band-mean Planck and Rosseland opacities are recalculated on the fly using local (or, in optically thin regimes, color-corrected) temperatures to ensure self-consistency between radiation and dust energy densities.

A Python-based pre-processing pipeline tabulates these means across a log-spaced grid in temperature for each band, circumventing floating-point issues (via asymptotic expansions) and accommodating custom band layouts. The codebase supports transparent switching (via runtime checks) between mean opacities evaluated at the local thermal or color temperature, essential for handling strong irradiation of cool optically thin media by hot stellar sources.

Radial Ray Enhancement

To precisely capture collimated stellar irradiation, especially for hydrostatic or nearly hydrostatic disks at large $r$, the angular quadrature is augmented with two radial rays (aligned to $+\hat{r}$ and $-\hat{r}$ at each cell), in addition to the standard local coordinate quadrature [Davis et al. 2012]. This preserves both the accuracy necessary for direct irradiation and isotropic reemission or scattering, with minor modifications ensuring quadrature weights and angular moment constraints remain satisfied.

Implicit Time Integration and Hydrostatic Equilibrium Calibration

The system supports fully implicit time integration, eliminating the Courant limitation associated with light-crossing times. For benchmarking and calibration, the authors consider only static density fields, evolving the radiation field to thermal equilibrium without solving for gas dynamics or radiation pressure feedback. Iterative convergence is greatly improved by reduction of the effective speed of light during radiative initialization, which does not affect hydrostatic results.

Benchmarking and Accuracy: Comparison with Monte Carlo and Hybrid Methods

The framework is stringently calibrated against high-fidelity Monte Carlo radiative transfer using RADMC-3D [Dullemond et al. 2012], as well as hybrid ray tracing methods in Athena++. Benchmarks utilize disk models with DSHARP dust and physical parameters chosen to probe a range of optical depth regimes and vertical temperature gradients relevant for ALMA-observed Class II disks.

Gray and Multigroup Regimes

Gray (frequency-independent) models recover analytic solutions with <5% average error at moderate optical depth and 10–20% maximum errors at high $\tau$ even for coarse resolution, consistent with known limitations of finite-volume approaches. Frequency-dependent models with $N_f=64$ converge to within 2–5% average error (maximum 8%) versus Monte Carlo in the disk atmosphere ($\tau_*<1$), and 7–11% in the midplane ($\tau_*\gg1$), demonstrating strong agreement across the critical parameter space. Decreasing $N_f$ to 3 reduces compute time by over an order of magnitude while inflating worst-case deviations only up to 19%, showing practical accuracy-performance tradeoffs.

Scattering and Vertical Thermal Structure

The inclusion of isotropic scattering is essential: neglecting it leads to systematic midplane and super-thermal atmosphere errors. The effect of anisotropic scattering (via $g$ parameterization), tested with Monte Carlo, induces only modest ($r$0 average, 9% max) deviations in thermal structure compared to isotropic cases, justifying the isotropic treatment for most disk regimes at current observational uncertainties.

Angular Resolution and Ray Geometry

The number of angular rays ($r$1) controls isotropy of disk emission and scattering. With $r$2 the method exhibits consistent convergence toward the Monte Carlo benchmarks in optically thick regions, with the higher resolution primarily affecting cooling efficiency and shielding near the midplane.

Computational Performance

Computation cost scales linearly with both $r$3 and $r$4. The optimized setup (3 custom frequency bands, $r$5) achieves accuracy comparable to typical ALMA modeling requirements at a compute expense appropriate for multidimensional parameter surveys, making self-consistent radiation hydrodynamics feasible in practice.

Implications, Limitations, and Future Directions

This framework quantifiably surpasses the fidelity of gray flux-limited diffusion, M1 closure, and hybrid ray tracing, especially for disks with strong irradiation vertical gradients or complex geometry (warps, flares, or planetary gaps). Unlike moment-based schemes, the discrete ordinates solver can accurately handle optically transitional layers and shadowing, and is generalizable to dynamic runs with full magnetohydrodynamics, chemical networks, and variable luminosity.

Practical implications include:

• Improved thermodynamics for planet formation models: Temperature structures from this approach correctly capture the steep UV heating of disk surfaces and the shielding of midplanes, which critically affects snowlines, dust coagulation, and migration torques.
• Chemistry and volatile evolution: Self-consistent, frequency-dependent irradiation can be coupled straightforwardly to time-dependent thermochemistry, enabling quantitatively robust predictions for the spatial distribution of molecular lines and condensation/sublimation fronts.
• Variability and episodic accretion: The framework's ability to handle time-dependent, spectral irradiation sources positions it well for studying accretion outbursts, dynamic circumbinary disks, and variable stellar phenomena.
• Application to transition and substructured disks: The softening of inner disk profiles (or opacity ramps) can model gaps, cavities, and shadowed regions, including self-consistent irradiation in disks with large-scale structure.

Limitations include the lack of inherent support for strongly anisotropic scattering in the discrete ordinates module (beyond what is tested here), possible performance bottlenecks in future full 3D, high-resolution dynamical runs, and the reliance on precomputed opacity tables that may not track complex dust evolution or compositional changes in situ. Furthermore, while the hydrostatic testbed is comprehensive, extension to full radiation-MHD (with dust-gas back-reaction, radiation pressure, and turbulence) will require additional calibration and potentially algorithmic tuning to maintain sub-percent energy conservation.

Conclusion

The presented framework constitutes a major technical advance for the simulation of irradiated protoplanetary disks, leveraging the flexibility and scalability of Athena++ with physical detail consistent with state-of-the-art Monte Carlo RHD codes. Its calibration and open-source release provide an immediate platform for high-fidelity, multidimensional studies of disk thermodynamics, evolution, and substructure formation under realistic, frequency-coupled irradiation and dust physics. As such, it will facilitate significant progress in our theoretical understanding of disk-driven planet formation and disk chemistry, and will support quantitative interpretation of ongoing high-resolution observations (ALMA, JWST) and future surveys.