Radiative transfer¶

Module: src/radiative_transfer/ — agb::RadiativeTransfer, agb::RadiationField, agb::ImpactParam.

The radiative transfer is solved with the Rybicki–Hummer variable-Eddington- factor (VEF) moment method for a spherically symmetric, static medium, as described in the diploma thesis of D. Kitzmann. The implementation matches that scheme term by term.

Method overview¶

The VEF method splits the angle-dependent transfer problem into two coupled pieces that are iterated to consistency:

A formal solution along tangent rays (impact parameters) gives the angular distribution of the specific intensity at each radius, and from it the closure quantities — the Eddington factor \(f_\nu=K_\nu/J_\nu\) and the sphericality factor \(q_\nu\).
A moment system for \(J_\nu(r)\), closed with the just-computed \(f_\nu\) and \(q_\nu\), is a cheap tridiagonal problem that yields a new mean intensity and hence a new source function.

Iterating (1) ↔ (2) converges because the Eddington factors are only weakly dependent on the source function, while the expensive angular detail is needed only to update them.

Geometry: impact parameters and tangent rays¶

For a spherical atmosphere the natural coordinates are impact parameters \(p\) (tangent rays). The set of rays comprises:

nb_core_impact_param rays that strike the stellar core (\(p<R_\star\)), carrying the inner boundary intensity, and
one tangent ray per radial shell (\(p=r_i\)), so that each shell is sampled.

Along each ray the geometric coordinate is \(z=\sqrt{r^2-p^2}\); the set of zPoint crossings of a ray with the radial shells defines the spatial grid for the Feautrier variables

\[u_\nu = \tfrac12 (I_\nu^+ + I_\nu^-), \qquad v_\nu = \tfrac12 (I_\nu^+ - I_\nu^-),\]

the mean and difference of the outgoing/incoming intensities. The ray geometry is built once and rebuilt (agb::RadiativeTransfer::rebuildGeometry()) only when the radial grid moves (movable grid).

Formal solution on a ray (Feautrier)¶

Along each impact parameter the second-order Feautrier equation for \(u_\nu\) is discretised on the ray’s optical-depth grid and solved as a tridiagonal system (agb::ImpactParam). Two discretisations are available:

Taylor finite differences (assembleSystemTaylor) — the default.
Cubic splines (assembleSystemSpline).

The difference variable \(v_\nu\) is then recovered (calcV), and angular integration over all rays passing through a shell yields the moments \(J_\nu, H_\nu, K_\nu\) and thus the closure factors at that shell.

The moment system for \(J_\nu\)¶

With \(f_\nu\) and \(q_\nu\) frozen from the formal solution, the zeroth/first moment equations combine into a single second-order ODE for \(J_\nu\). In the optical-depth-like coordinate \(X\) (built by generateXGrid from the extinction and the sphericality factor) this is

(1)¶\[\frac{\partial}{\partial X}\!\left(\frac{\partial (f_\nu q_\nu r^2 J_\nu)} {\partial X}\right) = q_\nu r^2 (J_\nu - S_\nu) ,\]

with the source function

\[S_\nu = \frac{\kappa_\nu B_\nu(T) + \sigma_\nu J_\nu}{\chi_\nu} \quad\Longrightarrow\quad \chi_\nu(J_\nu - S_\nu) = \kappa_\nu (J_\nu - B_\nu) .\]

Here \(\kappa_\nu B_\nu\) is the thermal emission; for a two-temperature medium it is the sum over gas and dust, \(\kappa_{\nu,\mathrm{gas}}B_\nu(T_\mathrm{gas})+ \kappa_{\nu,\mathrm{dust}}B_\nu(T_\mathrm{dust})\). This emission is iteration-invariant within one RT solve and is therefore precomputed once (precomputeEmission → planck_emission[nu][i]).

Eq. 1 is assembled as a tridiagonal system (assembleMomentSystemTaylor / …Spline) and solved with the Thomas algorithm (agb::aux::TriDiagonalMatrix). Iterating with the scattering term \(\sigma_\nu J_\nu\) updated each pass converges the field; solveMomentSystem drives this for every frequency.

Taylor vs. spline discretisation¶

The two discretisations differ in their stability:

Taylor finite differences keep the system matrix an M-matrix (positivity-preserving) unconditionally, so \(J_\nu\ge0\) always.
Cubic splines lose the M-matrix property once the optical depth per cell exceeds \(\sqrt 6\approx2.45\), producing negative mean intensities — exactly as the thesis (§A.1.3) warns.

Warning

Use the Taylor discretisation (use_spline_discretisation = false) for production. The spline option is retained for cross-checks only. (Historic negative-intensity crashes were ultimately traced to a density-collapse bug upstream in the hydrodynamics, not to the RT scheme: a transparent atmosphere makes both the moment system and the Feautrier solve degenerate.)

Two flux definitions and why they coexist¶

The frequency-integrated flux is needed in two roles that pull toward different discretisations, so the code keeps two integrated-flux fields:

Transport form (eq. 2.58), eddington_flux / eddington_flux_int — the per-frequency flux from the gradient of \(f_\nu q_\nu r^2 J_\nu\) (calcFlux). This is physical (non-negative) per frequency, so it is used for the emergent spectrum and as the denominator of the flux-mean extinction and other ratios.
Divergence form (eq. 2.59), eddington_flux_int_conservative — obtained by integrating

(2)¶\[\frac{\partial (r^2 H_\mathrm{int})}{\partial r} = r^2 \!\int \kappa_{\nu,\mathrm{abs}}\,(B_\nu - J_\nu)\,\mathrm{d}\nu ,\]

outward from the diffusion inner boundary condition \(r^2 H_\mathrm{int}(r_1)=r_1^2\,\mathrm{bfc}\cdot (\mathrm{d}B/\mathrm{d}T)/\chi\). Its divergence equals the local balance that the moment-equation \(J\)-solve enforces, so \(r^2 H_\mathrm{int}\) is conserved at radiative equilibrium. This is the quantity compared against the target luminosity in the flux-convergence test and used by the temperature corrector’s flux term.

Note

Keeping these separate is essential. eddington_flux_int is simultaneously a ratio denominator (which must stay consistent with its eq. 2.58 numerator) and a flux value (which wants the conservative eq. 2.59 form). Conflating them corrupts the flux-mean extinction, which feeds the wind, and makes the coupled loop diverge. The behaviour is controlled by flux_from_divergence (Configuration (config.toml)).

Outputs of the RT module¶

For each radial shell, agb::RadiationField stores the spectral moments (mean_intensity, eddington_flux, eddington_k), their wavelength integrals, the Eddington and sphericality factors, and the angle-grid intensities. It also provides the flux-weighted (flux-mean) extinction \(\chi_F\) (fluxWeightedExtinction), which is the single most important quantity handed to the hydrodynamics — it sets the radiative acceleration Eq. 1.

The linearised moment operator¶

For the full-linearisation temperature corrector the RT module can also expose its operator in linearised form (buildLinearisedMomentSystem): per frequency it returns the frozen Taylor moment matrix \(V_\nu\), the diagonal source derivative \(\mathrm{d}S/\mathrm{d}T\), and the moment-equation residual \(K_\nu=\mathrm{rhs}_\nu - V_\nu J_\nu\). The temperature corrector consumes these; see Radiative-equilibrium temperature correction.