Formal solution and Schwarzschild equation

The radiation transport equation (RTE) for time-independent plane-parallel atmosphere has well known form (e.g. Rutten, 2003):
\begin{equation}
\mu\,\frac{\mathrm{d}I_{\nu}}{\mathrm{d} \tau_{\nu}} = I_{\nu} - S_{\nu},
\end{equation}
 where $I_{\nu}$ is specific intensity, $S_{\nu}$ source function, $\tau_{\nu}$ is the radial optical depth and $\mu$ is the cosine of the angle between the line of sight and vertical axis. Mathematically, this is a linear first order differential equation with constant coefficients. Physically, it describes propagation of radiation through an interacting medium. Properties of the medium are hidden in the source function and the optical depth. However, to make the first step toward the solution of RTE, we will assume that the optical depth is independent variable and that the source function is a know function of it.


Formal solution

The formal solution of the equation (1) in the interval $[\tau_1, \tau_2]$ is ($\nu$ index is omitted for simplicity):
\begin{equation}
I(\tau_1, \mu) = I(\tau_2, \mu)\; \mathrm{e}^{-(\tau_2 - \tau_1)/\mu} + \int_{\tau_1}^{\tau_2}S(t)\;\mathrm{e}^{-(t - \tau_1)/\mu}\,\frac{\mathrm{d}t}{\mu}
\end{equation}
This equation says that the radiation intensity outgoing from the slab is equal to the incoming intensity diminished by the absorption plus the contribution from the slab. Particularly important solution is for the semi-infinite atmosphere ($\tau_1 = 0, \tau_2 \rightarrow \infty$):
\begin{equation}
I(0, \mu) = \int_{0}^{\infty}S(t)\;\mathrm{e}^{-t/\mu}\,\frac{\mathrm{d}t}{\mu}
\end{equation}
Example: If the source function is linear with depth, $S = a + b\tau$, we immediately get the well-known Eddington-Barbier relation for the emergent intensity at the surface: $I = a + b\mu$.

For an arbitrary point along the optical depth scale, we can separate the formal solution in two parts: for outgoing radiation $(\mu \ge 0)$ and for incoming part $(\mu \le 0)$. For the outgoing radiation $(\tau1 = \tau, \tau_2 \rightarrow \infty)$, the formal solution becomes:
\begin{equation*}
I(\tau, \mu) = \lim_{\tau_2 \rightarrow \infty}I(\tau_2, \mu)\; \mathrm{e}^{-(\tau_2 - \tau)/\mu} + \int_{\tau}^{\infty}S(t)\;\mathrm{e}^{-(t - \tau)/\mu}\,\frac{\mathrm{d}t}{\mu}.
\end{equation*}
The first term on the right-hand side is zero (a common boundary condition (see Mihalas, 19.., p..), so for the outgoing radiation we get
\begin{equation}
I^+(\tau, \mu) = \int_{\tau}^{\infty}S(t)\;\mathrm{e}^{-(t - \tau)/\mu}\,\frac{\mathrm{d}t}{\mu}\;\;\;\;(1 \ge \mu \ge 0).
\end{equation}
For incoming radiation at the optical depth $t$, in a similar manner, we get:
\begin{equation}
I^-(\tau, \mu) = - \int_{0}^{\tau}S(t)\;\mathrm{e}^{-(t - \tau)/\mu}\,\frac{\mathrm{d}t}{\mu}\;\;\;\;(0 \ge \mu \ge -1).
\end{equation}
Equations (4) and (5) represent the formal solution of RTE in a semi-infinite plane-parallel atmosphere.

Schwarzschild equation

The next step is to express the mean intensity $J$ in terms of known source function, so we integrate the specific intensity over the directions:
\begin{eqnarray*}
J(\tau) &=& \frac{1}{2}\,\int_{-1}^{1} I(\tau, \mu)\, \mathrm{d}\mu = \\
&=& \frac{1}{2}\,\int_{0}^{1}\mathrm{d}\mu \,\int_{\tau}^{\infty}S(t)\;\mathrm{e}^{-(t - \tau)/\mu}\,\frac{\mathrm{d}t}{\mu} - \frac{1}{2}\,\int_{-1}^{0}\mathrm{d}\mu \int_{0}^{\tau}S(t)\;\mathrm{e}^{-(t - \tau)/\mu}\,\frac{\mathrm{d}t}{\mu}.
\end{eqnarray*}
To rewrite this in a more concise form, we introduce substitution
$$w = \frac{1}{\mu} \;\;\;\;\;\;\; \frac{\mathrm{d}w}{w} = -\frac{\mathrm{d}\mu}{\mu}$$
in the first integral, and 
$$w = -\frac{1}{\mu} \;\;\;\;\;\;\; \frac{\mathrm{d}w}{w} = -\frac{\mathrm{d}\mu}{\mu}$$
in the second, to get:
\begin{equation}
J(\tau) = \int_{\tau}^{\infty}S(t)\,\mathrm{d}t\;\int_{1}^{\infty}\mathrm{e}^{-w(t+\tau)}\frac{\mathrm{d}w}{w}
+ \frac{1}{2}\int_{0}^{\tau}S(t)\,\mathrm{d}t\;\int_{1}^{\infty}\mathrm{e}^{-w(t+\tau)}\frac{\mathrm{d}w}{w}.
\end{equation}
Further compression of this expression can be made if we introduce exponential integrals. The first exponential integral is:
\begin{equation}
\mathrm{E}_1(x) = \int_1^{\infty}\frac{\mathrm{e}^{-xt}}{t}\,\mathrm{d}t,
\end{equation}
thus the mean intensity becomes:
\begin{equation}
J(\tau) = \frac{1}{2}\int_{\tau}^{\infty}S(\tau) \mathrm{E}_1(t-\tau)\,\mathrm{d}t + \frac{1}{2}\int_{0}^{\tau}S(\tau) \mathrm{E}_1(t-\tau)\,\mathrm{d}t.
\end{equation}
or equivalently:
\begin{equation}
J(\tau) = \frac{1}{2}\int_{0}^{\infty}S(\tau) \mathrm{E}_1|t-\tau|\,\mathrm{d}t.
\end{equation}
The last equation is known as Schwarzschild equation or Schwarzschild-Milne equation.

Lambda operator

Let's introduce lambda operator as:
\begin{equation}
\Lambda[f(t)] \equiv \frac{1}{2}\int_{0}^{\infty}f(t) \mathrm{E}_1|t-\tau|\,
\mathrm{d}t.
\end{equation}
Schwarzschild equation can now be rewritten as
\begin{equation}
J(\tau) = \Lambda[S(t)].
\end{equation}
This form of the equation is the starting point for many numerical methods for solving the RTE. 

In the forthcoming entries I'll present some properties of the lambda operator and demonstrate how to use in real computations.