Dust formation (Gail & Sedlmayr)

Module: src/dust/ — agb::GailSedlmayrDust (derived from the abstract agb::DustSpecies), with the thermochemistry in gail_sedlmayr_thermo_functions.cpp. An agb::AnalyticDust prescription is also available for testing.

Carbon condenses out of the cooling, expanding gas into small amorphous-carbon grains. The code follows the classical Gail & Sedlmayr moment method: homogeneous nucleation of seed clusters followed by grain growth, tracked through the moments of the grain-size distribution.

The moment equations

Let \(K_j(r)\) be the \(j\)-th moment of the grain-size distribution (number of monomers per grain to the power \(j/3\)), normalised per hydrogen nucleus. In the stationary, spherically symmetric wind they obey the coupled system

(1)\[\begin{split}v\,\frac{\mathrm{d}K_0}{\mathrm{d}r} &= \frac{J_\star}{n_{\langle H\rangle}}, \\ v\,\frac{\mathrm{d}K_j}{\mathrm{d}r} &= \frac{j}{3}\,\frac{1}{\tau}\,K_{j-1} + N_\ell^{\,j/3}\,\frac{J_\star}{n_{\langle H\rangle}}, \qquad j = 1,\dots,5 ,\end{split}\]

where

  • \(J_\star\) is the stationary nucleation rate (seed clusters per unit volume and time),

  • \(1/\tau\) is the net monomer growth rate (so \(\tau\) is the grain-growth timescale),

  • \(N_\ell\) (minimum_monomer_number = 1000) is the lower cluster-size bound of the integration, and

  • \(n_{\langle H\rangle}\) is the hydrogen number density.

The code integrates six moments (\(K_0\dots K_5\)). Physically: \(K_0\) is the grain number density, and \(K_3\) is the condensed volume per hydrogen nucleus — i.e. the carbon locked into dust. The mean grain radius follows from \(\langle a\rangle = a_0\,(K_3/K_0)^{1/3}\) with the monomer radius \(a_0=1.28\times10^{-8}\,\mathrm{cm}\).

Nucleation rate

The nucleation rate is built from classical nucleation theory (nucleationRate()):

\[J_\star = c_0\,\beta\,A_\star\,Z ,\]

assembled from

  • the supersaturation ratio \(S(T,n_C)\) and its logarithm (saturationRatio, saturationVapourPressure); nucleation is shut off (\(J_\star\to0\)) where \(\ln S\le0\);

  • the critical cluster size \(n_\star\) (criticalClusterSize);

  • the free energy of cluster formation \(\Delta F\) (freeEnergyOfFormation), using the surface tension and the temperature scale \(\theta_\infty=\sigma A_1/k_B\);

  • the monomer growth rate onto the critical cluster \(\beta\) (monomerGrowthRate);

  • the Zeldovich factor \(Z\) (zeldovichFactor);

  • the equilibrium cluster distribution \(c_0\) (equilibriumClusterDistribution); and

  • the critical-cluster surface area \(A_\star = A_1 n_\star^{2/3}\).

Growth is fed by collisions of the grains with the carbon-bearing species \(\mathrm{C}, \mathrm{C_2}, \mathrm{C_2H}, \mathrm{C_2H_2}\), weighted by their sticking coefficients and thermal velocities (growthRate()). The material constants (graphite/amorphous-carbon: surface tension \(1400\,\mathrm{erg\,cm^{-2}}\), critical saturation ratio 3, sticking coefficients 0.37/0.34) are the Gail & Sedlmayr (1984) values.

Carbon depletion in a single causal sweep

The decisive physics: the moment equations Eq. 1 are integrated outward in one coupled RK4 sweep (agb::GailSedlmayrDust::calcDistribution()), and at each radius the growth species that feed both \(1/\tau\) and \(J_\star\) are throttled by the running degree of condensation

\[f_c = \frac{K_3}{(\varepsilon_C-\varepsilon_O)\,n_{\langle H\rangle}/n_{\langle H\rangle}} \;\equiv\; \frac{K_3}{c_\mathrm{cond}} , \qquad \text{throttle factor } = 1 - f_c .\]

Both the (linear) growth rate and the (strongly nonlinear) nucleation rate are evaluated from the depleted densities \(n\cdot(1-f_c)\), so condensation self-limits as carbon is consumed. Because the sweep marches outward and \(K_3\) only grows, \(f_c\) rises monotonically, the throttle decreases smoothly, and there is no oscillation and no outer fixed-point iteration.

Important

This differential, in-sweep depletion (Gail & Sedlmayr book §14.3; Winters thesis Eq. 5.4) was the missing physics that earlier caused the dust opacity — and hence the radiative acceleration \(\alpha\) — to grow without bound, collapsing the wind into a spurious “converged” state. With it in place the self-consistent eigenvalue mass-loss rate for IRC+10216 settles at \(\dot{M}\approx7.8\times10^{-5}\,M_\odot\,\mathrm{yr}^{-1}\) (a ~3 % match to the observed value), with a degree of condensation \(f_c\approx0.13\).

The integration detail: gas-phase densities, the growth timescale and the hydrogen density are interpolated to the RK4 midpoints (logarithmically for positive-definite quantities, linearly for temperature/velocity); the moments are seeded to zero at the inner boundary (no dust there). After the sweep the moments are de-normalised by \(n_{\langle H\rangle}\) and floored to avoid underflow in the downstream RT.

Outputs and dust opacity

For each layer the module stores the (depleted) nucleation rate, growth timescale, grain number density, mean grain radius, and the degree of condensation, written to the dust output file. The grain optical properties are then obtained from Mie theory (LX-MIE) using the material’s complex refractive index (Opacities (transport coefficients)); calcTransportCoefficients() returns the per-layer dust absorption and scattering coefficients that enter the total opacity.

Coupling to the hydrodynamics

The frozen kernels needed by the Henyey structure solver are exposed through the agb::DustSpecies accessors nucleationRate(), growthTimescale() and secondMoment(), handed to the hydrodynamics via agb::Hydrodynamics::setDustState() (a no-op for the default shooting path).