Radiative Transfer
The step-by-step equation of ARTS radiative transfer calculation is described in this section.
The key variables to look out for are spectral_rad and spectral_rad_jac.
A single step
A single step of the radiative transfer calculation in ARTS solvers is assumed to follow from
where \(I_{1}\) is the spectral radiance after the step, \(I_{0}\) is the spectral radiance before the step, \(J_{0}\) is the spectral radiance source at the step, and \(T_{0}\) is the transmittance through the current step.
To get more of an intuitive understanding of this equation, a “step” could be considered as a single layer of the atmosphere. Imagine that the radiation is coming into the layer at some level is \(I_{0}\), and that the radiation is leaving at another level is \(I_{1}\). The radiation that is leaving is the sum of the radiation that is coming in, and the transmitted incoming radiation differences to the local source of the radiation \(T_{0} \left(I_{0} - J_{0}\right)\), and the radiation that is emitted by the layer itself \(J_{0}\)
The way to compute \(J_{0}\) and \(T_{0}\) is not discussed here. What is important is that they are computed based on the state of the model in the surrounding levels of the model. Noting that the radiation coming into the layer may also have been reflected and can thus depend on all model states, we can write this as
where \(\vec{x}_0\), \(\vec{x}_1\) is the state of the levels above and below the layer, and \(\mathbf{x}\) is the state of the model everywhere. Both \(\vec{x}_0\) and \(\vec{x}_1\) are covered by \(\mathbf{x}\), though they are not necessarily in \(\mathbf{x}\), but can be derived from it, e.g., via interpolation. This will matter when we compute the Jacobian - the Jacobian is on \(\mathbf{x}\) but radiative transfer quantities are on \(\vec{x}_0\) and \(\vec{x}_1\), etc, so there must always be a map back from \(\vec{x}_0\) and \(\vec{x}_1\), etc, to \(\mathbf{x}\).
Note
The notation of ARTS often calls the content of \(\vec{x}_i\) as targets of the Jacobian and a partially flattened version of \(\mathbf{x}\) as the state vector.
The unflattened dimension being the measurement dimension, e.g., frequency or channels.
The concept of targets comes from the fact that the Jacobian is only computed for user-selected variables in \(\mathbf{x}\). These selected variables are thus the users’ “target”.
Multiple steps
We are not always interested in the radiation at level 1, but quite often in the radiation at the end of the atmospheric path. This, and the radiation at intermittent levels, can be computed by simply taking multiple steps of the above equation. Thus, going from level 0 to level N, we get
Both \(J_i\) and \(T_i\) depend on values at levels \(\vec{x}_i\) and \(\vec{x}_{i+1}\), even though these have been omitted for brevity. There are here \(N-1\) steps or layers, and \(N\) levels. Note that these steps are non-commutative, since the step through a layer polarizes the radiation.
Note
\(I_0\) is often called the background radiation, as it is the radiation entering the atmosphere from space, from the ground, or elsewhere. The background, it other words.
If we expand the above expressions and collect them by transmittance, we get
which we can rewrite as a simple sum of terms:
and so on. Here, \(\vec{J}'\) is the source vector, but we use a prime-tick because it is 0-indexed whereas the resulting \(\vec{I}\) is 1-indexed.
In practice, it is often more convenient to just use the step-by-step equations for forward calculations, but the matrix above gives insight on optimizations for the Jacobian of the “measured” \(I_N\).
Propagation-matrix exponential
For polarized radiative transfer, \(K\) and \(T\) are matrices and the transmittance is evaluated as a matrix exponential. Both the ordinary and Magnus layer options use the same specialized Cayley–Hamilton evaluator; they differ only in the matrix supplied to that evaluator.
After separating the scalar diagonal part, the characteristic polynomial of the polarized part \(P\) has the form
With \(S=\sqrt{B^2-4C}\), define
The four eigenvalues of \(P\) are then \(+x\), \(-x\), \(+iy\), and \(-iy\). In particular, the expressions above already give \(x^2\) and \(y^2\); the square roots used to obtain \(x\) and \(y\) are taken only once afterward. The exponential can consequently be reduced to
where the coefficients contain combinations of \(\cosh(x)\),
\(\sinh(x)\), \(\cos(y)\), and \(\sin(y)\). ARTS evaluates the
zero and repeated-eigenvalue limits separately to avoid divisions by small
numbers. The derivatives of the exponential coefficients and their
degenerate limits are evaluated analytically and are used by all
rte_option choices described below.
Partial derivatives
The Jacobian after the last step is important for the retrieval of the state of the atmosphere. It is in this formalism defined as
We now introduce an alternative notation to make the expressions below more compact,
where for sake of keeping the expressions less cluttered, we add \(T_N=1\).
We can summarize the full Jacobian after N steps as
The first term is the background radiation contribution, which is lifted from the rest of the expression because we can then simplify the rest of the expression. Since we know that we can map from \(\vec{x}_i\) to \(\mathbf{x}\), we can write the partial derivatives in terms of \(\vec{x}_i\) instead of \(\mathbf{x}\). The source term then becomes
and for the transmittance, the term becomes
The expression is then mapped back to \(\mathbf{x}\) in the following:
where the last term is 0 for all but \(i=N\) and \(i=N-1\) and where the function \(f\) is defined as the map from \(\vec{x}_i\rightarrow\mathbf{x}\). Often in ARTS, \(f\) is just an inverse interpolation operator.
Layer interpolation options
Note that here, the indexing is not the same as above as
it is layer-oriented notation. ARTS provides several approximations for
how the source and propagation matrix vary between the two endpoints of a
layer. They are selected through rte_option.
The option names and their selection are described in Radiative-transfer layer options. Their mathematical definitions follow below.
Linear source function
This is based on the assumption of a linearly changing source function and constant propagation matrix through a layer. The indexing is the same for the transmittance and incoming background radiation as above - in layers. However, the source function indexing below is now shifted to represent levels - the same indexing as we use for the spectral radiance term.
Linear in source means that we assume \(J = J_0 + (J_1 - J_0) \frac{r}{r_x}\), where \(r\) is the distance through the layer and \(r_x\) is the total thickness of the layer. With constant propagation matrix \(K\), we define \(T_0 = \exp(-K_0 r_x)\) and get
or
This is a classical solution to the radiative transfer equation to improve numerical stability when the propagation matrix \(K\) is large. Going through the motion of extending this to multiple steps as for the constant source case is then just a matter of adding 1 extra term per layer in the dot products above. If we define
we get
and again, we can expand these expressions and collect them by transmittance, but this time with an additional term per layer
Which in shortened matrix form is:
The \(\vec{J}\) term is now without a prime-tick as it is 1-indexed like \(\vec{I}\), though \(J_0\) is still used in the RHS vector.
The partial derivative propagation is then derivable in the same way as above,
And again doing the expansion of the dot products and using \(T_N=1\) for brevity, we get for the source partial derivatives at the end of the path:
Following the same procedure as above, we can rewrite this in terms of \(\vec{x}_i\) instead of \(\mathbf{x}\). Here, a key difference from before is that the source term now only depends on \(\vec{x}_i\) and not on \(\vec{x}_{i+1}\).
The derivative contribution is then given by
Linear propagation matrix and source function
The linprop option assumes that both propagation and source vary linearly
through a layer. Let endpoint 0 be the incoming/background end, endpoint 1
the outgoing/observer end, and let \(0\leq s\leq r\) measure distance
between them:
The layer step has the same affine form as for lintau,
For scalar propagation the endpoint-average transmission
is exact, and ARTS evaluates the exact source operator
Equivalently, define the dimensionless average optical depth and endpoint difference
Then the same operator is
This form is evaluated without numerical quadrature: ARTS uses a convergent
endpoint-difference expansion near \(\delta=0\), stable Dawson or scaled
complementary-error-function forms for well-separated endpoints, and an
endpoint asymptotic expansion for sufficiently thick layers. Both partial
derivatives of the integral are evaluated analytically. At
\(K_0=K_1\), it reduces continuously to the lintau operator
\(\phi_1(-rK_0)\).
The scalar completion-of-the-square identity does not extend element by
element to a polarized propagation matrix. For polarized linprop, ARTS
instead defines
where \(\phi_1(X)=X^{-1}(\exp(X)-\mathbf{I})\) is evaluated through its regular power-series limit at singular arguments, \(\kappa_i\) is the scalar (diagonal) component of \(K_i\), and \(g_{\rm av}=-r(\kappa_0+\kappa_1)/2\). The scalar correction \(c\) makes this construction reduce exactly and continuously to the scalar linear-propagation integral above as polarization tends to zero. The remaining term is a commutator-free augmented-source approximation: the source-gradient correction \(Q\) is retained, but the transmission is still the ordinary endpoint-average exponential. The polarized derivatives of both \(T_{\rm av}\) and \(\Lambda_{\rm pol}\) are analytical.
Magnus propagation matrix
The Magnus options treat the propagation matrix as linear across a layer,
and retain the first commutator correction in the ordered exponential. With \([A,B]=AB-BA\), the effective optical-depth matrix is
When \(K_0=K_1\), or when the endpoint matrices commute, the commutator
vanishes and this reduces to the ordinary exponential of the endpoint
average. The commutator term is the distinction between magop and the
constant transmittance approximation.
The magop option uses the endpoint-average source
This option improves the ordered propagation matrix but does not model the source gradient explicitly.
Magnus propagation matrix and linear source
The magop_linsrc option augments the Magnus system with a source that is
linear between \(J_0\) and \(J_1\). Define the matrix function
where its power-series limit is used when \(X\) is singular. The linear-source operator used by ARTS is
and the layer step is
For a constant propagation matrix this reduces to the lintau source
operator. The polarized linprop and magop_linsrc options use the
same augmented-source factor
\(\mathbf{I}-r(K_1-K_0)/12\), with the scalar completion \(c\) above
added for linprop. The linprop option evaluates \(\phi_1\) at
the average generator \(G_{\rm av}\) and uses
\(\exp(G_{\rm av})\) for transmission. magop_linsrc instead
evaluates it at \(-\Omega\) and uses \(\exp(-\Omega)\), thereby
including the first ordering commutator in both operators. For varying
scalar propagation, linprop uses the exact integral above, whereas
magop_linsrc remains the truncated augmented-Magnus approximation.
The derivatives of both \(T\) and \(L\) with respect to either
endpoint and the layer length are evaluated analytically.
Both Magnus options truncate the Magnus series after the first commutator term. They therefore still require sufficiently short layers when the propagation matrix changes rapidly; reducing the ray-path step size is the appropriate convergence check in that regime.
Deriving the expressions
This is not an important section for users of ARTS, but is included for completeness to show how the expressions above are derived. It may also be useful if you want to implement other variations of the radiative transfer equation solution.
The integrating-factor derivations below apply directly to scalar transfer. They also apply to the constant-matrix cases because all functions of the single matrix \(K\) commute. For a varying matrix \(K(s)\), the derivation applies only when the matrices commute along the path; the polarized approximations used by ARTS are described above. The path variable is \(s \in [0,r]\) for the distance through a layer of understood physics, with \(I(0)=I_0\) and \(I(r)=I_1\).
The differential equation in all cases is
The integrating factor is
and then
with
All three cases are just different assumptions for \(K(s)\) and \(J(s)\).
Constant \(K\), constant \(J\).
Take
\[\begin{split}K(s)&=K_0 \\ J(s)&=J_0 \\ T_0 &= e^{-K_0 r}.\end{split}\]Then
\[\begin{split}\mu(s) &= e^{K_0 s} \\ T(r) &=e^{-K_0 r} \\ &=T_0.\end{split}\]The solution:
\[\begin{split}I(r) &= T_0 I_0 + T_0 \int_0^r e^{K_0 s}K_0J_0ds \\ &= T_0 I_0 + T_0 \left(\int_0^r e^{K_0 s}ds\right) K_0 J_0 \\ &= T_0 I_0 + T_0 \left(e^{K_0 r} - 1\right) K_0^{-1} K_0 J_0 \\ &= T_0 I_0 + T_0 \left(e^{K_0 r} - 1\right) J_0 \\ &= T_0 I_0 + \left(1 - T_0\right) J_0.\end{split}\]Writing \(I_1 = I(r)\), we get exactly:
\[I_1 = J_0 + T_0 (I_0 - J_0),\]which is the form used above. Note that no parts of this derivation depends on whether or not the expression is scalar or matrix.
Constant \(K\), linearly varying \(J\). Assume
\[\begin{split}K(s)&=K_0 \\ J(s)&=J_0 + \frac{s}{r}(J_1 - J_0).\end{split}\]Again \(T_0 = e^{-K_0 r}\), and \(\mu(s)=e^{K_0 s}\).
Use the general solution:
\[I(r) = T_0 I_0 + T_0 \int_0^r e^{K_0 s} K_0 J(s)ds.\]Insert \(J(s)\) and split:
\[\begin{split}I(r) &= T_0 I_0 + T_0 \int_0^r e^{K_0 s} K_0 \Bigl[J_0 + \frac{s}{r}(J_1-J_0)\Bigr] ds \\ &= T_0 I_0 + T_0 \left(\int_0^r e^{K_0 s} ds\right) K_0 J_0 + T_0 \left(\int_0^r s e^{K_0 s} ds\right) K_0 \frac{J_1-J_0}{r}.\end{split}\]Compute the two integrals:
\[\begin{split}\int_0^r e^{K_0 s} ds &= \left(e^{K_0 r}-1\right) K_0^{-1} \\ \int_0^r s e^{K_0 s} ds &= \Bigl[se^{K_0 s}K_0^{-1}\Bigr]_0^r - \int_0^r e^{K_0 s} K_0^{-1} ds\\ &= r e^{K_0 r} K_0^{-1} - \left(e^{K_0 r}-1\right)K_0^{-2}.\end{split}\]Insert:
\[I(r) = T_0 I_0 + T_0 \left(e^{K_0 r}-1\right)K_0^{-1} K_0 J_0 + T_0 \Bigl[r e^{K_0 r}K_0^{-1} - \left(e^{K_0 r}-1\right){K_0^{-2}}\Bigr] K_0 \frac{J_1-J_0}{r}.\]Simplify each term, using \(T_0=e^{-K_0 r}\) and cancelling where possible.
First integral term:
\[\begin{split}T_0 \left(e^{K_0 r}-1\right)K_0^{-1} K_0 J_0 &= T_0 \left(e^{K_0 r}-1\right) J_0 \\ &= \left( 1-T_0 \right) J_0.\end{split}\]Second integral term:
\[\begin{split}T_0 \Bigl[r e^{K_0 r}K_0^{-1} - \left(e^{K_0 r}-1\right)K_0^{-2}\Bigr] K_0 \frac{J_1-J_0}{r} &= \left[ T_0 r e^{K_0 r} + T_0 \left(e^{K_0 r}-1\right)K_0^{-1} \right] \frac{J_1-J_0}{r} \\ &= \left[ r - \left(1-T_0\right)K_0^{-1} \right] \frac{J_1-J_0}{r}.\end{split}\]So
\[I(r) = T_0 I_0 + \left( 1-T_0 \right) J_0 +\left[r - \left( 1-T_0 \right)K_0^{-1}\right] \frac{J_1-J_0}{r}.\]Rearrange as
\[I_1 = J_1 + T_0 (I_0 - J_0) + \underbrace{\frac{1}{r} K_0^{-1}(1-T_0)}_{\displaystyle \Lambda_0}(J_0 - J_1),\]which gives the form used above
\[\begin{split}I_1 &= J_1 + T_0 (I_0 - J_0) + \Lambda_0 (J_0 - J_1),\\ \Lambda_0 &= \frac{1}{r} K_0^{-1} \left( 1-T_0 \right) \\ &=-\left(\log T_0\right)^{-1} \left( 1-T_0 \right).\end{split}\]Again, no assumptions were made about scalar vs. matrix. The only scary step of matrix notation non-commutativity is that you have to be aware that you can write, e.g., \(A^3=A A^2=A^2 A\), which means with \(\exp(A) = \sum_{n=0}^\infty \frac{A^n}{n!}\), you can do \(\exp(A) A = \left(\sum_{n=0}^\infty \frac{A^{n}}{n!}\right) A = \sum_{n=0}^\infty \frac{A^{n+1}}{n!} = A \sum_{n=0}^\infty \frac{A^{n}}{n!} = A \exp(A)\).
See Radiative-transfer implementation for stable evaluation of this operator.
Linear \(K(s)\) and linear \(J(s)\).
Now set both \(K\) and \(J\) to be linear in \(s\) across the layer:
\[\begin{split}K(s) &= K_0 + \frac{s}{r}(K_1-K_0),\\ J(s) &= J_0 + \frac{s}{r}(J_1-J_0).\end{split}\]Let
\[\begin{split}\alpha &= \frac{K_1-K_0}{r},\\ \beta &= \frac{J_1-J_0}{r},\\ K(s) &= K_0 + \alpha s,\\ J(s) &= J_0 + \beta s.\\\end{split}\]3.1. Integrating factor and transmittance.
The integrating factor and transmittance become
\[\begin{split}\mu(s) &= \exp\Bigl(\int_0^s K(u)du\Bigr)\\ &= \exp\Bigl(\int_0^s (K_0 + \alpha u)du\Bigr) \\ &= \exp\bigl(K_0 s + \tfrac12 \alpha s^2\bigr)\\ T(r) &= \mu(r)^{-1}\\ &= \exp\bigl(-K_0 r - \tfrac12 \alpha r^2\bigr)\\ &= T_0.\end{split}\]That’s the same \(T_0\) as always, assuming that we only evaluate the expression at \(0\) and \(r\), and that the layer value of \(K\) above is the average of the endpoints (\(K_0\) and \(K_1\)).
3.2. General solution. The general solution is still
\[I(r) = T_0 I_0 + T_0 \int_0^r \mu(s) K(s) J(s)\,ds.\]Insert the forms:
\[\begin{split}\mu(s) &= e^{K_0 s + \frac12 \alpha s^2},\\ K(s) &= K_0 + \alpha s,\\ J(s) &= J_0 + \beta s\\ &= J_1 + (J_0 - J_1)\Bigl(1 - \frac{s}{r}\Bigr),\end{split}\]where in the last step we used that \(J(r)=J_1\).
Then
\[I(r) = T_0 I_0 + T_0 \int_0^r \mu(s) K(s) \Bigl[J_1 + (J_0-J_1)\Bigl(1 - \frac{s}{r}\Bigr)\Bigr]\,ds.\]Split the integral into the part proportional to \(J_1\) and the part proportional to \((J_0-J_1)\):
\[I(r) = T_0 I_0 + T_0 \left(\int_0^r \mu(s) K(s)\,ds\right) J_1 + T_0 \left( \int_0^r \mu(s) K(s) \Bigl(1 - \frac{s}{r}\Bigr)\,ds \right) (J_0-J_1).\]The first integral is the same as for a constant source \(J(s)\equiv J_1\). For any \(K(s)\), that problem has the known solution
\[I(r) = J_1 + T_0 (I_0 - J_1),\]so comparing with the integral form gives
\[\begin{split}T_0 \left(\int_0^r \mu(s) K(s)\,ds\right) J_1 &= (1 - T_0) J_1 \quad\Rightarrow\quad\\ \int_0^r \mu(s) K(s)\,ds &= T_0^{-1} - 1.\end{split}\]Insert this back:
\[I(r) = T_0 I_0 + (1-T_0) J_1 + T_0 \left( \int_0^r \mu(s) K(s) \Bigl(1 - \frac{s}{r}\Bigr)\,ds \right) (J_0-J_1).\]Rearranging the first two terms:
\[T_0 I_0 + (1-T_0) J_1 = J_1 + T_0 (I_0 - J_1) = J_1 + T_0 (I_0 - J_0) + T_0 (J_0 - J_1),\]so
\[I(r) = J_1 + T_0 (I_0 - J_0) + \Biggl[ T_0 + T_0 \int_0^r \mu(s) K(s)\Bigl(1 - \frac{s}{r}\Bigr)\,ds \Biggr] (J_0-J_1).\]This shows that the solution has the affine form
\[I_1 = J_1 + T_0 (I_0 - J_0) + \Lambda_0 (J_0 - J_1),\]with
\[\Lambda_0 = T_0 \Biggl[ 1 + \int_0^r \mu(s) K(s)\Bigl(1 - \frac{s}{r}\Bigr)\,ds \Biggr].\]To obtain a more compact expression for \(\Lambda_0\), use integration by parts on the remaining integral. Note that
\[\bigl[\mu(s)\bigl(1 - \tfrac{s}{r}\bigr)\bigr]' = \mu'(s)\Bigl(1 - \frac{s}{r}\Bigr) - \frac{1}{r}\mu(s) = \mu(s)K(s)\Bigl(1 - \frac{s}{r}\Bigr) - \frac{1}{r}\mu(s),\]hence
\[\mu(s)K(s)\Bigl(1 - \frac{s}{r}\Bigr) = \bigl[\mu(s)\bigl(1 - \tfrac{s}{r}\bigr)\bigr]' + \frac{1}{r}\mu(s).\]Integrating this from \(0\) to \(r\) gives
\[\begin{split}\int_0^r \mu(s)K(s)\Bigl(1 - \frac{s}{r}\Bigr) ds &= \Bigl[\mu(s)\Bigl(1 - \frac{s}{r}\Bigr)\Bigr]_0^r + \frac{1}{r}\int_0^r \mu(s)\,ds \\ &= -1 + \frac{1}{r}\int_0^r \mu(s)\,ds,\end{split}\]because \(\mu(0)=1\) and \(1-r/r=0\). Inserting this into the expression for \(\Lambda_0\) gives
\[\Lambda_0 = T_0\Biggl[1 - 1 + \frac{1}{r}\int_0^r \mu(s)\,ds\Biggr] = \frac{1}{r} T_0 \int_0^r \mu(s)\,ds.\]Finally, for the linear \(K(s)\) assumed above,
\[\mu(s) = \exp\Bigl(\int_0^s K(u)du\Bigr) = \exp\Bigl(K_0 s + \frac{K_1-K_0}{2r} s^2\Bigr),\]so we obtain exactly
\[\Lambda_0 = \frac{1}{r} T_0 \int_0^r \exp\Bigl(K_0 s + \frac{K_1-K_0}{2r} s^2\Bigr)\,ds,\]which is the expression used above (with index \(i\)):
\[\Lambda_i = \frac{1}{r_i} T_i \int_0^{r_i} \exp\left(K_i s + \frac{K_{i+1}-K_i}{2 r_i} s^2\right) ds.\]In other words, \(\Lambda_i\) is precisely the coefficient multiplying \((J_i - J_{i+1})\) that comes from solving the ODE with both \(K\) and \(J\) linear in \(s\).
3.3. Stable scalar evaluation and polarized approximation.
For scalar \(K_0\) and \(K_1\), introduce
\[\begin{split}\tau &= \frac{r}{2}(K_0+K_1),\\ \delta &= \frac{r}{2}(K_1-K_0),\\ q &= 1-\frac{s}{r}.\end{split}\]Combining \(T_0=\exp(-\tau)\) with the integral above gives
\[\Lambda_0(\tau,\delta) = \int_0^1 \exp\left[-\tau q-\delta q(1-q)\right]\,dq.\]This representation remains regular when the propagation gradient vanishes. It also gives the two partial derivatives directly:
\[\begin{split}\frac{\partial\Lambda_0}{\partial\tau} &= -\int_0^1 q \exp\left[-\tau q-\delta q(1-q)\right]\,dq,\\ \frac{\partial\Lambda_0}{\partial\delta} &= -\int_0^1 q(1-q) \exp\left[-\tau q-\delta q(1-q)\right]\,dq.\end{split}\]ARTS evaluates these three quantities together. An expansion in \(\delta\) uses stable exponential moments near equal endpoints. Away from that regime, the increasing-gradient branch uses a scalar Dawson representation and the decreasing-gradient branch uses the scaled complementary error function. A Watson endpoint expansion avoids cancellation for sufficiently thick layers. The derivatives with respect to \(K_0\), \(K_1\), and \(r\) then follow analytically by the chain rule.
ARTS does not apply this scalar construction component by component to a polarized matrix. Instead, polarized
linpropuses \(\Lambda_{\rm pol}=\phi_1(G_{\rm av})Q+c\mathbf{I}\) with ordinary endpoint-average transmission, as defined above. Themagop_linsrcoption uses the same augmented-source factor \(Q\) but replaces \(G_{\rm av}\) by the ordered generator \(-\Omega\) in both the source operator and transmission.