Simulating the thermodynamics of gas, radiation, and multispecies dust: Method and applications to protoplanetary disks

The thermodynamics of protoplanetary disks, governed by a complex interplay between gas, dust, and radiation, exerts a strong influence on their morphology and dynamics at both large and small scales. Historically, hydrodynamical simulations have treated thermodynamics with approximate prescriptions such as local isothermality or parametrized β-cooling. The authors devised a radiation-hydrodynamics scheme for the PLUTO code that includes energy exchange (via absorption, emission, and collisions) for gas, radiation, and multiple species of dust. Dust-gas dynamics are handled by treating each dust species as a pressureless, diffusive fluid to represent changes in disk illumination caused by grain settling and trapping.

The paper addresses the limitations of existing modeling approaches. Thermochemical models like ProDiMo or DALI are too computationally expensive to couple with gas dynamics because energy exchange must be recomputed at each timestep. Simulations have historically assumed locally isothermal equations of state or adiabatic equations with exponential thermal relaxation to a fixed background. In contrast, the authors advocate for physically-motivated thermodynamic prescriptions based on detailed calculations of equilibrium dust distribution, temperature structure, and collisional/radiative relaxation timescales. They introduce a three-temperature (3T) scheme based on the M1 method for the PLUTO code that self-consistently evolves the separate energies of gas, dust, and radiation via thermal accommodation, absorption/emission of radiation, and stellar irradiation to allow dust and gas temperatures to decouple in upper layers.

The authors compare modeling dust as superparticles (which suffers from Poisson noise) versus a pressureless-fluid approach (which loses realism when Stokes number St ≡ tsomegaK ≳ 1). They also note that using multiple pressureless fluids with collisional β-cooling rates, or local dust abundance to compute opacities without considering collisional thermal relaxation, are existing approaches.

The paper develops a numerical method to compute the dynamical and thermodynamical coupling between gas, radiation, and multiple species of dust within the grid-based PLUTO code. The solution strategy involves an implicit-explicit (IMEX) scheme (IMEX1) to combine implicitly integrated source terms with non-stiff explicit transport terms. The simulation timestep is determined by the CFL criterion. A Strang-split approach is used for radiation blocks, where energy and momentum exchange between gas, dust, and radiation are computed implicitly, while radiation transport is computed explicitly within substeps of length ∆trad to ensure numerical stability.

The energy and momentum updates are performed using a Newton-Raphson technique. The linearized energy update (Equation A.1a) shows the coupling between the radiation field, gas internal energy, and dust internal energy:

