Radiative-equilibrium temperature correction¶

Module: src/temperature_correction/ — agb::TemperatureCorrection.

The temperature structure is fixed by radiative equilibrium, not by an energy equation. Two correctors are implemented and can be combined in a two-phase scheme:

the Unsöld–Lucy correction — a robust, approximate fixed-point update;
the full-linearisation (Newton) correction — an exact-Jacobian step that removes the Unsöld–Lucy accuracy floor.

Both correct the gas and the dust temperatures independently, each toward its own local radiative-equilibrium condition Eq. 2.

Why two temperatures¶

Gas and dust exchange energy with the radiation field through different opacities (molecular/atomic lines and continua for the gas; the amorphous-carbon Mie opacity for the dust) and are not collisionally locked in the thin outer wind. Each component therefore has its own radiative-equilibrium temperature: the energy-balance condition Eq. 2 is solved separately for \(T_\mathrm{gas}\) (with \(\kappa_{\nu,\mathrm{gas}}\)) and \(T_\mathrm{dust}\) (with \(\kappa_{\nu,\mathrm{dust}}\)), sharing the same radiation field \(J_\nu\).

Unsöld–Lucy correction¶

Method: agb::TemperatureCorrection::calculate().

The Unsöld–Lucy scheme produces a Planck-integrated correction \(\Delta B\) per layer from a combination of the local energy-balance residual and the non-local flux-constancy residual, weighted by the absorption-, Planck- and flux-mean opacities (calcIntegratedQuantities assembles \(f q J\), \(q\chi H\), \(\kappa_J\), \(\kappa_B\) and the integrated Planck function). The correction is converted to a temperature change through an exact \(T^4\) (Planck) update rather than a linearised one.

The scheme is only first-order accurate and leaves a residual accuracy floor (a few × \(10^{-3}\) in the energy balance, ~2 % in flux conservation for this problem). It is, however, very robust far from convergence, which is exactly where it is used.

The decisive subtlety: the ratio form of the constraint¶

Both correctors enforce local radiative equilibrium in ratio form (thesis 3.40/3.41), per species \(s\):

(1)¶\[g_s = \frac{\displaystyle\sum_\nu w_\nu\,\kappa_{\nu,s}\,J_\nu} {\displaystyle\sum_\nu w_\nu\,\kappa_{\nu,s}\,B_\nu(T_s)} - 1 = 0 ,\]

not the difference form \(\sum_\nu w_\nu\kappa_{\nu,s}(J_\nu-B_\nu)\).

Warning

This choice is essential for the linearisation. The difference-form integrals span many decades with depth, so the dense Newton matrix becomes catastrophically ill-scaled (deep rows huge, thin rows tiny) and produces garbage — uniform unresolved residuals, wild dust-front swings, and a flux that bulges in the interior. Dividing each species row by its denominator \(\mathrm{den}_{s,i}=\sum_\nu w_\nu\kappa_{\nu,s}B_\nu\) makes every row :math:`O(1)`, and the Newton solve is well conditioned.

Full-linearisation (Newton) correction¶

Method: agb::TemperatureCorrection::linearisedCorrection() (temperature_correction_linearisation.cpp), with the RT half supplied by agb::RadiativeTransfer::buildLinearisedMomentSystem().

One Newton step is taken on \((T_\mathrm{gas}, T_\mathrm{dust})\) with the converged RT operator frozen (Eddington/sphericality factors, opacities, geometry). The key point is that the intensity response is the exact moment- system response, not the \(\tau\)-weighted Unsöld–Lucy heuristic.

Per frequency, the moment system \(V_\nu J_\nu = \mathrm{rhs}_\nu(T)\) is re-assembled with the same Taylor stencil used in the RT solve, so the intensity response to a temperature change is

\[\delta J_\nu = V_\nu^{-1}\bigl(K_\nu + S_g\,\delta T_g + S_d\,\delta T_d\bigr), \qquad S_{s,i} = \frac{\partial\,\mathrm{rhs}_i}{\partial T_{s,i}} = -c_i\,\kappa_{\mathrm{abs},s}\,\frac{\mathrm{d}B}{\mathrm{d}T},\]

with \(c_i=r_i^2/(q_i\chi_i)\) and the residual \(K_\nu=\mathrm{rhs}_\nu - V_\nu J_\nu\). The diagonal (in the Taylor discretisation) source derivatives make \(S_{s}\) cheap.

The mean intensities are eliminated Rybicki-style — solving \(V_\nu\) against each unit node vector gives the inverse columns — collapsing the problem to a dense \(2D\times2D\) temperature system (\(D=\) number of grid points, factor 2 for gas+dust), solved by Gaussian elimination.

Note

The frequency reduction into the dense Newton matrix uses a deterministic (static-scheduled, ascending-thread-order) reduction. An earlier arrival-order reduction summed the ~5000 frequencies in a different order each run, perturbing the matrix at the \(10^{-15}\) level; amplified by the stiff dust-nucleation coupling this changed the iteration count run-to-run (same solution). The reduction is now bit-reproducible at a fixed thread count.

Local RE alone is not enough: the flux-constancy term¶

