Overview

Scientific context

Carbon-rich asymptotic giant branch (AGB) stars lose mass through a slow, massive, dust-driven wind. In the cool, extended atmosphere, carbon that is not locked into carbon monoxide condenses into small amorphous-carbon grains. These grains are extremely opaque to the stellar radiation field; the momentum they absorb is shared with the gas through collisions, and the resulting radiative acceleration, once it exceeds gravity, drives a stationary outflow. The wind is therefore an intrinsically coupled problem: the temperature structure sets the chemistry and the dust nucleation; the dust sets the opacity; the opacity sets both the temperature (through radiative transfer) and the radiative acceleration (through the wind dynamics); and the wind density feeds back on the chemistry and dust.

This code solves the stationary, spherically symmetric version of that problem. All time derivatives are dropped; the model seeks the steady-state structure \(\{T_\mathrm{gas}(r),\,T_\mathrm{dust}(r),\,\rho(r),\,v(r)\}\) together with the mass-loss rate \(\dot{M}\) that is consistent with the radiation field, the chemistry, and the dust.

The scheme follows the coupled-model approach of Winters et al. and Gail & Sedlmayr, with the singularity-free wind formulation of Melia.

What the code computes

Given a small number of stellar parameters (radius, mass, luminosity, C/O ratio, and either a prescribed or eigenvalue mass-loss rate) and a starting structure, the model produces:

  • the radial temperature structure of gas and dust, in radiative equilibrium with the frequency-dependent radiation field;

  • the gas-phase chemical composition at every radius (FastChem);

  • the dust distribution — nucleation rate, grain-growth timescale, grain number density, mean grain radius, and the degree of carbon condensation — from the Gail & Sedlmayr moment method;

  • the wind structure — velocity \(v(r)\), density \(\rho(r)\), the location of the sonic (critical) point, and the mass-loss rate \(\dot{M}\);

  • the emergent spectrum and the radial run of the radiation moments.

The solution strategy in one paragraph

The model is solved by fixed-point (Picard) iteration over two nested loops, as sketched in Winters’ coupled-model diagram:

  1. Outer loop (global iteration): update the chemistry, dust and wind structure at the current temperature, then perform one radiative-equilibrium temperature correction. Repeat until the temperature has converged.

  2. Inner loop (chemistry–hydrodynamics iteration): at fixed temperature, alternate the chemistry+dust calculation, the radiative transfer, and the wind solve until the radiative acceleration parameter \(\alpha\) is self-consistent.

Each physics block is a separate module with a narrow interface; see Architecture. The numerical robustness machinery (adaptive under-relaxation, Anderson acceleration, the two-phase Unsöld–Lucy → linearisation corrector, the optional movable grid) is layered on top of this loop and is described in Numerical methods and robustness.

Status and validation

  • The radiative-transfer module reproduces the thesis scheme term by term; the Taylor (finite-difference) moment discretisation is unconditionally stable, whereas the cubic-spline discretisation loses positivity for optically thick cells and is provided mainly for cross-checks.

  • With self-consistent carbon depletion in the dust moment sweep, the eigenvalue mass-loss rate for IRC+10216 settles at \(\dot{M}\approx 7.8\times10^{-5}\,M_\odot\,\mathrm{yr}^{-1}\), a ~3 % match to the observed value.

  • The default wind solver is the Melia Φ shooting method with an eigenvalue \(\dot{M}\); the Henyey global-Newton structure solver is implemented and available but is not required for a converged model (see Wind hydrodynamics).