Analytic temperature Jacobians

For thermal (longwave) problems the C++ backend can return the analytic derivatives of every output field with respect to every temperature degree of freedom, at a small fraction of the cost of the forward solve. This is enabled by compute_temperature_jacobian and follows the linearisation of the interaction principle of Spurr & Christi (2006), specialised to the temperature variable.

What is differentiated

Four output fields are differentiated:

\[\frac{\partial F^{\uparrow}_k}{\partial T_m},\quad \frac{\partial F^{\downarrow}_k}{\partial T_m},\quad \frac{\partial J_k}{\partial T_m},\quad \frac{\partial (\nabla\!\cdot\!F)_k}{\partial T_m},\]

for the upward flux, downward flux, actinic mean intensity, and flux divergence (heating rate). They fill the four RTOutput arrays flux_up_temperature_jac, flux_down_temperature_jac, mean_intensity_temperature_jac, and flux_divergence_temperature_jac.

Each array is indexed [interface][dof]:

interface \(k = 0, \dots, N_{\text{lyr}}\) — the output level.
dof \(m\) — the temperature it is differentiated against: dof 0 .. num_layers are the level temperatures \(\partial/\partial T_m\); dof num_layers+1 is the independent surface skin temperature \(\partial/\partial T_s\).

When the solver is driven by planck_levels instead of temperatures, the arrays hold \(\partial/\partial B_m\) (the Planck chain rule is skipped).

Why it is cheap and exact

Temperature enters the adding–doubling equations only through the Planck emission source. Every reflection and transmission operator is temperature-independent,

\[\frac{\partial \mat{R}}{\partial B_m} = \frac{\partial \mat{T}}{\partial B_m} = \mat{0}.\]

The Jacobian is therefore obtained by re-running the same forward linear operators with differentiated right-hand sides. It is exact to machine precision within the discrete-ordinate approximation itself. The direct solar beam carries no temperature derivative, so the solar source vectors are ignored throughout the Jacobian pass.

The implementation works in two stages: it first forms derivatives with respect to the level Planck values \(B_m\), then applies a single per-column Planck factor \(\dd B_m / \dd T_m\) at the end.

Source seeds

Each layer’s thermal source is written compactly as

\[\vect{s}^{\uparrow}_l = \vect{a}_l\,\bar B_l + \vect{b}_l\,B'_l, \qquad \vect{s}^{\downarrow}_l = \vect{a}_l\,\bar B_l - \vect{b}_l\,B'_l,\]

with layer mean \(\bar B_l = \tfrac12(B_l + B_{l+1})\) and gradient \(B'_l = (B_{l+1} - B_l)/\tau_l\). The coefficient vectors \(\vect{a}_l, \vect{b}_l\) are temperature-independent (they are the doubling vectors \(\vect{y}_l, \vect{z}_l\) in the scattering branch, and the closed-form absorption coefficients otherwise). Because \(\bar B_l\) and \(B'_l\) are linear in the two bounding level values, defining

\[\vect{p}_l \equiv \tfrac12\vect{a}_l - \frac{\vect{b}_l}{\tau_l}, \qquad \vect{q}_l \equiv \tfrac12\vect{a}_l + \frac{\vect{b}_l}{\tau_l},\]

gives the complete layer-level seed as just two \(N\)-vectors:

\[\frac{\partial\vect{s}^{\uparrow}_l}{\partial B_l} = \vect{p}_l,\quad \frac{\partial\vect{s}^{\uparrow}_l}{\partial B_{l+1}} = \vect{q}_l,\quad \frac{\partial\vect{s}^{\downarrow}_l}{\partial B_l} = \vect{q}_l,\quad \frac{\partial\vect{s}^{\downarrow}_l}{\partial B_{l+1}} = \vect{p}_l.\]

Layer \(l\) seeds only columns \(m = l\) and \(m = l+1\). The surface seed is \(\partial\vect{s}^{\uparrow}/\partial B_s = (1-\alpha)\vect{1}\); the boundary intensities contribute \(\partial\vect{I}^{\downarrow}_{\text{top}}/\partial B_0 = \vect{1}\) and, for the diffusion lower boundary, \(\partial I^{\uparrow}_{\text{bot},i}/\partial B_L = 1 + \mu_i/\tau_{L-1}\) and \(\partial I^{\uparrow}_{\text{bot},i}/\partial B_{L-1} = -\mu_i/\tau_{L-1}\).

Propagation through adding

The source-combination law (3) is linear in the sources with temperature-independent operators, so it applies verbatim to the derivatives:

\[\begin{split}\frac{\partial\vect{S}^{+}}{\partial B_m} &= \frac{\partial\vect{s}^{+}_1}{\partial B_m} + \mat{T}_{ba}\mat{D}_1\!\left( \frac{\partial\vect{s}^{+}_2}{\partial B_m} + \mat{R}_{bc}\frac{\partial\vect{s}^{-}_1}{\partial B_m}\right), \\ \frac{\partial\vect{S}^{-}}{\partial B_m} &= \frac{\partial\vect{s}^{-}_2}{\partial B_m} + \mat{T}_{bc}\mat{D}_2\!\left( \frac{\partial\vect{s}^{-}_1}{\partial B_m} + \mat{R}_{ba}\frac{\partial\vect{s}^{+}_2}{\partial B_m}\right).\end{split}\]

All \(N_{\text{lyr}}+2\) derivative columns are propagated together as an \(N \times (N_{\text{lyr}}+2)\) block, so each operator application is a single matrix product against all degrees of freedom at once. One bottom-up sweep produces the source derivatives of every rbase composite and one top-down sweep every rtop composite.