i+1 δXg /∆t = − δEgi+1 /∆t + cκi+1 − 4ar Tg3,i+1 Wgi+1 δEgi+1 g ρg (δEr −4ar Tgj=1

The momentum update involves solving a system of equations where the dust species do not directly interact with one another but influence each other through their effects on gas and dust energy and momenta. The resulting reduced Jacobian matrices allow for a significant reduction in computational complexity from O((nd + 2)3) to O(nd).

Conservation laws are verified, showing that the implicit step conserves the sum of all internal energies and momenta, modified only by the injection of energy to the system via stellar irradiation. The conservation law for total energy (Equation A.11) shows that the residuals cancel, and terms of the form Exi+1,∗k + δExi+1,∗k have been identified with the energy value at the next iteration.

The paper presents several test problems to verify functionality:

Energy and momentum coupling:

"In Figure 2, we plot the evolution of the temperatures and velocities of all species in the top and middle panels. Over the lifetime of the simulation, all species relax to the same temperature and velocity, at which point the energy and momentum exchange terms net to zero."

Ray-traced irradiation:

The authors test ray-traced starlight by assigning frequency-dependent absorption opacities (Equation 33) and initial fluxes based on a sum over frequency bands (Equation 34). They show that the resulting gas-pressure gradient drives motions, and the radiation flux couples to dust momentum.

Ray-traced force on dust:

They verify the radiative force by setting all stopping times to a large value and decoupling dynamics, showing that the velocity of each grain species would be expected to evolve as vd,j (t) = Fr σd,j t, with relative deviation less than 10−3 at all times.

Multispecies dusty shock:

The method is tested in the presence of shocks, where analytical solutions for dust velocities are compared to numerical simulations, finding close agreement between simulation outputs (open circles, sampled every 5 grid cells) and analytical theory (solid curves).

Diffusion in a homogeneous medium:

For pure diffusion, the analytical evolution of dust density is solved using Equation (39), which yields the classical Gaussian vertical profile of protoplanetary disks. The numerical solution agrees well with this analytical form across many orders of magnitude in density.

Vertical dust settling:

The vertical structure is analyzed by balancing gravity and pressure support, leading to a characteristic dust scale height derived from Equations (40) and (42). They find that the dust scale height is Hd,j = p / (1 + Stj /α).

Radial drift:

The radial drift velocity is governed by Equation (45), which shows that for large bodies (St ≫ 1), angular momentum loss and thus radial drift are very slow. They verify the correct implementation of geometric source terms in curvilinear coordinates, finding that the numerical results agree with analytical predictions, with the error decreasing linearly with increasing Stokes number St.

Radial dust trapping:

They test a scenario where pressure bumps create a local maximum in radial drift velocity (Equation D.11). The equilibrium between dust-gas drag and diffusion sets the ring width, yielding a characteristic dust concentration width of wdg ≡ −α(St + St−1)/(p−1g ∂R pg)Rmax, which shows that both small and large grains would not concentrate at the pressure bump.

Application: Dust drift and sublimation:

In the final application, they simulate dust drifting inward toward the star. They incorporate a novel numerical method where dust opacity and stopping time are functions of temperature (Equation 52), allowing for a phase transition from solid to gas at sublimation temperatures (e.g., Tsub ≈ 150 K for water ice). This treatment demonstrates that the simulation accurately reproduces the true steadystate solution even when local thermodynamic equilibrium (LTE) between gas, dust, and radiation cannot be assumed.

The authors conclude by summarizing their method as a novel numerical method that solves for thermodynamic coupling between gas, dust, and radiation. They state that this enables the computation of variables like temperatures, velocities, and densities self-consistently within a single simulation. This capability will significantly benefit investigations of how large-scale disk substructure is driven by planets or shadows. The methods presented are a starting point for future improvements, such as incorporating photoelectric heating or gas line transitions. Additionally, they suggest generalizing frictional heating to apply to each dust species separately and directly incorporating dynamical coupling into the solution of the Riemann problem.

The paper concludes with an extensive list of references covering topics from hydrodynamics and radiative transfer to specific dust dynamics applications. The simulations were carried out on VERA, FREYA, and ORION clusters at MPCDF in Germany. The authors acknowledge helpful discussions on scientific applications with various researchers and thank the anonymous referee for suggestions that improved the quality of the paper.

The final conclusion emphasizes that by incorporating multiple species of dust, their individual self-consistent dynamics (drag and diffusion), this method builds upon its precursor, which only used an explicit step in the temperature scheme. They can therefore define dust opacity as κd,j = Sj (Td,j)κ0 and consider a single dust species perfectly coupled to the gas. These changes make it possible to reproduce spa[cial] variations in grain opacity and size distributions by incorporating phenomena such as radial drift and dust settling. This reproduces more [realistic] observational features of protoplanetary disks. Monte Carlo radiative transfer (MCRT) tools can then be used to convert hydrodynamical simulations into multiwavelength mock images for comparison with observations, where the new method allows these variables to be computed self-consistently within a single simulation, greatly improving the predictive power of mock observations. This would significantly benefit investigations of how large-scale disk substructure, whether driven by planets (Speedie et al. 2022), shadows (Zhang & Zhu 2024), or infalling streamers (Kuznetsova et al. 2022), and their observational appearance in near-infrared scattered light (tracing small grains), millimeter/submillimeter continuum (tracing large grains), and various molecular lines (tracing the gas). Circumplanetary disks are also a potential application, where thermal relaxation timescales depend on the dynamically evolving size distribution of dust in the circumplanetary region. The opacities and stopping times which mediate thermal relaxation are, in turn, determined by the dynamically evolving size distribution and abundances of dust in the circumplanetary region."

The paper is structured with sections detailing equations (2) through (53), solution strategies (3), and test problems (6). Figures 1 through 13 provide graphical summaries of the simulation results. The Appendix details the linearized energy and momentum updates, including the full set of Jacobian matrices for each species. Appendix D provides derivations for key source terms like irradiation heating and ring-trapping tests. The overall structure is designed to demonstrate a comprehensive numerical method capable of handling the stiffness inherent in coupling gas, radiation, and multiple dust species across various physical regimes.

The authors explicitly state that the new method allows them to reproduce the twotemperature limit used in most radiation-hydrodynamics schemes when dust and gas are tightly coupled. They also confirm that they have designed a novel numerical method where Tsub,j = 2000 K + j × (500 K) is a smooth planetary disks.

The final summary of the paper is:

"We have designed and implemented a novel numerical method that solves for thermodynamic coupling between gas, dust, and radiation, and in a number of test problems, where Tsub,j = 2000 K + j × (500 K) is the sublimation temperature. By incorporating multiple species of dust, our method builds upon its precursor [which only used an explicit step] to ensure numerical stability in the iterative each with its own self-consistent dynamics (drag and diffusion), and δsub = 10−10 is a small factor to pre-fusion), our method builds upon its precursor, which only used an explicit step. We can therefore define κd,j = Sj (Td,j)κ0 and consider a single dust species perfectly coupled to the gas. These changes make it possible to reproduce spa[cial] variations in grain opacity and size distributions by incorporating phenomena such as radial drift (e.g., Dullemond et al. 2018) and dust settling (e.g., Bae et al. 2021). These comprehensive treatments (e.g., Wang et al. 2023) in the set the local frequency-integrated opacities and gas-grain thermal accommodation timescales, and in turn the heating, cooling, and background temperature structure of the disk as a whole."

The paper is essentially a detailed description of a new radiation-hydrodynamics scheme for simulating protoplanetary disks that self-consistently couples gas dynamics with multiple species of dust by treating dust as pressureless fluids. It details the numerical implementation using an IMEX1 scheme and Newton-Raphson iterations to solve the coupled energy and momentum equations, ensuring conservation laws are satisfied in each iteration. The paper rigorously tests this method across various physical scenarios, including thermal relaxation times, radiative heating from starlight (ray-traced irradiation), collisional coupling, shocks, diffusion in homogeneous media (leading to Gaussian profiles), vertical settling of dust grains to the midplane, and radial drift toward pressure maxima or sublimation fronts. The key innovation lies in treating the dust energy and temperature evolution self-consistently with radiation and gas via coupled source terms derived from detailed physical processes like collisional coupling (Equation 8) and frictional heating (Equation 9), while also incorporating a gradient-diffusion prescription for dust transport (Equation 11). This allows for the computation of variables like temperatures, velocities, and densities self-consistently within a single simulation. The method is shown to be numerically stable and accurate across various resolutions, with errors decreasing as the timestep is refined. It successfully reproduces the twotemperature limit in tightly coupled regimes and accurately models complex phenomena such as dust ring formation at pressure bumps and sublimation fronts, demonstrating its utility for understanding planet formation processes in protoplanetary disks."

The paper is essentially a detailed description of a new radiation-hydrodynamics scheme for simulating protoplanetary disks that self-consistently couples gas dynamics with multiple species of dust by treating dust as pressureless fluids. It details the numerical implementation using an IM

Based on this scientific paper, here are specific improvements that could be made to AI systems, categorized by the capabilities they would gain:

---

)AI System Improvement 1: Multiphysics Simulation Engine for Protoplanetary Disks (The Core Method)