A frozen-structure test revealed that driving the local RE residual Eq. 1 to zero does not conserve the flux: in the optically thick interior local RE is degenerate (\(J\approx B\) regardless of \(T\)), so the Newton converges without pinning the flux (\(r^2 H\) bulged by ~13 %).

The fix is the composite constraint (thesis §3.2.4, eq. 3.64/3.65): per species,

(2)¶\[\xi\,\underbrace{g_{\text{local},s}}_{\text{local-RE ratio}} \;+\; \zeta(r)\,\underbrace{g_{\text{flux},s}}_{\text{flux constancy}} = 0 ,\]

with the flux-constancy residual written in the integrated (eq. 2.59) form

\[g_{\text{flux},i} = \frac{1}{r_\star^2 H_\star} \int_{r_1}^{r_i} r'^2\,\Bigl[\textstyle\sum_\nu w_\nu\bigl( \kappa_g (J_\nu-B_g) + \kappa_d (J_\nu-B_d)\bigr)\Bigr]\,\mathrm{d}r' .\]

Because the divergence of this integral is exactly the local-balance source that solveMomentSystem discretises (eq. 2.61 = eq. 2.59 in \(X\)-coordinates), the flux term is consistent with the \(J\)-solve and does not fight local RE — and the radial integral with the \(1/(r_\star^2 H_\star)\) normalisation keeps it \(O(1)\)-scaled.

The depth weight is \(\zeta(r)=\tau/(\tau+\tau_\mathrm{scale})\), built from a grey radial optical depth (flux-mean extinction integrated inward): \(\zeta \to1\) deep where flux constancy must dominate, \(\zeta\to0\) at the thin outer edge where the flux derivative is numerically noisy (and forced to 0 at the two boundary nodes). The weights are linearisation_xi (\(\xi\), ≈1) and linearisation_zeta_tau_scale.

Important

Only the integrated (eq. 2.59) flux operator works. Linearising the node-centred eq. 2.58 (transport) flux fails — its discrete divergence is inconsistent with the \(J\)-solve stencil and the composite Jacobian cancels where the two constraints pull oppositely (determinant → 0, divergence). The flux term is applied on the gas row only (dust keeps pure local RE) to avoid gas/dust collinearity.

Two-phase corrector and coupled-loop stability¶

Driver: agb::AGBStarModel::temperatureIteration().

The frozen-opacity Newton step over-corrects far from convergence — the dust nucleation rate \(J_\star\sim\exp(-/T)\) is so temperature-sensitive that a frozen-opacity step badly mispredicts the dust response. Two mechanisms keep the coupled iteration stable:

Two-phase start. Unless told otherwise (linearisation_start_unsoeld_lucy), the iteration begins with Unsöld–Lucy and latches onto the linearisation only once, for linearisation_switch_count consecutive steps, both the per-step max|dT|/T has settled below linearisation_switch_dt_fraction × max_relative_change and the maximum RE residual is below linearisation_switch_re_residual. The switch is one-way (no revert).
Per-layer adaptive under-relaxation of the Newton step. In the coupled loop the structure (dust, hydro) updates between temperature steps, so a full Newton step can drive a period-2 ±cap limit cycle just beyond the dust front. The Newton \(\Delta T\) is damped with the same adaptive scheme as Unsöld–Lucy: a per-layer factor \(\omega\) halved on a sign flip and grown back when the step keeps its sign, applied after the per-step cap so the damping is visible even while the cap binds. A global linearisation_relaxation (default 0.5) further under-relaxes the step.

The step order in temperatureIteration() is, deliberately:

compute the raw correction (UL or linearisation);
(linearisation) per-layer adaptive damping;
(UL only, and only in the settling regime) Anderson acceleration;
optional monotonicity clip;
cap \(|\Delta T|\le\) max_relative_change \(\times T\) last, and floor \(T>0\).

Capping last guarantees that nothing (an Anderson overshoot, a monotonicity clip) produces a step larger than the bound the hydro–dust cycle tolerates.

Note

Two historical oscillation bugs are worth recording, as they recur if the ordering is disturbed:

force_monotonic fought the correction (the clip cascaded the just-applied heating back down every step) → ±cap limit cycle. Default off.
The under-relaxation \(\omega\) was applied before the cap, so while the cap was binding the step stayed pinned at ±cap regardless of \(\omega\) and an oscillating layer never damped. Fix: cap first, then multiply by \(\omega\).

Anderson acceleration¶

For the Unsöld–Lucy phase, Anderson mixing (agb::AGBStarModel::andersonStep()) accelerates the slow creep toward radiative equilibrium. It is gated to the settling regime only (prev_max_rel_change < anderson_activation_fraction × max_relative_change), so the early large-change transient — to which the hydro cycle is sensitive — stays on plain damped steps, and the history is restarted fresh on activation so it is not seeded by the transient. Anderson is bypassed once the linearisation is active (Newton already mixes).

Accuracy achieved¶

With the integrated flux term, the linearisation corrector reaches flux conservation at the \(\sim10^{-5}\) level on a frozen structure (versus the ~2.3 % Unsöld–Lucy floor) and converges the coupled IRC+10216 model in ~20 outer iterations from a warm start.