At each interface, the intensity solve \((\mat{I} - \mat{R}\mat{R})\,\vect{I} = \text{rhs}\) uses a temperature-independent matrix, so the forward LU factorisation is cached and reused to back-substitute \(\partial\vect{I}/\partial B_m\) against \(\partial(\text{rhs})/\partial B_m\). The flux and mean-intensity derivatives are then the same quadrature reductions as the forward pass.

Heating-rate Jacobian

Differentiating (6) gives

(13)\[\frac{\partial(\nabla\!\cdot\!F)_k}{\partial B_m} = 4\pi\,(1-\omega_{l(k)})\left(\frac{\partial J_k}{\partial B_m} - \delta_{km}\right),\]

i.e. the mean-intensity Jacobian scaled by \(4\pi(1-\omega_{l(k)})\), with a diagonal Planck-cooling correction \(-4\pi(1-\omega_{l(k)})\) on the level-\(k\) column (from \(\partial B_k/\partial B_m = \delta_{km}\)). The surface column carries no diagonal term because \(B_k\) is a level quantity.

Planck chain rule

Converting from \(B\)-Jacobians to \(T\)-Jacobians is a single per-column scaling,

\[\frac{\partial(\cdot)}{\partial T_m} = \frac{\partial(\cdot)}{\partial B_m}\,\frac{\dd B_m}{\dd T_m},\]

using the band-integrated Planck temperature derivative

\[\frac{\dd B}{\dd T} = \frac{4B}{T} - \frac{\sigma}{\pi}\frac{15}{\pi^4}\,T^3 \bigl(v_1 f(v_1) - v_0 f(v_0)\bigr), \qquad v = \frac{c_2 \nu}{T},\ f(v) = \frac{v^3}{\ee^v - 1}.\]

Level columns use \(\dd B_m/\dd T_m\) at temperature[m]; the surface column uses \(\dd B_s/\dd T_s\) at surface_temperature.

Depth limit of the heating Jacobian

The heating Jacobian \(\mat{H} = \partial(\nabla\!\cdot\!F)/\partial T\) is built from the even (zeroth) angular moment \(J\) and carries a local Planck-cooling diagonal, which makes it diagonally dominant in the mid optical-depth range. This is in deliberate contrast to the net-flux Jacobian \(\mat{K} = \partial F/\partial T\), built from the odd (first) moment, whose diagonal cancels by symmetry (it is identically zero on a uniform optical-depth grid) and which is therefore dense and long-ranged. When a temperature corrector needs a well-conditioned operator, it should drive the collocated heating residual \(\mat{H}\), not the flux residual.

However, \(\mat{H}\) degenerates at large optical depth. Using the linear-in-\(\tau\) source, the self-response of the mean intensity to the local level temperature has the closed form

\[\frac{\partial J_i}{\partial B_i} = \tfrac12\bigl[G(a) + G(b)\bigr], \qquad G(c) = 1 - \frac{1}{2c} + \frac{E_3(c)}{c},\]

where \(a\) and \(b\) are the optical thicknesses of the layers just above and below level \(i\), and \(E_n\) are the exponential integrals. In the two limits,

\(G \to 0\) as \(c \to 0\) — optically thin: \(J\) is insensitive to the local temperature (the “top collapse”);
\(G \to 1\) as \(c \to \infty\) — optically thick: \(J\) is slaved to the local Planck function.

The heating diagonal is therefore

\[H_{ii} = 4\pi(1-\omega)\,\frac{\partial B_i}{\partial T_i} \left(\tfrac12[G(a)+G(b)] - 1\right) \;\xrightarrow[\text{deep}]{}\; -\frac{2\pi(1-\omega)}{\Delta\tau}\,\frac{\partial B}{\partial T} \to 0.\]

The Planck-cooling diagonal is cancelled by the self-response \(\partial J_i/\partial B_i \to 1\), and the whole operator collapses as \(1/(2\Delta\tau)\). The crossover is governed by the local layer optical thickness, not the cumulative optical depth: by \(\Delta\tau \approx 3\) the diagonal is roughly 16 % of its ceiling, and by \(\Delta\tau \approx 10\) the heating Jacobian is effectively zero.

Note

This is transport physics, not a numerical defect: at large optical depth the field is in local thermodynamic equilibrium with the source (\(J \to B\)), so there is no local radiative restoring force on temperature. A radiative-equilibrium temperature corrector must combine the heating Jacobian (which “owns the middle of the optical-depth range”) with a deep, flux-based diffusion treatment and a thin-top Planck inversion.

Cost, caching, and scope

No forward work is repeated. The one-time \(O(N_{\text{lyr}} N^3)\) doubling and adding costs are reused. The extra work is two source-only adding sweeps plus one differentiated output pass, each \(O(N_{\text{lyr}}^2 N^2)\) for the full block — typically a single-digit percentage on top of the forward solve.
Caching. On the templated CPU path the forward build captures the LU factorisations of the adding matrices and interface matrices, so the Jacobian never re-factorises.
Validation. The feature is checked against central finite differences to \(\sim 10^{-3}\) relative error and cross-checked against the sibling DisORT++ solver to \(\sim 10^{-5}\), over the absorption and scattering branches, both boundary conditions, both surface modes, and delta-M on/off.
Scope. Implemented on the templated CPU path (with caching) and the dynamic fallback (Eigen-based, no caching). A CUDA port is deferred. Only temperature (Planck-source) Jacobians are provided; optical-depth, SSA, and phase-function Jacobians are out of scope.