This paper introduces a novel radiation-hydrodynamics scheme coupled with multiple pressureless dust fluids, treating gas, radiation, and multispecies dust self-consistently.

Improvements to the AI system would involve implementing this scheme as a core module within a general-purpose astrophysical simulation framework.

1. The AI system should be capable of solving complex conservation laws (Equations 1a-1d) incorporating:

2. A multi-fluid approach (gas + multiple pressureless dust species, treated as tracers/diffusive fluids).

3. Self-consistent energy exchange mechanisms including absorption, emission, and collisions between gas, radiation, and all dust species (Equations 7 through 9).

4. Dynamic radiation transport using a frequency-integrated moment method with the M1 closure (Equation 3), allowing for the evolution of background temperature profiles.

5. Explicit treatment of collisional coupling terms (momentum exchange and frictional heating) that depend on local dust abundance and grain size distribution (Equations 7-9).

6. Modeling dust advection and diffusion using a gradient-diffusion prescription, including the calculation of momentum transport due to this flux (Equations 11-12).

The improved AI system could:

• Accurately predict the morphology and dynamics of protoplanetary disks (e.g., spiral structures, gaps) under conditions where dust settling and trapping significantly alter radiation profiles.

• Simulate the thermal decoupling between gas and dust in different disk regions (atmosphere vs. midplane) based on grain size distribution and stopping times.

• Model the effects of planetary accretion or shadowing on radial temperature gradients by treating radiation transport dynamically rather than using fixed background profiles.

)AI System Improvement 2: Self-Consistent Thermodynamic Modeling for Chemical/Physical Properties

The paper demonstrates how dust properties (opacity, specific heat capacity, sublimation temperature) evolve based on their local thermodynamic state (temperature).

Improvements to the AI system would involve integrating this feedback loop into property calculation routines.

1. The AI system should incorporate a dynamic relationship where dust properties are functions of the local gas/dust temperature:

2. Dust opacity and specific heat capacity should be calculated using temperature-dependent functional forms (Equations 5, 33).

3. The sublimation front for various dust species must be dynamically determined based on their evolving temperature profile (Equation 52).

The improved AI system could:

• Simulate the chemical evolution of disk material by tracking which dust species are sublimated and entrained in the gas due to local heating events (e.g., near a star or planet).

• Predict the resulting changes in dust size distribution and abundance profiles as a function of disk evolution, crucial for understanding planetesimal formation pathways.

)AI System Improvement 3: Advanced Numerical Solver for Stiff Systems

The paper utilizes an Implicit-Explicit (IMEX) scheme combined with Newton-Raphson iteration and Jacobian matrix reduction to solve the coupled, potentially stiff equations efficiently.

Improvements to the AI system would involve optimizing its numerical backbone for long-term stability and speed.

1. The solver should employ an IMEX time integration scheme (Equation 16) where radiation/source terms are treated implicitly and hydrodynamical advection/diffusion is treated explicitly.

2. It must utilize algebraic rowreduction of the Jacobian matrices for energy and momentum exchange to reduce computational complexity from cubic to quadratic in the number of unknowns (Equations 4.3, 4.4).

3. The system should be able to handle stiff dynamics where collisional coupling times are much shorter than hydrodynamic timescales (i.e., when dust and gas temperatures rapidly equilibrate) using appropriate implicit treatments for the source terms (Appendix A).

The improved AI system could:

• Perform long-term, stable simulations of protoplanetary disks over many dynamical times without accumulating significant numerical errors in total energy or momentum conservation.

• Achieve high computational efficiency by solving the massive coupled system through reduced Jacobian methods, making complex multi-species simulations tractable for large grids.

)AI System Improvement 4: Predictive Modeling for Astrophysical Observables (Synthetic Data Generation)

The paper provides a framework to compute self-consistent thermodynamic variables (T, v, density) without relying on pre-assumptions like local thermal equilibrium (LTE).

Improvements to the AI system would involve using the simulation output as a high-fidelity data generator.

1. The AI system should be capable of running simulations that yield full spatio-temporal fields for gas temperature, dust temperature, and radiation energy density (Figures 2, 3).

2. It must integrate these fields with Monte Carlo radiative transfer (MCRT) tools to generate synthetic observational data across various wavelengths (e.g., near-infrared scattered light for small grains, millimeter continuum for large grains).

The improved AI system could:

• Generate realistic mock observational data sets that accurately reflect the complex interplay between gas, dust, and radiation in protoplanetary disks.

• Test and validate other AI models or physical theories by comparing their predictions against these high-fidelity synthetic datasets, thereby accelerating discovery in astrophysical modeling.